diff options
author | Karl Berry <karl@freefriends.org> | 2009-07-03 00:28:16 +0000 |
---|---|---|
committer | Karl Berry <karl@freefriends.org> | 2009-07-03 00:28:16 +0000 |
commit | fc2b1bcc12c2d8f670d19f9d8a2baef03e9cd63a (patch) | |
tree | 35138f2dd9b288a4ca15c4772861766f826434fe /Master/texmf/asymptote/ode.asy | |
parent | 47058e8e026a4d06115016a53af65d0d8921d606 (diff) |
asymptote 1.79
git-svn-id: svn://tug.org/texlive/trunk@14080 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf/asymptote/ode.asy')
-rw-r--r-- | Master/texmf/asymptote/ode.asy | 314 |
1 files changed, 237 insertions, 77 deletions
diff --git a/Master/texmf/asymptote/ode.asy b/Master/texmf/asymptote/ode.asy index 643c56feec7..77ab8d6e589 100644 --- a/Master/texmf/asymptote/ode.asy +++ b/Master/texmf/asymptote/ode.asy @@ -1,40 +1,151 @@ -real stepfactor=2.0; // Maximum dynamic step size adjustment factor. +real stepfactor=2; // 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; + 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.steps=steps; + a.factors=array(a.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; } } +real[] Coeff={1,1/2,1/6,1/24,1/120,1/720,1/5040,1/40320,1/362880,1/3628800, + 1/39916800.0,1/479001600.0,1/6227020800.0,1/87178291200.0, + 1/1307674368000.0,1/20922789888000.0,1/355687428096000.0, + 1/6402373705728000.0,1/121645100408832000.0, + 1/2432902008176640000.0,1/51090942171709440000.0, + 1/1124000727777607680000.0}; + +real phi1(real x) {return x != 0 ? expm1(x)/x : 1;} + +real phi2(real x) +{ + real x2=x*x; + if(fabs(x) > 1) return (exp(x)-x-1)/x2; + real x3=x2*x; + real x5=x2*x3; + if(fabs(x) < 0.1) + return Coeff[1]+x*Coeff[2]+x2*Coeff[3]+x3*Coeff[4]+x2*x2*Coeff[5] + +x5*Coeff[6]+x3*x3*Coeff[7]+x5*x2*Coeff[8]+x5*x3*Coeff[9]; + else { + real x7=x5*x2; + real x8=x7*x; + return Coeff[1]+x*Coeff[2]+x2*Coeff[3]+x3*Coeff[4]+x2*x2*Coeff[5] + +x5*Coeff[6]+x3*x3*Coeff[7]+x7*Coeff[8]+x8*Coeff[9] + +x8*x*Coeff[10]+x5*x5*Coeff[11]+x8*x3*Coeff[12]+x7*x5*Coeff[13]+ + x8*x5*Coeff[14]+x7*x7*Coeff[15]+x8*x7*Coeff[16]+x8*x8*Coeff[17]; + } +} + +real phi3(real x) +{ + real x2=x*x; + real x3=x2*x; + if(fabs(x) > 1.6) return (exp(x)-0.5*x2-x-1)/x3; + real x5=x2*x3; + if(fabs(x) < 0.1) + return Coeff[2]+x*Coeff[3]+x2*Coeff[4]+x3*Coeff[5] + +x2*x2*Coeff[6]+x5*Coeff[7]+x3*x3*Coeff[8]+x5*x2*Coeff[9] + +x5*x3*Coeff[10]; + else { + real x7=x5*x2; + real x8=x7*x; + real x16=x8*x8; + return Coeff[2]+x*Coeff[3]+x2*Coeff[4]+x3*Coeff[5] + +x2*x2*Coeff[6]+x5*Coeff[7]+x3*x3*Coeff[8]+x5*x2*Coeff[9] + +x5*x3*Coeff[10]+x8*x*Coeff[11] + +x5*x5*Coeff[12]+x8*x3*Coeff[13]+x7*x5*Coeff[14] + +x8*x5*Coeff[15]+x7*x7*Coeff[16]+x8*x7*Coeff[17]+x16*Coeff[18] + +x16*x*Coeff[19]+x16*x2*Coeff[20]; + } +} + +void expfactors(real x, coefficients a) +{ + for(int i=0; i < a.steps.length; ++i) + a.factors[i]=exp(x*a.steps[i]); + a.factors[a.steps.length]=exp(x); +} + // First-Order Euler RKTableau Euler=RKTableau(1,new real[][], new real[] {1}); +// First-Order Exponential Euler +RKTableau E_Euler=RKTableau(1,new real[][], new real[] {1}, + new void(real h, real c, coefficients a) { + real x=-c*h; + expfactors(x,a); + a.highOrderWeights[0]=phi1(x); + }); + // Second-Order Runge-Kutta RKTableau RK2=RKTableau(2,new real[][] {{1/2}}, - new real[] {0,1}); + new real[] {0,1}, // 2nd order + new real[] {1,0}); // 1st order + +// Second-Order Exponential Runge-Kutta +RKTableau E_RK2=RKTableau(2,new real[][] {{1/2}}, + new real[] {0,1}, // 2nd order + new real[] {1,0}, // 1st order + new void(real h, real c, coefficients a) { + real x=-c*h; + expfactors(x,a); + a.weights[0][0]=1/2*phi1(x/2); + real w=phi1(x); + a.highOrderWeights[0]=0; + a.highOrderWeights[1]=w; + a.lowOrderWeights[0]=w; + }); // Second-Order Predictor-Corrector RKTableau PC=RKTableau(2,new real[][] {{1}}, - new real[] {1/2,1/2}); + new real[] {1/2,1/2}, // 2nd order + new real[] {1,0}); // 1st order + +// Second-Order Exponential Predictor-Corrector +RKTableau E_PC=RKTableau(2,new real[][] {{1}}, + new real[] {1/2,1/2}, // 2nd order + new real[] {1,0}, // 1st order + new void(real h, real c, coefficients a) { + real x=-c*h; + expfactors(x,a); + real w=phi1(x); + a.weights[0][0]=w; + a.highOrderWeights[0]=w/2; + a.highOrderWeights[1]=w/2; + a.lowOrderWeights[0]=w; + }); // Third-Order Classical Runge-Kutta RKTableau RK3=RKTableau(3,new real[][] {{1/2},{-1,2}}, @@ -45,17 +156,41 @@ 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 +// Third-Order Exponential Bogacki-Shampine Runge-Kutta +RKTableau E_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 + new void(real h, real c, coefficients a) { + real x=-c*h; + expfactors(x,a); + real w=phi1(x); + real w2=phi2(x); + a.weights[0][0]=1/2*phi1(x/2); + real a11=9/8*phi2(3/4*x)+3/8*phi2(x/2); + a.weights[1][0]=3/4*phi1(3/4*x)-a11; + a.weights[1][1]=a11; + real a21=1/3*w; + real a22=4/3*w2-2/9*w; + a.highOrderWeights[0]=w-a21-a22; + a.highOrderWeights[1]=a21; + a.highOrderWeights[2]=a22; + a.lowOrderWeights[0]=w-17/12*w2; + a.lowOrderWeights[1]=w2/2; + a.lowOrderWeights[2]=2/3*w2; + a.lowOrderWeights[3]=w2/4; + }); + // 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}}, +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,9 +219,9 @@ 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) { + if(initial != 0 && lowOrder != initial) { static real epsilon=realMin/realEpsilon; real denom=max(abs(norm),abs(initial))+epsilon; return max(error,max(abs(diff)/denom)); @@ -94,130 +229,154 @@ 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 dtmin, real dtmax, RKTableau tableau, bool verbose=true) +void report(real old, real h, real t) { - 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); - } + write("Time step changed from "+(string) old+" to "+(string) h+" at t="+ + (string) t+"."); +} + +real adjust(real h, real error, real tolmin, real tolmax, RKTableau tableau) +{ + if(error > tolmax) + h *= max((tolmin/error)^tableau.pshrink,1/stepfactor); + else if(error > 0 && error < tolmin) + h *= min((tolmin/error)^tableau.pgrow,stepfactor); 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[] Y={y}; + if(h == 0) { - if(b == a) return y; + if(b == a) return Y; if(n == 0) abort("Either n or h must be specified"); else h=(b-a)/n; } + + real f(real t, real y)=(c == 0 || tableau.exponential) ? g : + new real(real t, real y) {return g(t,y)-c*y;}; + + 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; + real dt=h; while(t < b) { + h=min(h,b-t); + if(t+h == t) break; + if(h != dt) { + if(verbose) report(dt,h,t); + tableau.stepDependence(h,c,tableau.a); + dt=h; + } + 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))); + for(int i=0; i < tableau.a.steps.length; ++i) + predictions.push(f(t+h*tableau.a.steps[i], + 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 y0=tableau.a.factors[tableau.a.steps.length]*y; if(dynamic) { real f1; if(fsal) { - f1=f(t+h,y+highOrder); + f1=f(t+h,y0+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); - real dt=h; - h=adjust(h,error,t,tolmin,tolmax,dtmin,min(dtmax,b-t-h),tableau,verbose); + error=error(error,y,y0+lowOrder,y0+highOrder,highOrder-lowOrder); + h=adjust(h,error,tolmin,tolmax,tableau); if(h >= dt) { t += dt; - y += highOrder; + y=y0+highOrder; + Y.push(y); f0=f1; } + h=min(max(h,dtmin),dtmax); } else { t += h; - y += highOrder; + y=y0+highOrder; + Y.push(y); } - h=min(h,b-t); - if(t >= b || t+h == t) break; } - return y; + 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) +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[][] Y={copy(y)}; + if(h == 0) { - if(b == a) return y; + if(b == a) return Y; if(n == 0) abort("Either n or h must be specified"); else h=(b-a)/n; } - 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; + real dt=h; while(t < b) { + h=min(h,b-t); + if(t+h == t) break; + if(h != dt) { + if(verbose) report(dt,h,t); + dt=h; + } + 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)); + for(int i=0; i < tableau.a.steps.length; ++i) + predictions.push(f(t+h*tableau.a.steps[i], + 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,tolmin,tolmax,tableau); if(h >= dt) { t += dt; y += highOrder; + Y.push(y); f0=f1; } + h=min(max(h,dtmin),dtmax); } else { t += h; y += highOrder; + Y.push(y); } - h=min(h,b-t); - if(t >= b || t+h == t) break; } - return y; + return Y; } real[][] finiteDifferenceJacobian(real[] f(real[]), real[] t, @@ -253,7 +412,8 @@ real[] solveBVP(real[] f(real, real[]), real a, real b=a, real h=0, int n=0, real[] guess, RKTableau tableau, int iterations=100) { real[] g(real[] t) { - return discrepancy(integrate(initial(t),f,a,b,h,n,tableau)); + real[][] y=integrate(initial(t),f,a,b,h,n,tableau); + return discrepancy(y[y.length-1]); } real[][] jacobian(real[] t) {return finiteDifferenceJacobian(g,t);} return initial(newton(iterations,g,jacobian,guess)); |