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 steph,
and estimates the error.
The interface could be like this,
static void rkstepXY( Func<double,vector,vector> f, /* the right-hand-side of 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.
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.
The interface could be like this,
static vector driver(
Func<double,vector,vector> f, /* right-hand-side of dy/dt=f(t,y) */
double a, /* the start-point a */
vector y, /* y(a) */
double b, /* the end-point of the integration */
double h /* initial step-size */
double acc, /* absolute accuracy goal */
double eps, /* relative accuracy goal */
){ /* return y(b) */
Demonstrate that your routines work as intended by solving some interesting systems of ordinary differential equations, for example u''=-u.
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).