Ordinary Differential Equations

The CCATSL library provides several routines for solving ODEs of first order,  $\dot x=a(x,t),$ and also for solving the more general system of ODEs of the form:

$$\renewedcommand{arraycolsep}{0.12em} \begin{array}{rcl} \dot... ...ts & \\ \dot x_n & = & a_n(t,x_1,\ldots{,}\,x_n). \end{array}$$

To handle a second order ODE (or a more complex system), you must first re-write the problem in the form above. For example, suppose you want to solve

$$\ddot x = t\dot x - 3t - x^2 + 1,$$

subject to  $x(t_0)=\alpha,$  $\dot x(t_0)=\beta.$ By introducing the variable $z(t),$ defined by $z(t)=\dot x(t),$ we can find the solution to the ODE above by finding the $x$-part of the solution to

$$\renewedcommand{arraycolsep}{0.12em} \begin{array}{rcl} \dot x & = & z \\ \dot z & = & tz - 3t - x^2 + 1 \end{array}$$

subject to $x(t_0)=\alpha,$ $z(t_0)=\beta.$

There are three CCATSL routines for ODE solving. The simplest, Rk4CL, is fast and easy to use, but employs a fixed stepsize and gives no indication of the error in the solution. The second, RkfCL, is more sophisticated; it adjusts the stepsize when necessary and allows control over the error introduced at each step. The final routine, Rkf45CL, also adjusts the stepsize automatically, but allows control of the error over the whole integration interval.

Fourth order Runge-Kutta, Rk4CL

This is the simplest method provided by CCATSL to solve ODEs and assumes that problem has been written in the form of a system of first-order ODEs. The CCATSL routine implementing Fourth order Runge-Kutta is called Rk4CL, and is declared as follows:

	void Rk4CL (	int n,
		double *ti,
		double dt,
		ODEFunctionCT F,
		double *x,
		double *xdot);

The arguments have the following meanings:

	n		The size of the system.
	ti		Pointer to a variable holding the initial time. When Rk4CL returns, it will have been updated to the time at the end of the integration interval.
	dt		The stepsize.
	F		A user-defined function which implements the system of ODEs above by computing  $\dot x_1,\ldots{,}\,\dot x_n$ from the values of $t$ and  $x_1,\ldots{,}\,x_n$). If something goes wrong, the function should set ErrorFlagCD=true before returning, which will cause Rk4CL to abort. (See below for an illustration.)
	x		Array used to hold the values of  $x_1,\ldots{,}\,x_n.$ Before calling Rk4CL for the first time, these should be set to the initial conditions. When Rk4CL returns, they will have been updated to hold the values of  $x_1,\ldots{,}\,x_n$ at the end of the integration interval.
	xdot		Array which will be used by the user-defined function F to pass back the values of  $\dot x_1,\ldots{,}\,\dot x_n$ to the library.

