←to practical programming

Homework "ODE"

Objective

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

Tasks

  1. (6 points) Embedded rule Runge-Kutta ODE integrator

    1. 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 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);
}

    2. Implement an adaptive-step-size driver routine wchich advances the solution from start-point a to end-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,    /* (start-point,end-point) */
	vector ystart,               /* y(start-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=ystart.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);
		}
	h *= Min( Pow(tol/err,0.25)*0.95 , 2); // readjust stepsize
        }while(true);
}//driver

    3. Debug your routines by solving some interesting systems of ordinary differential equations, for example u''=-u.

    4. Reproduce (with your routines) the example from the scipy.integrate.odeint manual (oscillator with friction) and/or the example from the scipy.integrate.solve_ivp manual (Lotka-Volterra system).

  2. (3 points) Alternative interface

    1. Change the interface of your integrator: make it return the interpolant of the solution based on the path recorded by your driver. For that you might need to make an interpolation routine (based on you own interpolation routines) that interpolate vector-valued functions.

      Something like this (using linear interpolation, you might want to do better: quadratic, for example), 

      public static Func<double,vector> make_linear_interpolant(genlist<double> x,genlist<vector> y)
{
	Func<double,vector> interpolant = delegate(double z){
		int i=binsearch(x,z);
		double Δx=x[i+1]-x[i];
		vector Δy=y[i+1]-y[i];
		return y[i]+Δy/Δx*(z-x[i]);
	};
	return interpolant;
}

      public static Func<double,vector> make_ode_ivp_interpolant
(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_linear_interpolant(xlist,ylist);
}

    2. Consider the equation of equatorial motion (in certain units) 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 φ.

      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 y0=u and y1=u', which gives the following ODE, 
        y0' = y1
y1' = 1-y0+ε*y0*y0
      • 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
  3. (1 points) 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: