summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/base/ode.asy
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/utils/asymptote/base/ode.asy')
-rw-r--r--Build/source/utils/asymptote/base/ode.asy125
1 files changed, 82 insertions, 43 deletions
diff --git a/Build/source/utils/asymptote/base/ode.asy b/Build/source/utils/asymptote/base/ode.asy
index 643c56feec7..fe65f2f1c8e 100644
--- a/Build/source/utils/asymptote/base/ode.asy
+++ b/Build/source/utils/asymptote/base/ode.asy
@@ -1,24 +1,39 @@
real stepfactor=2.0; // Maximum dynamic step size adjustment factor.
-struct RKTableau
+struct coefficients
{
- int order;
- real[] steps;
+ real[] factors;
real[][] weights;
real[] highOrderWeights;
real[] lowOrderWeights;
+}
+
+struct RKTableau
+{
+ int order;
+ real[] steps;
+ coefficients a;
+ void stepDependence(real h, real c, coefficients a) {}
+
real pgrow;
real pshrink;
+ bool exponential;
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)) {
+ return sum(weights[i]);},weights.length),
+ void stepDependence(real, real, coefficients)=null) {
this.order=order;
this.steps=steps;
- this.weights=weights;
- this.highOrderWeights=highOrderWeights;
- this.lowOrderWeights=lowOrderWeights;
+ a.factors=array(steps.length+1,1);
+ a.weights=weights;
+ a.highOrderWeights=highOrderWeights;
+ a.lowOrderWeights=lowOrderWeights;
+ if(stepDependence != null) {
+ this.stepDependence=stepDependence;
+ exponential=true;
+ }
pgrow=(order > 0) ? 1/order : 0;
pshrink=(order > 1) ? 1/(order-1) : pgrow;
}
@@ -28,6 +43,13 @@ struct RKTableau
RKTableau Euler=RKTableau(1,new real[][],
new real[] {1});
+RKTableau E_Euler=RKTableau(1,new real[][], new real[] {1},
+ new void (real h, real c, coefficients a) {
+ real x=-c*h;
+ a.factors[0]=exp(x);
+ a.highOrderWeights[0]=x != 0 ? expm1(x)/x : 1;
+ });
+
// Second-Order Runge-Kutta
RKTableau RK2=RKTableau(2,new real[][] {{1/2}},
new real[] {0,1});
@@ -50,12 +72,12 @@ 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}},
+RKTableau RK5=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,
@@ -84,7 +106,7 @@ RKTableau RK5DP=RKTableau(5,new real[][] {{1/5},
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)
+real error(real error, real initial, real lowOrder, real norm, real diff)
{
if(initial != 0.0 && lowOrder != initial) {
static real epsilon=realMin/realEpsilon;
@@ -94,73 +116,89 @@ real error(real error, real initial, real norm, real lowOrder, real diff)
return error;
}
-real adjust(real h, real error, real t, real tolmin, real tolmax,
+real adjust(real h, real error, real t, real c, 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(h != dt) {
+ tableau.stepDependence(h,c,tableau.a);
+ if(verbose)
+ 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);
+ report(t);
return h;
}
if(error > 0 && error < tolmin) {
h=min(h*min((tolmin/error)^tableau.pgrow,stepfactor),dtmax);
- if(verbose) report(t+dt);
+ report(t+dt);
}
return h;
}
-// Integrate dy/dt=f(t,y) from a to b using initial conditions y,
+// Integrate dy/dt+cy=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)
+real integrate(real y, real c=0, real g(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)
{
+ real f(real t, real y)=(c == 0 || tableau.exponential) ? g :
+ new real(real t, real y) {return g(t,y)-c*y;};
+
if(h == 0) {
if(b == a) return y;
if(n == 0) abort("Either n or h must be specified");
else h=(b-a)/n;
}
+
+ tableau.stepDependence(h,c,tableau.a);
+
real t=a;
real f0;
- bool fsal=tableau.lowOrderWeights.length > tableau.highOrderWeights.length;
+ if(tableau.a.lowOrderWeights.length == 0) dynamic=false;
+ bool fsal=dynamic &&
+ (tableau.a.lowOrderWeights.length > tableau.a.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)));
+ tableau.a.factors[i]*y+h*dot(tableau.a.weights[i],
+ predictions)));
- real highOrder=h*dot(tableau.highOrderWeights,predictions);
+ real highOrder=h*dot(tableau.a.highOrderWeights,predictions);
+ real Y=tableau.a.factors[tableau.steps.length]*y;
if(dynamic) {
real f1;
if(fsal) {
- f1=f(t+h,y+highOrder);
+ f1=f(t+h,Y+highOrder);
predictions.push(f1);
}
- real lowOrder=h*dot(tableau.lowOrderWeights,predictions);
+ real lowOrder=h*dot(tableau.a.lowOrderWeights,predictions);
real error;
- error=error(error,y,y+highOrder,y+lowOrder,highOrder-lowOrder);
+ error=error(error,y,Y+lowOrder,Y+highOrder,highOrder-lowOrder);
real dt=h;
- h=adjust(h,error,t,tolmin,tolmax,dtmin,min(dtmax,b-t-h),tableau,verbose);
+ h=adjust(h,error,t,c,tolmin,tolmax,dtmin,min(dtmax,b-t-h),tableau,verbose);
if(h >= dt) {
t += dt;
- y += highOrder;
+ y=Y+highOrder;
f0=f1;
}
} else {
t += h;
- y += highOrder;
+ y=Y+highOrder;
+ }
+ real remain=b-t;
+ if(h > remain) {
+ h=remain;
+ tableau.stepDependence(h,c,tableau.a);
}
- h=min(h,b-t);
if(t >= b || t+h == t) break;
}
return y;
@@ -181,30 +219,31 @@ real[] integrate(real[] y, real[] f(real t, real[] y), real a, real b=a,
real[] y=copy(y);
real t=a;
real[] f0;
- bool fsal=tableau.lowOrderWeights.length > tableau.highOrderWeights.length;
+ if(tableau.a.lowOrderWeights.length == 0) dynamic=false;
+ bool fsal=dynamic &&
+ (tableau.a.lowOrderWeights.length > tableau.a.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));
+ y+h*tableau.a.weights[i]*predictions));
- real[] highOrder=h*tableau.highOrderWeights*predictions;
+ real[] highOrder=h*tableau.a.highOrderWeights*predictions;
if(dynamic) {
real[] f1;
if(fsal) {
f1=f(t+h,y+highOrder);
predictions.push(f1);
}
- real[] lowOrder=h*tableau.lowOrderWeights*predictions;
+ real[] lowOrder=h*tableau.a.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],
+ error=error(error,y[i],y[i]+lowOrder[i],y[i]+highOrder[i],
highOrder[i]-lowOrder[i]);
real dt=h;
- h=adjust(h,error,t,tolmin,tolmax,dtmin,min(dtmax,b-t-h),tableau,verbose);
+ h=adjust(h,error,t,0,tolmin,tolmax,dtmin,min(dtmax,b-t-h),tableau,verbose);
if(h >= dt) {
t += dt;
y += highOrder;