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) of your choice, which advances the solution of 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 a to b (by calling your rkstepX routine with appropriate step-sizes) keeping the specified relative, eps, and absolute, acc, precision. Your driver must store 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).

      The interface could be like this,

      void driver(
      	double a, double b, double h_start, vector*y_a, 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),
      	matrix* t_y_table)
      

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

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

   Use the same driver as in part-A, but implement a new stepper which uses the higher-order method from part-A 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.)

3. (1 point) A definite integral as an ODE

   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.

4. (0 points) Compare your implementations with your favourite library routines.