diff options
Diffstat (limited to 'Master/texmf/asymptote/ode.asy')
-rw-r--r-- | Master/texmf/asymptote/ode.asy | 457 |
1 files changed, 0 insertions, 457 deletions
diff --git a/Master/texmf/asymptote/ode.asy b/Master/texmf/asymptote/ode.asy deleted file mode 100644 index 0b88060e081..00000000000 --- a/Master/texmf/asymptote/ode.asy +++ /dev/null @@ -1,457 +0,0 @@ -real stepfactor=2; // Maximum dynamic step size adjustment factor. - -struct coefficients -{ - 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), - void stepDependence(real, real, coefficients)=null) { - this.order=order; - 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}, // 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}, // 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}}, - 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 - -// 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 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, - 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 lowOrder, real norm, real diff) -{ - 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)); - } - return error; -} - -void report(real old, real h, real t) -{ - 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; -} - -struct solution -{ - real[] t; - real[] y; -} - -void write(solution S) -{ - for(int i=0; i < S.t.length; ++i) - write(S.t[i],S.y[i]); -} - -// 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. -solution 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) -{ - solution S; - S.t=new real[] {a}; - S.y=new real[]{y}; - - if(h == 0) { - if(b == a) return S; - 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; - 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); - - 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.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.a.highOrderWeights,predictions); - real y0=tableau.a.factors[tableau.a.steps.length]*y; - if(dynamic) { - real f1; - if(fsal) { - f1=f(t+h,y0+highOrder); - predictions.push(f1); - } - real lowOrder=h*dot(tableau.a.lowOrderWeights,predictions); - real error; - error=error(error,y,y0+lowOrder,y0+highOrder,highOrder-lowOrder); - h=adjust(h,error,tolmin,tolmax,tableau); - if(h >= dt) { - t += dt; - y=y0+highOrder; - S.t.push(t); - S.y.push(y); - f0=f1; - } - h=min(max(h,dtmin),dtmax); - } else { - t += h; - y=y0+highOrder; - S.y.push(y); - } - } - return S; -} - -struct Solution -{ - real[] t; - real[][] y; -} - -void write(Solution S) -{ - for(int i=0; i < S.t.length; ++i) { - write(S.t[i],tab); - for(real y : S.y[i]) - write(y,tab); - write(); - } -} - -// 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. -Solution 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) -{ - Solution S; - S.t=new real[] {a}; - S.y=new real[][] {copy(y)}; - - if(h == 0) { - if(b == a) return S; - if(n == 0) abort("Either n or h must be specified"); - else h=(b-a)/n; - } - real t=a; - real[] f0; - 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); - - 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.a.steps.length; ++i) - predictions.push(f(t+h*tableau.a.steps[i], - y+h*tableau.a.weights[i]*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.a.lowOrderWeights*predictions; - real error; - for(int i=0; i < y.length; ++i) - error=error(error,y[i],y[i]+lowOrder[i],y[i]+highOrder[i], - highOrder[i]-lowOrder[i]); - h=adjust(h,error,tolmin,tolmax,tableau); - if(h >= dt) { - t += dt; - y += highOrder; - S.t.push(t); - S.y.push(y); - f0=f1; - } - h=min(max(h,dtmin),dtmax); - } else { - t += h; - y += highOrder; - S.t.push(t); - S.y.push(y); - } - } - return S; -} - -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, - bool dynamic=false, real tolmin=0, real tolmax=0, real dtmin=0, - real dtmax=realMax, RKTableau tableau, bool verbose=false, - real[] initial(real[]), real[] discrepancy(real[]), - real[] guess, int iterations=100) -{ - real[] g(real[] t) { - real[][] y=integrate(initial(t),f,a,b,h,n,dynamic,tolmin,tolmax,dtmin,dtmax, - tableau,verbose).y;return discrepancy(y[y.length-1]); - } - real[][] jacobian(real[] t) {return finiteDifferenceJacobian(g,t);} - return initial(newton(iterations,g,jacobian,guess)); -} |