Implement an embedded Runge-Kutta stepper rkstepXY
(where XY describes the order of the imbedded method used, for
example XY=12—"one-two"—for the midpoint-Euler
method) of your choice, which advances the solution to the equation
dy/dx = f(x,y)
(where y and f are vectors) by a given steph,
and estimates the error.
Something like this,
public static (vector,vector) rkstep12(
Func<double,vector,vector> f /* the f from dy/dx=f(x,y) */
double x, /* the current value of the variable */
vector y, /* the current value y(x) of the sought function */
double h /* the step to be taken */
){ // Runge-Kutta Euler/Midpoint method (probably not the most effective)
vector k0 = f(x,y); /* embedded lower order formula (Euler) */
vector k1 = f(x+h/2,y+k0*(h/2)); /* higher order formula (midpoint) */
vector yh = y+k1*h; /* y(x+h) estimate */
vector er = (k1-k0)*h; /* error estimate */
return (yh,er);
}
Implement an adaptive-step-size driver routine wchich
advances the solution from a to b (by calling
your rkstepXY routine with appropriate step-sizes)
keeping the specified relative, eps, and absolute,
acc, precision.
The interface could be like this,
public static vector driver(
Func<double,vector,vector> f /* the f from dy/dx=f(x,y) */
double a, /* the start-point a */
vector ya, /* y(a) */
double b, /* the end-point of the integration */
double h=0.01, /* initial step-size */
double acc=0.01, /* absolute accuracy goal */
double eps=0.01 /* relative accuracy goal */
){
if(a>b) throw new Exception("driver: a>b");
double x=a; vector y=ya;
do {
if(x>=b) return y; /* job done */
if(x+h>b) h=b-x; /* last step should end at b */
var (yh,erv) = rkstep12(F,x,y,h);
double tol = Max(acc,yh.norm()*eps) * Sqrt(h/(b-a));
double err = erv.norm();
if(err<=tol){ x+=h; y=yh; } // accept step
h *= Min( Pow(tol/err,0.25)*0.95 , 2); // reajust stepsize
}while(true);
}//driver
Debug your routines by solving some interesting systems of ordinary differential equations, for example u''=-u.
• Check the acceptence condition for every component of y and investigate the tolerance/error ratio also for every component (pick the most difficult: the smallest, I suppose). Something like
tol[i]=Max(acc,Abs(yh[i])*eps)*Sqrt(h/(b-a));
bool ok=true;
for(int i=0;i<tol.size;i++) ok = (ok && err[i]<tol[i])
if(ok){ x+=h; y=yh; }
double factor = tol[0]/Abs(err[0]);
for(int i=1;i<tol.size;i++) factor=Min(factor,tol[i]/Abs(err[i]));
h *= Min( Pow(factor,0.25)*0.95 ,2);
• Make your driver store (optionally) the points {xi, y(xi)} it calculated along the way in two generic lists, say xlist and ylist. Use your own implementation of generic lists ;)