Exercise ODEs

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

   1. Implement an embedded Runge-Kutta stepper rkstepX (where X describes the order of the imbedded method used, for example X=12 ("one two", not "twelve") for the midpoint-Euler method) of your choice, which advances the solution to the equation

      dy/dt = f(t,y) ,

      (where y and f are vectors) by a given step h, and estimates the error.

      The interface could be like this,

      void rkstepX(
      	double t, double h, vector* y,
      	void f(double t, vector* y, vector* dydt),
      	vector* yh, vector* err
      )
      

      where vector*y contains the current values of the vector-function {yi(t)}, and the user supplied right-hand-side function void f(double t, vector*y, vector*dydt) fills-in the values dydt[k]=fk(t,y).

      The function rkstepX has to estimate the values yk(t+h) and store them in the array vector*yh and also estimate the errors δyk and store them in the array vector*err.

   2. Implement an adaptive-step-size driver routine wchich advances the solution from the given value of t to t=b (by calling your rkstepX routine with appropriate step-sizes) keeping the specified relative, eps, and absolute, acc, precision.

      The interface could be like this,

      void driver(
      	double* t, double b, double* h, vector*y, double acc, double eps,
      	void stepper(
      		double t, double h, vector*y,
      		void f(double t,vector*y,vector*dydt),
      		vector*yh, vector*err
      		),
      	void f(double t,double*y,double*dydt),
      )
      

      In the beginning
      • double*t contains the startpoint a;
      • double*h contains an estimate of the initial stepsize;
      • double*y contains y(a).
      In the end
      • double*t equals the endpoint b;
      • double*h contains the last accepted stepsize;
      • double*y contains y(b).

   3. Demonstrate that your routines work as intended by solving some interesting systems of ordinary differential equations.

2. (3 points) Storing the path

   Modify your driver such that it stores the points {ti,y1(ti),...,yn(ti)} it calculated along the way in, for example, a matrix [t,y1,y2,...,yn] provided by the user (or in any other convenient object).

3. (1 points) Non-embedded Runge-Kutta stepper with Runge's error estimate

   Use the same driver as in part-A/B, but implement a new stepper that uses the higher-order method from part-A/B and estimates the error using the Runge's principle (making a full step and two half-steps).

   Compare it with the your embedded method from exercise-A. Comparing here means counting the number of right-hand-side evaluations each of the methods used to achieve the same tolerance. The method with fewest evaluations wins. (Of course if the routine failed to reach the given tolerance, it is a failure -- ikke bestået.)

4. (0 point) A definite integral as an ODE

   1. Compare your implementations with your favourite library routines.

   2. A definite integral

      
      I = ∫abf(x)dx
      

      can be reformulated as an ODE,

      
      y'=f(x), y(a)=0 : I = y(b),
      

      which can be solved with your adaptive ODE solver.

      Implement this idea in a subroutine which calculates definite integrals based on your favourite ODE routine (to be later compared with your dedicated adaptive integrators).

      Calculate some interesting integrals.