DDESolve Help

Defining Models

To define a model you need to edit various functions in a standard template file. You may also need to use some standard functions. This section explains these. The following table summarises what you need to know about for each class of models that the package can deal with (SDDE stands for delay differential equations with switches):

Functions to Define	ODE	DDE	SDDE
initcons()	Y	Y	Y
switchfunctions()	-	-	Y
map()	-	-	Y
grad()	Y	Y	Y
storehistory()	-	Y	Y
initialstate()	Y	Y	Y
statescale()	Y	Y	Y
initout()	Y	Y	Y
Functions to use
pastvalue()	-	Y	Y
pastgradient()	-	Y	Y

State variables and constants

The programme aims to solve problems of the general form:

s₀' = g₀(s(t), s(t-t₁), s(t-t₂), ..., s(t-t_m), t)
s₁' = g₁(s(t), s(t-t₁), s(t-t₂), ..., s(t-t_m), t)
...
s_n' = g_n(s(t), s(t-t₁), s(t-t₂), ..., s(t-t_m), t)

subject to initial conditions on the s_i's, and with the possible addition of some discontinuities in the s_i's. The s_i's are the state variables of the problem and s(t) is used to mean the vector of all state variables at time t, so s(t-t_j) is the vector of state variables time t_j ago. To define a model the user must specify the functions g_i which return the gradients of s_i w.r.t. time given the state variables at time t and at a series of lagged times. The user must also specify the starting values for all state variables.

Within the program the state variables are passed to and from functions as an array s[]. The functions defining the gradients of the state variables may depend on various user defined constants (as may the initial values) - these constants are passed around in an array c[]. All arrays start at element 0, not element 1. Elements of s[] must not be modified anywhere, except in map() and initstate(). If you modify elements of c[] then the changes will be permanent until you next change them.

initcons(no_vars, no_cons)

The arguments of this function are all pointers to double, except t, the time, which is a double. sw is an array in which to return switchfunction values; s[] is the array of state variables; c[] is the array of constants. Elements of s[] must not be altered in this function.

Switches are discontinuous changes in state variables. Switchfunctions are functions that pass through zero, with a negative time derivative at the time that a switch is to occur. Switches have an index (starting at zero). Users define the switchfunctions in the switchfunctions() routine by setting the values of the switchfunctions in the array sw[]. For example, if the 0th switch is to occur every 33.2 time units starting at time 11.1 the following line could be used:

sw[0]= - sin(6.28318530718*(t-11.1)/33.2);

As another example, if switch 1 is to occur just once at time 55.67 then the following line is appropriate:

sw[1]=55.67-t;

Switchfunctions can depend on constants and state variables, although some care may be required to ensure a well defined model when switchfunctions are state variable dependent.

When the i^th switchfunction passes through zero, from positive to negative, then integration is suspended (in an orderly fashion) and the routine map() is called with i passed as argument swno. map() is where the user defines switches.

map(s,c,t,swno)

The arguments of this function are (pointer to double) s[] the array of state variables; (pointer to double) c[] the array of constants; (double)t the time; (integer) swno, the index of the switch whose switchfunction is passing through zero (with negative slope).

map() is where you must define the switches, i.e. the discontinuous changes in state variables s[] triggered by a switchfunction descending through zero ( swichfunctions() above). For example, in the following code state variables 0 and 1 are doubled by switch 0 every time it is triggered, while state variable 0 is multiplied by c[0]/t by switch 1 whenever it is triggered:

if (swno==0) {s[0]*=2;s[1]*=2;} else
if (swno==1) {s[0]*=c[0]/t;}

grad(g,s,c,t,pastvalue)

The arguments of this function are all pointers to double, except time t, which is a double; s[] is the state variable array; c[] is the array of model constants; g[] is the array in which the user must return the time derivatives of the state variables. Elements of s[] must not be altered in this function. It is very unlikely to be a good idea to change elements of c[] in this function.

There is an additional argument passed to the DDESolve version of grad, a pointer to the pastvalue() function. This is because all symbols must be defined in the dynamic library.

This is the routine in which the user specifies the gradients of the state variables at time t, given the the state variables at time t (supplied in s[]), the model constants, and (optionally) lagged values of the state variables (or lagged functions of state variables). For each state variable s[i] you must supply a value in g[i]. For example if the equation for the gradient of state variable 0 is:

s₀' = c₀ s₀ - c₁ s₀ s₁

then the corresponding line in the routine grad() would be:

g[0]=c[0]*s[0]-c[1]*s[0]*s[1];

Lagged variables are accessed using the routinespastvalue() and (occasionally) pastgradient(). These are described later, but for completeness an example of their use is given here. Suppose that you want to code the equation:

s₀' = { c₀s₀(0) t< c₁

c₀s₀(t-c₁) t &ge c₁

subject to the initial condition s₀(0)=c₂. Assuming that the 0^th history variable (see storehistory() and pastvalue(), below) contains the lagged values of s[0] then the appropriate code would be:

