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
Something like this (Runge-Kutta Euler/Midpoint 12-method, you might want to do better),
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 */
)
{
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 δy = (k1-k0)*h; /* error estimate */
return (yh,δy);
}
Implement an adaptive-step-size driver routine wchich
advances the solution from initial-point a to final-point b
(by calling your rkstepXY routine with appropriate step-sizes)
keeping the specified relative, eps, and absolute,
acc, precision. You driver should record the solution in two
generic lists,
"genlist<double> x"
and
"genlist<vector> y" and then return the two lists.
The interface could be like this,
public static (genlist<double>,genlist<vector>) driver(
Func<double,vector,vector> F,/* the f from dy/dx=f(x,y) */
(double,double) interval, /* (initial-point,final-point) */
vector yinit, /* y(initial-point) */
double h=0.125, /* initial step-size */
double acc=0.01, /* absolute accuracy goal */
double eps=0.01 /* relative accuracy goal */
){
var (a,b)=interval; double x=a; vector y=yinit.copy();
var xlist=new genlist<double>(); xlist.add(x);
var ylist=new genlist<vector>(); ylist.add(y);
do{
if(x>=b) return (xlist,ylist); /* job done */
if(x+h>b) h=b-x; /* last step should end at b */
var (yh,δy) = rkstep12(F,x,y,h);
double tol = (acc+eps*yh.norm()) * Sqrt(h/(b-a));
double err = δy.norm();
if(err<=tol){ // accept step
x+=h; y=yh;
xlist.add(x);
ylist.add(y);
}
if(err>0) h *= Min( Pow(tol/err,0.25)*0.95 , 2); // readjust stepsize
else h*=2;
}while(true);
}//driver
Debug your routines by solving some interesting systems of ordinary differential equations, for example u''=-u.
Consider the equation (in certain units) of equatorial motion of a planet around a star in General Relativity,
u''(φ) + u(φ) = 1 + εu(φ)2 ,
where u(φ) ≡ 1/r(φ) , r is the (reduced-circumference) radial coordinate, φ is the azimuthal angle, ε is the relativistic correction (on the order of the star's Schwarzschild radius divided by the radius of the planet's orbit), and primes denote the derivative with respect to φ.
Integrate (for several rotations) this equation with ε=0 and initial conditions u(0)=1, u'(0)=0 : this should give a Newtonian circular motion.
Integrate (for several rotations) this equation with ε=0 and initial conditions u(0)=1, u'(0)≈-0.5 : this should give a Newtonian elliptical motion. Hint: u'(0) shouldn't bee too large or you will lose your planet.
Integrate (for several rotations) this equation with ε≈0.01 and initial conditions u(0)=1, u'(0)≈-0.5 : this should illustrate the relativistic precession of a planetary orbit.
Hints:
y0' = y1 y1' = 1-y0+ε*y0*y0
plot "data" using (1/$2)*cos($1):(1/$2)*sin($1) with lines notitle
Implement an imbedded stepper of order 23. Test whether it intergates the equation y''=2x exactly: solve this equation numerically, and check that the driver keeps doubling the step-size (as the error from the stepper should be close to zero) still producing the exact result.
Implement a quadratic spline interpolation routine that interpolates tables, {xi,yi}, of vector-valued functions y=f(x). Something like this,
public static Func<double,vector> make_qspline(genlist<double> x,genlist<vector> y)
{
/* calculate bi and ci of the quadratic spline */
Func<double,vector> qspline = delegate(double z){
int i=binsearch(x,z);
return y[i]+(z-x[i])*(b[i] + c[i]*(z-x[i]));;
};
return qspline;
}
public static Func<double,vector> make_ode_ivp_qspline
(Func<double,vector,vector> F,(double,double)interval,vector y,
double acc=0.01,double eps=0.01,double hstart=0.01 )
{
var (xlist,ylist) = driver(F,interval,y,acc,eps,hstart);
return make_qspline(xlist,ylist);
}
The dynamics of the Newtonian gravitational three-body problem is generally chaotic. However, there apparently exists a remarkable stable planar periodic solution in the shape of figure-8 [Wikipedia: Special-case solutions; Alain Chenciner, Richard Montgomeryi, arXiv:math/0011268].
Using your numerical ODE integrator reproduce this solution.
Hints: