Homework "ODE"

Objective

Implement an embedded Runge-Kutta stepper with error estimate and an adaptive step-size driver for solving Ordinary Differential Equation (ODE) Initial Value Problems (IVP).

Tasks

(6 points) Embedded Runge-Kutta ODE integrator

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 step h, and estimates the error.

Something like this (Runge-Kutta Euler/Midpoint method),

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 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. 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 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 ArgumentException("driver: a>b");
double x=a; vector y=ya.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,erv) = rkstep12(F,x,y,h);
        double tol = (acc+eps*yh.norm()) * Sqrt(h/(b-a));
        double err = erv.norm();
        if(err<=tol){ // accept step
		x+=h; y=yh;
		xlist.add(x);
		ylist.add(y);
		}
	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.
Reproduce (with your routines) the example from the scipy.integrate.odeint manual (oscillator with friction).

(3 points) Driver improvements
• Quite often one is only interested in the solution-vector at the end-point of the integration interval. In this case recording the path is a waste of time. To address this case, modify the interface to your driver such that the driver records the path only if the user provides initialized generic lists for both "x" and "y". Otherwise, that is, if the lists are "null", the driver only returns the "y" at the end-point.
Something like
```
public static vector driver(
	Func<double,vector,vector> F,
	double a, vector ya, double b,
	double acc=1e-2, double eps=1e-2, double h=0.01,
	genlist<double> xlist=null, genlist<vector> ylist=null
	){...}
```
• Change the step acceptence condition: investigate the tolerance/error ratio separately for every component of the vector y and pick the smallest. Something like
```
for(int i=0;i<y.size;i++)
	tol[i]=(acc+eps*Abs(yh[i]))*Sqrt(h/(b-a));
bool ok=true;
for(int i=0;i<y.size;i++)
	if(! err[i]<tol[i]) ok=false;
if(ok){ x+=h; y=yh; }
double factor = tol[0]/Abs(erv[0]);
for(int i=1;i<y.size;i++) factor=Min(factor,tol[i]/Abs(erv[i]));
h *= Min( Pow(factor,0.25)*0.95 ,2);
```
• At your choice, either
- Reproduce the example from the scipy.integrate.solve_ivp manual (Lotka-Volterra system).
or solve the following problem:
- Consider the equation of equatorial motion (in certain units) of a planet around a star in General Relativity,
  u''(φ) + u(φ) = 1 + εu(φ)² ,
  where u(φ) ≡ ¹/_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 φ.
  1. Integrate this equation with ε=0 and initial conditions u(0)=1, u'(0)=0 : this should give a Newtonian circular motion.
  2. Integrate 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.
  3. Integrate 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:
  - You can rewrite this second order equation as a system of ordinary differential equations (ODE) by introducing the functions y₀=u and y₁=u', which gives the following ODE,
```
y₀' = y₁
y₁' = 1-y₀+ε*y₀*y₀ 
```
  - If you have calculated the data-file as two columns with φ and u, you can plot the orbit with the following gnuplot/pyxplot command,
```
plot "data" using (1/$2)*cos($1):(1/$2)*sin($1) with lines notitle 
```
(1 point) Newtonian gravitational three-body problem
• 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: the Newtonian equations of motion for several gravitating bodies with masses m_i is given as
m_i^{^dv_i}/_{_dt} = ∑_j≠i ( ^{^Gm_im_j} /_{_{|r_j-r_i|²}} ) ( ^{^r_j-r_i} /_{_{|r_j-r_i|}} )