if (t < c[1]) g[0]=c[0]*c[2];
else g[0]=c[0]*pastvalue(0,t-c[1],0);

Notes

If the gradient changes discontinuously, then it is often good practice to specify a switch at the time of the discontinuity: this forces the integration to step to exactly the point of discontinuity, before continuing, which ensures that the continuity assumptions of the integrator are not violated (the specified switch should leave all state variables unchanged).
The specification of your model should never involve writing to global variables from within grad(): grad() is called multiple times per step of the integrator, and t may increase or decrease between calls. Similarly writing to static variables is very rarely appropriate as part of model specification.
g[] contains meaningless values on entry to the function.

storehistory(his,ghis,g,s,c,t)

The arguments of this function are all pointers to double except for time t which is a double. g[] and s[] contain the current values of the gradients of the state variables (w.r.t. time) and the state variables; their elements must not be altered by this function. c[] is the array of model constants. his[] and ghis[] are arrays in which you should store the value and time derivative, at time t, of any quantity that you want to use as a lagged variable within your model.

If your model equations depend on lagged quantities then these must be stored. storehistory() is the function that allows you to do this. The lagged variables (or ``history variables'') are stored only at discrete times. Interpolation is necessary to estimate the values of lagged variables between these storage times. In order to ensure correct adaptive stepping for integration this interpolation must be performed to a higher order of accuracy than the integration. To achieve this requires that both the values and time derivatives of the state variables be stored¹. The lagged values and gradients for each history variable are stored in a ringbuffer for access by pastvalue() or occasionally pastgradient() (see below).

As an example, if your model is the simple DDE:

s₀' = { c₀s₀(0) t< c₁

c₀s₀(t-c₁) t &ge c₁

then history variable 0 would be defined as s₀. The appropriate piece of code within storehistory() would be:

his[0]=s[0];
ghis[0]=g[0];

and the within grad() the 0^th state variable would be used to specify the gradient of s₀ as follows:

if (t<c[1]) g[0]=c[0]*c[2];
else g[0]=c[0]*pastvalue(0,t-c[1],0);

History variables do not have to be state variables themselves (although it's usually easiest to set models up that way), for example if history variable 1 was to be exp(s₀) + 2s₁ then the appropriate code in storehistory() would be:

his[1]=exp(s[0])+2*s[1];
ghis[1]=g[0]*exp(s[0])+2*g[1];

If you must write to global variables, this routine is a reasonable place to do it, but be aware that you may compromise the accuracy of numerical model solution (since the adaptive stepping algorithm has no way of `seeing' this variable).

initst(s,c,t)

This function is called initstate in solv95: however, this clashes with a Unix library routine. Hence, it has been renamed initst, which is in line with the usage in ddefit.

The arguments of this routine are: (pointer to double) s, the state variable array; (pointer to double) c the constant array; (double) t the time at start of model integration.

This function is where the initial values of all state variables must be set by the user. For example to set state variable 0 to c₃ and state variable 1 to 3.45, you would use the following code:

s[0]=c[3];
s[1]=3.45;

It is quite safe to modify global variables and elements of c[] within this function.

statescale(double *scale)

DDESolve controls integration error by specifying that it should be less than some specified proportion of each state variable. This can cause difficulties when some state variables are in the vicinity of zero. For example, if a state variable is at zero, but about to leave that state, then the integrator may attempt to estimate that variable with zero error - which requires a zero timestep! The solution adopted is to get the user to specify a small number to be added to the magnitude of each state variable when estimating proportianal integration error. These numbers (one for each state variable) must be supplied by the user in array scale[].

Any elements of scale[] that are unspecified are taken as zero.

initout(usercontrol *out)

This routine is where you provide various initialization information for setting up the model and the display of the model's solutions. The routine is used to fill out the structure out: there is a large amount to specify, but it is straightforward to do so. The following sections cover the elements of out: any element marked by a * is optional.

Model control constants and defaults

out->nhv the (integer) number of history variables in the model.
out->hbsize the (integer) number of lagged values to store for each history variable, the larger this number, the larger the model time lags can be.
out->nlag this is the (integer) number of distinct lags at which each particular history variable is to be accessed. For example if you want to access history variable 0 at lag t-12.0, and history variable 1 at t-3.0, t-20.0 and t-13.2, then you would set out->nlag=3. All out->nlag does is to control the number of distinct place markers used in the history buffer for each state variable - the place markers are used to save time when searching for the right place in the history buffer. This feature is purely to allow efficient access of the history buffer - so you can always set out->nlag=1, it's just that this may be a bit inefficient.
out->nsw the number of switches (which is also the number of switch functions).
out->dout the output timestep - determines the approximate time between outputting results to windows and any specified file.
out->tol the integration tolerance - that is the maximum error tolerated at each timestep (as a proportion of the state variable concerned).
out->dt the default initial timestep this isn't too crucial, as the timestep will be rejected if it's too big. But note that the maximum timestep is set to 100*out->dt and the minimum to 1*10^-9*out->dt (I know, I know - you can change it in dde() in file ddeq.c.)
out->t0 the default integration start time.
out->t1 the default integration end time.

constants c[]

out->c[i] the default value for the i^th constant c[i] i (starts at zero).
out->cname[i] name for the i^th constant. e.g. to set the name of constant 0 to ``fred'' use the line out->cname[0]="fred";.
out->cinfo[i] extra information about a constant that can be displayed at run time by clicking on the constant's name. e.g.
```
  out->cinfo[0]="This constant measures the rate of accumulation of \
frustration when writing documentation";
```

State variable labels and other information

out->label[i] label for the i^th state variable. e.g. out->label[1]="population of adults". It's only worth supplying these for state variables that you intend to output.
out->initialtitle the title to be displayed in the About box.
out->initialtext explanatory text to be displayed in the About box.

Windowed output

DDESolve allows the user to control how many windows output will be displayed in, as well as how many and which state variables will be displayed in each window.

out->xlabel the label to be attached to the `time' axis of each output plot.
out->no_windows the (integer) number of windows for displaying output.
out->lines[i] the number of state variables to be plotted in the i^th window.
out->index[i].win the index of the window in which the i^th state variable is to be displayed - only supply this for those state variables that are to be displayed.
out->index[i].cur the index of the curve within a particular window which will display state variable i. e.g. to display state variable 0 as curve 1 of window 2 use the lines:
```
out->index[0].win=2;
out->index[0].cur=1;
```
out->wname[i] the name of the i^th window (displayed in the window header bar). \item out->range[i].y0 the initial minimum y axis value for window i. \item out->range[i].y1 the initial maximum y axis value for window i.

File output

A text output file can be produced, it's columns are time (or whatever variable is being integrated with respect to) followed by those state variables specified for output.

out->fileno the integer number of state variables to output.
out->fout[i] the index of the state variable to be output in column i+1 of the output file (i starts at zero, column zero of the output file is time).

When you run the model, you will be presented wth a save file panel to select the output file name. If this annoys you every time you run a model, set out->fileno = 0. You can still save model results using the Model->Write Results... menu in DDESolve.

pastvalue(i,t,j)

This function is called by the user to access history variables at lagged times. i is the index of the history variable (starting at zero); t is the (double) time at which the history variable is to be evaluated; j is the index of the history buffer location marker being used at this call. The function returns a double. Examples of the use of this function are given in the sections on storehistory() and grad() above, so here a slightly more complicated example is given. Consider the contrived situation of a population that produces half its offspring in an environment promoting fast maturation and half in an environment promoting slow maturation, suppose also that the maturation time varies throughout the year. Suitable equations might be:

s₀' =  (1/2)c₀s₀(t-T₀) + (1/2)c₀s₀(t-T₁) - c₁s₀(y) 
T₀ = 20 + 10 sin(2pi(t-90)/365)
T₁ = 2T₀

Suitable code for this in grad() might be:

. 
. 
    T0 = 20.0+10.0*sin(6.28318530718*(t-90.0)/365.0);
    T1 = 2*T0;
    if (t>T1) {
        g[0]=pastvalue(0,t-T1,0)*0.5*c[0] + pastvalue(0,t-T2,1)*0.5*c[0] - c[1]*s[0];
    } else {
.
.

Note the important point that you can not request lagged variables from before the start of integration - instead you have to define your model in such a way that it does not require such values: this is always possible for a well defined model.

pastgradient(i,t,j)

This is exactly the same as pastvalue() except that it returns the time derivative of the lagged variable - its use is somewhat dubious, since the order of approximation of the interpolated gradients is lower than the order of approximation of the integrator which may lead to infelicities in the numerical solution of models using this function.

Differences between the Mac and Windows versions

The "windows.h" file should not be included.
The grad function takes an extra argument pastvalue, a pointer to a function.

Coding hints

One use of switches, which is not obvious, is to make sure that the integrator ``sees'' short events - for example, if a population is at equilibrium it is possible for the integrator to be striding along taking such big steps that it steps right over a short pulse of immigrants - to get around this, you either need to modify the maximum timestep (see above) or put a switch in place somewhere during the pulse (which will force the integrator to stop at that point, thereby ``noticing'' the pulse).
You may want to write information to the screen from within your model - you can use NSLog() to do this. As a rather contrived example:
```
	NSLog(@"Initial Title = %@", [NSString stringWithUTF8String:out->initialtitle]);
```
Another way is to use a Cocoa ``NSRunAlertPanel''. Here is an example of a snippet of code that does just that:
```
	NSRunAlertPanel(@"Info!", @"Total host pop = %g", @"Dismiss", @"", @"", totpop);
```
You must set up your project to be a Cocoa dynamic library to use NSLog & NSRunAlertPanel. For pure C models, use printf();

Back to index

This consistency of integrator and interpolator is solv95's chief advantage over other DDE solving packages

s₀' = {	c₀s₀(0)	t< c₁
s₀' = {	c₀s₀(t-c₁)	t &ge c₁