summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/base/ode.asy
diff options
context:
space:
mode:
authorKarl Berry <karl@freefriends.org>2009-06-21 23:42:52 +0000
committerKarl Berry <karl@freefriends.org>2009-06-21 23:42:52 +0000
commit8542b3da905deccc50b1bd6d21c7474a29286b3f (patch)
treec9a21ed9311ad8f967c67d17efad475eef50c5c1 /Build/source/utils/asymptote/base/ode.asy
parent95be153f845a361eb22745faa0d1ce509eb18c39 (diff)
asymptote 1.77 sources
git-svn-id: svn://tug.org/texlive/trunk@13865 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/utils/asymptote/base/ode.asy')
-rw-r--r--Build/source/utils/asymptote/base/ode.asy257
1 files changed, 251 insertions, 6 deletions
diff --git a/Build/source/utils/asymptote/base/ode.asy b/Build/source/utils/asymptote/base/ode.asy
index 5135dd9443f..643c56feec7 100644
--- a/Build/source/utils/asymptote/base/ode.asy
+++ b/Build/source/utils/asymptote/base/ode.asy
@@ -1,15 +1,260 @@
-real euler(real y, real f(real x, real y), real a, real b=a, int n=0,
- real h=0, bool dynamic=false, real tolmin=0, real tolmax=0)
+real stepfactor=2.0; // Maximum dynamic step size adjustment factor.
+
+struct RKTableau
+{
+ int order;
+ real[] steps;
+ real[][] weights;
+ real[] highOrderWeights;
+ real[] lowOrderWeights;
+ real pgrow;
+ real pshrink;
+
+ void operator init(int order, real[][] weights, real[] highOrderWeights,
+ real[] lowOrderWeights=new real[],
+ real[] steps=sequence(new real(int i) {
+ return sum(weights[i]);},weights.length)) {
+ this.order=order;
+ this.steps=steps;
+ this.weights=weights;
+ this.highOrderWeights=highOrderWeights;
+ this.lowOrderWeights=lowOrderWeights;
+ pgrow=(order > 0) ? 1/order : 0;
+ pshrink=(order > 1) ? 1/(order-1) : pgrow;
+ }
+}
+
+// First-Order Euler
+RKTableau Euler=RKTableau(1,new real[][],
+ new real[] {1});
+
+// Second-Order Runge-Kutta
+RKTableau RK2=RKTableau(2,new real[][] {{1/2}},
+ new real[] {0,1});
+
+// Second-Order Predictor-Corrector
+RKTableau PC=RKTableau(2,new real[][] {{1}},
+ new real[] {1/2,1/2});
+
+// Third-Order Classical Runge-Kutta
+RKTableau RK3=RKTableau(3,new real[][] {{1/2},{-1,2}},
+ new real[] {1/6,2/3,1/6});
+
+// Third-Order Bogacki-Shampine Runge-Kutta
+RKTableau RK3BS=RKTableau(3,new real[][] {{1/2},{0,3/4}},
+ new real[] {2/9,1/3,4/9}, // 3rd order
+ new real[] {7/24,1/4,1/3,1/8}); // 2nd order
+
+// Fourth-Order Classical Runge-Kutta
+RKTableau RK4=RKTableau(4,new real[][] {{1/2},{0,1/2},{0,0,1}},
+ new real[] {1/6,1/3,1/3,1/6});
+
+// Fifth-Order Cash-Karp Runge-Kutta
+RKTableau RK5CK=RKTableau(5,new real[][] {{1/5},
+ {3/40,9/40},
+ {3/10,-9/10,6/5},
+ {-11/54,5/2,-70/27,35/27},
+ {1631/55296,175/512,575/13824,
+ 44275/110592,253/4096}},
+ new real[] {37/378,0,250/621,125/594,
+ 0,512/1771}, // 5th order
+ new real[] {2825/27648,0,18575/48384,13525/55296,
+ 277/14336,1/4}); // 4th order
+
+// Fifth-Order Fehlberg Runge-Kutta
+RKTableau RK5F=RKTableau(5,new real[][] {{1/4},
+ {3/32,9/32},
+ {1932/2197,-7200/2197,7296/2197},
+ {439/216,-8,3680/513,-845/4104},
+ {-8/27,2,-3544/2565,1859/4104,
+ -11/40}},
+ new real[] {16/135,0,6656/12825,28561/56430,-9/50,2/55}, // 5th order
+ new real[] {25/216,0,1408/2565,2197/4104,-1/5,0}); // 4th order
+
+// Fifth-Order Dormand-Prince Runge-Kutta
+RKTableau RK5DP=RKTableau(5,new real[][] {{1/5},
+ {3/40,9/40},
+ {44/45,-56/15,32/9},
+ {19372/6561,-25360/2187,64448/6561,
+ -212/729},
+ {9017/3168,-355/33,46732/5247,49/176,
+ -5103/18656}},
+ new real[] {35/384,0,500/1113,125/192,-2187/6784,
+ 11/84}, // 5th order
+ new real[] {5179/57600,0,7571/16695,393/640,
+ -92097/339200,187/2100,1/40}); // 4th order
+
+real error(real error, real initial, real norm, real lowOrder, real diff)
+{
+ if(initial != 0.0 && lowOrder != initial) {
+ static real epsilon=realMin/realEpsilon;
+ real denom=max(abs(norm),abs(initial))+epsilon;
+ return max(error,max(abs(diff)/denom));
+ }
+ return error;
+}
+
+real adjust(real h, real error, real t, real tolmin, real tolmax,
+ real dtmin, real dtmax, RKTableau tableau, bool verbose=true)
+{
+ real dt=h;
+ void report(real t) {
+ if(h != dt)
+ write("Time step changed from "+(string) dt+" to "+(string) h+" at t="+
+ (string) t+".");
+ }
+ if(error > tolmax) {
+ h=max(h*max((tolmin/error)^tableau.pshrink,1/stepfactor),dtmin);
+ if(verbose) report(t);
+ return h;
+ }
+ if(error > 0 && error < tolmin) {
+ h=min(h*min((tolmin/error)^tableau.pgrow,stepfactor),dtmax);
+ if(verbose) report(t+dt);
+ }
+ return h;
+}
+
+// Integrate dy/dt=f(t,y) from a to b using initial conditions y,
+// specifying either the step size h or the number of steps n.
+real integrate(real y, real f(real t, real y), real a, real b=a, real h=0,
+ int n=0, bool dynamic=false, real tolmin=0, real tolmax=0,
+ real dtmin=0, real dtmax=realMax, RKTableau tableau,
+ bool verbose=false)
+{
+ if(h == 0) {
+ if(b == a) return y;
+ if(n == 0) abort("Either n or h must be specified");
+ else h=(b-a)/n;
+ }
+ real t=a;
+ real f0;
+ bool fsal=tableau.lowOrderWeights.length > tableau.highOrderWeights.length;
+ if(fsal) f0=f(t,y);
+ if(tableau.lowOrderWeights.length == 0) dynamic=false;
+
+ while(t < b) {
+ real[] predictions={fsal ? f0 : f(t,y)};
+ for(int i=0; i < tableau.steps.length; ++i)
+ predictions.push(f(t+h*tableau.steps[i],
+ y+h*dot(tableau.weights[i],predictions)));
+
+ real highOrder=h*dot(tableau.highOrderWeights,predictions);
+ if(dynamic) {
+ real f1;
+ if(fsal) {
+ f1=f(t+h,y+highOrder);
+ predictions.push(f1);
+ }
+ real lowOrder=h*dot(tableau.lowOrderWeights,predictions);
+ real error;
+ error=error(error,y,y+highOrder,y+lowOrder,highOrder-lowOrder);
+ real dt=h;
+ h=adjust(h,error,t,tolmin,tolmax,dtmin,min(dtmax,b-t-h),tableau,verbose);
+ if(h >= dt) {
+ t += dt;
+ y += highOrder;
+ f0=f1;
+ }
+ } else {
+ t += h;
+ y += highOrder;
+ }
+ h=min(h,b-t);
+ if(t >= b || t+h == t) break;
+ }
+ return y;
+}
+
+// Integrate a set of equations, dy/dt=f(t,y), from a to b using initial
+// conditions y, specifying either the step size h or the number of steps n.
+real[] integrate(real[] y, real[] f(real t, real[] y), real a, real b=a,
+ real h=0, int n=0, bool dynamic=false,
+ real tolmin=0, real tolmax=0, real dtmin=0, real dtmax=realMax,
+ RKTableau tableau, bool verbose=false)
{
if(h == 0) {
if(b == a) return y;
if(n == 0) abort("Either n or h must be specified");
else h=(b-a)/n;
}
- real x=a;
- for(int i=0; i < n; ++i) {
- y += h*f(x,y);
- x += h;
+ real[] y=copy(y);
+ real t=a;
+ real[] f0;
+ bool fsal=tableau.lowOrderWeights.length > tableau.highOrderWeights.length;
+ if(fsal) f0=f(t,y);
+ if(tableau.lowOrderWeights.length == 0) dynamic=false;
+
+ while(t < b) {
+ real[][] predictions={fsal ? f0 : f(t,y)};
+ for(int i=0; i < tableau.steps.length; ++i)
+ predictions.push(f(t+h*tableau.steps[i],
+ y+h*tableau.weights[i]*predictions));
+
+ real[] highOrder=h*tableau.highOrderWeights*predictions;
+ if(dynamic) {
+ real[] f1;
+ if(fsal) {
+ f1=f(t+h,y+highOrder);
+ predictions.push(f1);
+ }
+ real[] lowOrder=h*tableau.lowOrderWeights*predictions;
+ real error;
+ for(int i=0; i < y.length; ++i)
+ error=error(error,y[i],y[i]+highOrder[i],y[i]+lowOrder[i],
+ highOrder[i]-lowOrder[i]);
+ real dt=h;
+ h=adjust(h,error,t,tolmin,tolmax,dtmin,min(dtmax,b-t-h),tableau,verbose);
+ if(h >= dt) {
+ t += dt;
+ y += highOrder;
+ f0=f1;
+ }
+ } else {
+ t += h;
+ y += highOrder;
+ }
+ h=min(h,b-t);
+ if(t >= b || t+h == t) break;
}
return y;
}
+
+real[][] finiteDifferenceJacobian(real[] f(real[]), real[] t,
+ real[] h=sqrtEpsilon*abs(t))
+{
+ real[] ft=f(t);
+ real[][] J=new real[t.length][ft.length];
+ real[] ti=copy(t);
+ real tlast=ti[0];
+ ti[0] += h[0];
+ J[0]=(f(ti)-ft)/h[0];
+ for(int i=1; i < t.length; ++i) {
+ ti[i-1]=tlast;
+ tlast=ti[i];
+ ti[i] += h[i];
+ J[i]=(f(ti)-ft)/h[i];
+ }
+ return transpose(J);
+}
+
+// Solve simultaneous nonlinear system by Newton's method.
+real[] newton(int iterations=100, real[] f(real[]), real[][] jacobian(real[]),
+ real[] t)
+{
+ real[] t=copy(t);
+ for(int i=0; i < iterations; ++i)
+ t += solve(jacobian(t),-f(t));
+ return t;
+}
+
+real[] solveBVP(real[] f(real, real[]), real a, real b=a, real h=0, int n=0,
+ real[] initial(real[]), real[] discrepancy(real[]),
+ real[] guess, RKTableau tableau, int iterations=100)
+{
+ real[] g(real[] t) {
+ return discrepancy(integrate(initial(t),f,a,b,h,n,tableau));
+ }
+ real[][] jacobian(real[] t) {return finiteDifferenceJacobian(g,t);}
+ return initial(newton(iterations,g,jacobian,guess));
+}