←to practical programming

Exercise ODE

  1. (6 points) Embedded 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/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 rkstepXY(
	void f(double t,vector*y,vector*dydt) /* the f from dy/dt=f(t,y) */
	double t,              /* the current value of the variable */
	vector* yt,            /* the current value y(t) of the sought function */
	double h,              /* the step to be taken */
	vector* yh,             /* output: y(t+h) */
	vector* err             /* output: error estimate */
)

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

    2. 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.

      The interface could be like this, 

      void driver(
	void f(double,vector*,vector*), /* right-hand-side of dy/dt=f(t,y) */
	double a,                     /* the start-point a */
	vector* ya,                     /* y(a) */
	double b,                     /* the end-point of the integration */
	vector* yb,                     /* y(b) to be calculated */
	double h,                     /* initial step-size */
	double acc,                   /* absolute accuracy goal */
	double eps                    /* relative accuracy goal */
){ /* calculate y(b) */

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

    4. Consider the SIR model of epidemic development:
      • The population is divided into three compartments: Susceptible (S) – those who can be infected, Infectious (I), and Removed (R) – either recovered with aquired immunity, or dead.
      • The dynamics of the population (disregarding natural birth and death) is then described by the following set of ODE,

        dS/dt = - IS/NTc ,

        dI/dt = IS/NTc - I/Tr ,

        dR/dt = I/Tr ,

        where N is the population size, Tc is the typical time between contacts and Tr is the typical recovery time (Tr/Tc is the typical number of new infected by one infectious individual).

      Choose some reasonable model parameters for a country like Denmark, solve the ODE numerically with your routines, and investigate the effect of increasing Tc on the dynamics of the epidemic.

  2. (3 points) Store the path

    Make your driver store the points {ti, y(ti)} it calculated along the way. You can

    1. store the path in a matrix provided by the user;
    2. write the path to a file on the disk;
    3. build a suitable object, like a matrix where you can add extra rows at the bottom (man realloc), and store the path there;
    4. use a container from C++ (like std::vector): you need to compile you program with c++ (simply call your file *.cc and make should figure out) and LDLIBS+=-lstdc++.

  3. (1 point) Newtonian gravitational three-body problem