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;}

// phi2(x)=(exp(x)-1-x)/(x^2);
// Use the identity phi2(2x)=0.25*(x*phi2(x)+1)^2+0.5*phi2(x);
real phi2(real x)
{
 if(fabs(x) > 1) return (exp(x)-x-1)/(x^2);
 x *= 0.125;
 real x2=x*x;
 real x3=x2*x;
 real x5=x2*x3;
 real y=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];
 y=0.25*(x*y+1.0)^2+0.5*y;
 y=(x*y+0.5)^2+0.5*y;
 return (2.0*x*y+0.5)^2+0.5*y;
}

// phi3(x)=(exp(x)-1-x-x^2/2)/(x^3)
// Use the identity phi3(2x)=0.125*phi2(x)*(x*phi2(x)+2)+0.25*phi3(x)
// where phi2(x)=x*phi3(x)+0.5
real phi3(real x)
{
 if(fabs(x) > 1.6) return (exp(x)-0.5*x^2-x-1)/x^3;
 x *= 0.125;
 real x2=x*x;
 real x3=x2*x;
 real x5=x2*x3;
 real y=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];
 real y2=x*y+0.5;
 y=0.125*y2*(x*y2+2)+0.25*y;
 y2=2*x*y+0.5;
 y=0.25*y2*(x*y2+1)+0.25*y;
 y2=4*x*y+0.5;
 return 0.25*y2*(2*x*y2+1)+0.25*y;
}

// phi4(x)=(exp(x)-1-x-x^2/2-x^3/6)/(x^4)
// Use the identity phi4(2x)=0.0625*(x*phi3(x)+0.5)^2+0.125*(phi3(x)+phi4(x));
// where phi3(x)=x*phi4(x)+1/6
real phi4(real x)
{
 if(fabs(x) > 1.6) return (exp(x)-Coeff[2]*x^3-0.5*x^2-x-1)/x^4;
 x *= 0.125;
 real x2=x*x;
 real x3=x2*x;
 real x4=x2*x2;
 real x5=x2*x3;
 real y=Coeff[3]+x*Coeff[4]+x2*Coeff[5]+x3*Coeff[6]+
   x4*Coeff[7]+x5*Coeff[8]+x3*x3*Coeff[9]+x5*x2*Coeff[10]+
   x4*x4*Coeff[11];
 real y3=x*y+Coeff[2];
 y=0.0625*(x*y3+0.5)^2+0.125*(y3+y);
 y3=2*x*y+Coeff[2];
 y=(0.5*x*y3+0.125)^2+0.125*(y3+y);
 y3=4*x*y+Coeff[2];
 return (x*y3+0.125)^2+0.125*(y3+y);
}

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 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[] {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) ? f :
   new real(real t, real y) {return f(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.t.push(t);
     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));
}

You are viewing proxied material from ftp.icm.edu.pl. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.