Exercise ODE

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,                                  /* the current value of the variable */
      	double h,                                  /* the step to be taken */
      	vector* yt,                                /* the current value y(t) of the sought function */
      	void f(double t, vector* y, vector* dydt), /* the right-hand-side, dydt = f(t,y) */
      	vector* yth,                               /* output: y(t+h) */
      	vector* err                                /* output: error estimate dy */
      )
      

      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,                             /* the current value of the variable */
      	double b,                              /* the end-point of the integration */
      	double* h,                             /* the current step-size */
      	vector*yt,                             /* the current y(t) */
      	double acc,                            /* absolute accuracy goal */
      	double eps,                            /* relative accuracy goal */
      	void stepper(                          /* the stepper function to be used */
      		double t, double h, vector*yt,
      		void f(double t,vector*y,vector*dydt),
      		vector*yth, vector*err
      		),
      	void f(double t,double*y,double*dydt), /* right-hand-side */
      )
      

      In the beginning
      • double*t contains the startpoint a;
      • double*h contains an estimate of the initial stepsize;
      • double*yt contains y(a).
      In the end
      • double*t equals the endpoint b;
      • double*h contains the last accepted stepsize;
      • double*yt 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 point) A definite integral as an ODE

   A definite integral

   I = ∫abf(x)dx
   

   can be reformulated as an ODE with initial condition,

   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 ODE routine (to be later compared with your dedicated adaptive integrators).

   Calculate some interesting integrals.