(If you're confused by pointers, arrays, or arguments declared as double *, don't worry. You can find the details here, but they are easier to use than to explain. If a routine asks for a `pointer to a double', you just setup a variable of type double as normal, and then put a & symbol in front of the variable name. The following example should help to clarify this.)

Here's a simple example using Rk4CL to solve an ODE and XYCurveCL to plot the solution:

/* /examples/chapter2/rk4demo.c */ #include <catam.h> void a(double t, double *x, double *xdot) { xdot[0]=x[1]; /* x-dot = z */ xdot[1]=t*x[1]-3*t-x[0]*x[0]+1; /* z-dot = tz -3t -x^2 +1 */ } int MainCL(void) { double dt=0.001; /* integration stepsize */ double x[2]; double xp[2]; double xdata[1000]; double tdata[1000]; double t; int i; t=0.0; /* initial time is t=0 */ x[0]=0.0; /* initial condition x=0 */ x[1]=0.0; /* initial condition z=0 */ for (i=0; i<1000; i++) { Rk4CL(2,&t,dt,a,x,xp); /* perform a single Rk4 step */ xdata[i]=x[0]; /* store x */ tdata[i]=t; /* store t */ } /* plot the solution */ XYCurveCL(tdata,xdata,1000,1,JOIN,RedCC,AUTOAXES); return 0; }

Another way of writing this program, which avoids storing the values of x and t, is to use XYDrawCL.

Runge-Kutta-Fehlberg, RkfCL

This is a more powerful technique for solving ODEs than fourth-order Runge-Kutta. It automatically changes the stepsize when necessary and allows control of the local error (the error introduced at each timestep).

	boolean RKfCL (	int n,
		double aberr,
		double relerr,
		double *t,
		double *dt,
		double dtmin,
		ODEFunctionCT F,
		double *x,
		double *xdot,
		int *nleft);

	n		The size of the system.
	aberr		The acceptable absolute error over each integration step.
	relerr		The acceptable relative error at each step.
	t		Pointer to a variable initially set by the user to be the time at the start of the integration interval. When RkfCL returns, it will have been incremented to reflect the time at the end of the integration interval.
	dt		Pointer to a variable holding the desired step length. This value may be changed if RkfCL finds that a smaller stepsize is necessary to achieve the required tolerances. The value held by dt when RKfCL returns may differ from the length of the integration interval.
	dtmin		The minimum stepsize. If RkfCL cannot meet the required tolerances without using a stepsize below than this, RkfCL aborts.
	F		A user-defined function which implements the system of ODEs above by computing  $\dot x_1,\ldots{,}\,\dot x_n$ from the values of $t$ and  $x_1,\ldots{,}\,x_n$). If something goes wrong, the function should set ErrorFlagCD=true before returning, which will cause RkfCL to abort. (See below for an illustration.)
	x		Array used to hold the values of  $x_1,\ldots{,}\,x_n.$ Before calling RkfCL for the first time, these should be set to the initial conditions. When RkfCL returns, they will have been updated to hold the values of  $x_1,\ldots{,}\,x_n$ at the end of the integration interval.
	xdot		Array used by the user-defined function F to pass back the values of  $\dot x_1,\ldots{,}\,\dot x_n$ to the library.
	nleft		Pointer to a variable used to indicate how many integration steps remain until some desired time is reached, based on the current dt. Usually you want to integrate up to some fixed time $T$; you would choose $N$ to be the desired number of integration steps, passing this via nleft and set the suggested step size to be $T/N$. RkfCL will then adjust $N$ appropriately whenever an integration step is completed, or when it changes dt.

The return value is true if the integration tolerances can be met and false if they cannot.

(If you're confused by pointers, arrays, or arguments declared as double *, don't worry. You can find the details here, but they are easier to use than to explain. If a routine asks for a `pointer to a double', you just setup a variable of type double as normal, and then put a & symbol in front of the variable name. The following example should help to clarify this.)

Here is an example of how to use RkfCL to solve the ODE of the Van der Pol oscillator:

$$\m \ddot x+ \mu \dot x (x^2-1)+x=0.$$

/* /examples/chapter2/rkfdemo.c */ #include <catam.h> double mu=5.0; void a(double t, double *x, double *xdot) { xdot[0]=x[1]; xdot[1]=-x[0]-mu*x[1]*(x[0]*x[0]-1); } int MainCL(void) { double aberr=1e-5; double relerr=1e-5; double dtmin=0.001; double tfinal=25.0; double x[2]; /* will store x and z */ double xp[2]; /* will store x-dot and z-dot */ double t; double xdata[1000]; double tdata[1000]; int nres; /* number of results stored */ double dt; /* integration stepsize */ int nleft; /* number of steps left */ t=0.0; /* setup initial conditions */ x[0]=-1.42; x[1]=0.26; dt=0.2; nleft=(int)(tfinal/dt); /* workout a reasonable stepsize */ dt=tfinal/nleft; nres=0; while (nleft>0) { if (!RkfCL(2,aberr,relerr,&t,&dt,dtmin,a,x,xp,&nleft)) { printf("tolerances can't be met, sorry"); HaltCL(); } if (nres<1000) { /* careful not to overflow arrays */ xdata[nres]=x[0]; tdata[nres]=t; nres=nres+1; } } /* plot the solution */ XYCurveCL(tdata,xdata,nres,1,JOIN,RedCC,AUTOAXES); return 0; }

Rkf45CL

This is the most advanced ODE solving routine in CCATSL. It automatically chooses, and if necessary varies, the integration timestep and allows for error control over the whole integration interval.

	void Rkf45CL (	int n,
		double aberr,
		double relerr,
		double tend,
		double *t,
		ODEFunction F,
		double *x,
		double *xdot,
		int *control);

	n		The size of the system.
	aberr		The absolute error tolerance. You can set aberr to zero and Rkf45CL will report back if this leads to a problem
	relerr		The relative error tolerance. relerr should be more than about 1e-13 but you can try smaller values.
	tend		The time at the end of the integration interval.
	t		Pointer to a variable holding the time at the start of the integration interval. When Rkf45CL returns, this will have been updated appropriately.
	F		A user-defined function which implements the system of ODEs above by computing  $\dot x_1,\ldots{,}\,\dot x_n$ from the values of $t$ and  $x_1,\ldots{,}\,x_n$). If something goes wrong, the function should set ErrorFlagCD=true before returning, which will cause Rkf45CL to abort. (See below for an illustration.)
	x		Array used to hold the values of  $x_1,\ldots{,}\,x_n.$ Before calling Rkf45CL for the first time, these should be set to the initial conditions. When Rkf45CL returns, they will have been updated to hold the values of  $x_1,\ldots{,}\,x_n$ at the end of the integration interval.
	xdot		Array used by the user-defined function F to pass back the values of  $\dot x_1,\ldots{,}\,\dot x_n$ to the library.
	control		See below.

The parameter control is a pointer to a variable used to control the behaviour of Rkf45CL. The first call to Rkf45CL should set the variable to 1 or -1. In the first case, Rkf45CL will try to integrate up to time tend in one go, meeting the required tolerances. The second case causes Rkf45CL to return after each substep, with t and the array x set accordingly. When Rkf45CL returns, the variable acts as a status indicator. Values of 2 and -2 correspond to successful integrations in the two cases just described (in the step-by-step mode, you should leave control at -2 when you next call Rkf45CL). Other possible values are:

3

relerr was too small, but has been changed to a more appropriate value, you may continue if you like.

4

More than 5000 steps have been taken and tend has still not been reached, you may continue if you like.

5

Solution has vanished and aberr=0, you probably want to set aberr to a non-zero value and do a single step here (by setting the variable pointer to by control to -2).

6

Unable to achieve the desired tolerances, change relerr or aberr and try again.

7

Rkf45CL thinks it wants a stepsize greater than the whole integration interval, the problem is probably stiff (CCATSL has no suitable routines for such problems).

8

The input settings were invalid.

An example showing the use of Rkf45CL is c01orbit.c