johnpye/dopri5/dopri5.h

/*
    This is a modified form of the DOPRI5 code, originally written in Fortran 77
    by E Hairer and G Wanner, then translated to C by J Colinge.
    Changes have been made to suite ASCEND requirements. Minor only.
    -- John Pye, May 2007.
*//*
    DOPRI5
    Copyright (c) 2004, Ernst Hairer
 
    Redistribution and use in source and binary forms, with or without 
    modification, are permitted provided that the following conditions are 
    met:
 
    - Redistributions of source code must retain the above copyright 
    notice, this list of conditions and the following disclaimer.
 
    - Redistributions in binary form must reproduce the above copyright 
    notice, this list of conditions and the following disclaimer in the 
    documentation and/or other materials provided with the distribution.
 
    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS 
    IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED 
    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 
    PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR 
    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 
    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 
    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 
    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*//** @file
    This code computes the numerical solution of a system of first order ordinary
    differential equations y'=f(x,y). It uses an explicit Runge-Kutta method of
    order (4)5 due to Dormand & Prince with step size control and dense output.
 
    Authors : E. Hairer & G. Wanner
          Universite de Geneve, dept. de Mathematiques
          CH-1211 GENEVE 4, SWITZERLAND
          E-mail : HAIRER@DIVSUN.UNIGE.CH, WANNER@DIVSUN.UNIGE.CH
 
    The code is described in : E. Hairer, S.P. Norsett and G. Wanner, Solving
    ordinary differential equations I, nonstiff problems, 2nd edition,
    Springer Series in Computational Mathematics, Springer-Verlag (1993).
 
    Version of April 28, 1994.
 
    Remarks about the C version : this version allocates memory by itself, the
    iwork array (among the initial FORTRAN parameters) has been splitted into
    independant initial parameters, the statistical variables and last step size
    and x have been encapsulated in the module and are now accessible through
    dedicated functions, the variable names have been kept to maintain a kind
    of reading compatibility between the C and FORTRAN codes; adaptation made by
    J.Colinge (COLINGE@DIVSUN.UNIGE.CH).
 
    INPUT PARAMETERS
    ----------------
 
    n        Dimension of the system (n < UINT_MAX).
 
    fcn      A pointer the the function definig the differential equation, this
         function must have the following prototype
 
           void fcn (unsigned n, double x, double *y, double *f)
 
         where the array f will be filled with the function result.
 
    x        Initial x value.
 
    *y       Initial y values (double y[n]).
 
    xend     Final x value (xend-x may be positive or negative).
 
    *rtoler  Relative and absolute error tolerances. They can be both scalars or
    *atoler  vectors of length n (in the scalar case pass the addresses of
         variables where you have placed the tolerance values).
 
    itoler   Switch for atoler and rtoler :
           itoler=0 : both atoler and rtoler are scalars, the code keeps
                  roughly the local error of y[i] below
                  rtoler*abs(y[i])+atoler.
           itoler=1 : both rtoler and atoler are vectors, the code keeps
                  the local error of y[i] below
                  rtoler[i]*abs(y[i])+atoler[i].
 
    solout   A pointer to the output function called during integration.
         If iout >= 1, it is called after every successful step. If iout = 0,
         pass a pointer equal to NULL. solout must must have the following
         prototype
 
           solout (long nr, double xold, double x, double* y, unsigned n, int* irtrn)
 
         where y is the solution the at nr-th grid point x, xold is the
         previous grid point and irtrn serves to interrupt the integration
         (if set to a negative value).
 
         Continuous output : during the calls to solout, a continuous solution
         for the interval (xold,x) is available through the function
 
           contd5(i,s)
 
         which provides an approximation to the i-th component of the solution
         at the point s (s must lie in the interval (xold,x)).
 
    iout     Switch for calling solout :
           iout=0 : no call,
           iout=1 : solout only used for output,
           iout=2 : dense output is performed in solout (in this case nrdens
                must be greater than 0).
 
    fileout  A pointer to the stream used for messages, if you do not want any
         message, just pass NULL.
 
    icont    An array containing the indexes of components for which dense
         output is required. If no dense output is required, pass NULL.
 
    licont   The number of cells in icont.
 
 
    Sophisticated setting of parameters
    -----------------------------------
 
         Several parameters have a default value (if set to 0) but, to better
         adapt the code to your problem, you can specify particular initial
         values.
 
    uround   The rounding unit, default 2.3E-16 (this default value can be
         replaced in the code by DBL_EPSILON providing float.h defines it
         in your system).
 
    safe     Safety factor in the step size prediction, default 0.9.
 
    fac1     Parameters for step size selection; the new step size is chosen
    fac2     subject to the restriction  fac1 <= hnew/hold <= fac2.
         Default values are fac1=0.2 and fac2=10.0.
 
    beta     The "beta" for stabilized step size control (see section IV.2 of our
         book). Larger values for beta ( <= 0.1 ) make the step size control
         more stable. dopri5 needs a larger beta than Higham & Hall. Negative
         initial value provoke beta=0; default beta=0.04.
 
    hmax     Maximal step size, default xend-x.
 
    h        Initial step size, default is a guess computed by the function hinit.
 
    nmax     Maximal number of allowed steps, default 100000.
 
    meth     Switch for the choice of the method coefficients; at the moment the
         only possibility and default value are 1.
 
    nstiff   Test for stiffness is activated when the current step number is a
         multiple of nstiff. A negative value means no test and the default
         is 1000.
 
    nrdens   Number of components for which dense outpout is required, default 0.
         For 0 < nrdens < n, the components have to be specified in icont[0],
         icont[1], ... icont[nrdens-1]. Note that if nrdens=0 or nrdens=n, no
         icont is needed, pass NULL.
 
 
    Memory requirements
    -------------------
 
         The function dopri5 allocates dynamically 8*n doubles for the method
         stages, 5*nrdens doubles for the interpolation if dense output is
         performed and n unsigned if 0 < nrdens < n.
 
 
 
    OUTPUT PARAMETERS
    -----------------
 
    y       numerical solution at x=xRead() (see below).
 
    dopri5 returns the following values
 
         1 : computation successful,
         2 : computation successful interrupted by solout,
        -1 : input is not consistent,
        -2 : larger nmax is needed,
        -3 : step size becomes too small,
        -4 : the problem is probably stff (interrupted).
 
 
    Several functions provide access to different values :
 
    xRead   x value for which the solution has been computed (x=xend after
        successful return).
 
    hRead   Predicted step size of the last accepted step (useful for a
        subsequent call to dopri5).
 
    nstepRead   Number of used steps.
    naccptRead  Number of accepted steps.
    nrejctRead  Number of rejected steps.
    nfcnRead    Number of function calls.
*/

#include <stdio.h>
#include <limits.h>

typedef void FcnEqDiff(unsigned n, double x, double *y, double *f);
typedef void SolTrait(long nr, double xold, double x, double* y, unsigned n, int* irtrn);

extern int dopri5(
        unsigned n,      /* dimension of the system <= UINT_MAX-1*/
        FcnEqDiff *fcn,   /* function computing the value of f(x,y) */
        double x,        /* initial x-value */
        double* y,       /* initial values for y */
        double xend,     /* final x-value (xend-x may be positive or negative) */
        double* rtoler,  /* relative error tolerance */
        double* atoler,  /* absolute error tolerance */
        int itoler,      /* switch for rtoler and atoler */
        SolTrait *solout, /* function providing the numerical solution during integration */
        int iout,        /* switch for calling solout */
        FILE* fileout,   /* messages stream */
        double uround,   /* rounding unit */
        double safe,     /* safety factor */
        double fac1,     /* parameters for step size selection */
        double fac2,
        double beta,     /* for stabilized step size control */
        double hmax,     /* maximal step size */
        double h,        /* initial step size */
        long nmax,       /* maximal number of allowed steps */
        int meth,        /* switch for the choice of the coefficients */
        long nstiff,     /* test for stiffness */
        unsigned nrdens, /* number of components for which dense outpout is required */
        unsigned* icont, /* indexes of components for which dense output is required, >= nrdens */
        unsigned licont  /* declared length of icon */
    );

extern double contd5(
        unsigned ii,     /* index of desired component */
        double x         /* approximation at x */
    );

extern long nfcnRead (void);   /* encapsulation of statistical data */
extern long nstepRead (void);
extern long naccptRead (void);
extern long nrejctRead (void);
extern double hRead (void);
extern double xRead (void);

1	/*
2	This is a modified form of the DOPRI5 code, originally written in Fortran 77
3	by E Hairer and G Wanner, then translated to C by J Colinge.
4	Changes have been made to suite ASCEND requirements. Minor only.
5	-- John Pye, May 2007.
6	//
7	DOPRI5
8	Copyright (c) 2004, Ernst Hairer
9
10	Redistribution and use in source and binary forms, with or without
11	modification, are permitted provided that the following conditions are
12	met:
13
14	- Redistributions of source code must retain the above copyright
15	notice, this list of conditions and the following disclaimer.
16
17	- Redistributions in binary form must reproduce the above copyright
18	notice, this list of conditions and the following disclaimer in the
19	documentation and/or other materials provided with the distribution.
20
21	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS
22	IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
23	TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
24	PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR
25	CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
26	EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
27	PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
28	PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
29	LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
30	NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
31	SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32	//* @file
33	This code computes the numerical solution of a system of first order ordinary
34	differential equations y'=f(x,y). It uses an explicit Runge-Kutta method of
35	order (4)5 due to Dormand & Prince with step size control and dense output.
36
37	Authors : E. Hairer & G. Wanner
38	Universite de Geneve, dept. de Mathematiques
39	CH-1211 GENEVE 4, SWITZERLAND
40	E-mail : HAIRER@DIVSUN.UNIGE.CH, WANNER@DIVSUN.UNIGE.CH
41
42	The code is described in : E. Hairer, S.P. Norsett and G. Wanner, Solving
43	ordinary differential equations I, nonstiff problems, 2nd edition,
44	Springer Series in Computational Mathematics, Springer-Verlag (1993).
45
46	Version of April 28, 1994.
47
48	Remarks about the C version : this version allocates memory by itself, the
49	iwork array (among the initial FORTRAN parameters) has been splitted into
50	independant initial parameters, the statistical variables and last step size
51	and x have been encapsulated in the module and are now accessible through
52	dedicated functions, the variable names have been kept to maintain a kind
53	of reading compatibility between the C and FORTRAN codes; adaptation made by
54	J.Colinge (COLINGE@DIVSUN.UNIGE.CH).
55
56	INPUT PARAMETERS
57	----------------
58
59	n Dimension of the system (n < UINT_MAX).
60
61	fcn A pointer the the function definig the differential equation, this
62	function must have the following prototype
63
64	void fcn (unsigned n, double x, double y, double f)
65
66	where the array f will be filled with the function result.
67
68	x Initial x value.
69
70	*y Initial y values (double y[n]).
71
72	xend Final x value (xend-x may be positive or negative).
73
74	*rtoler Relative and absolute error tolerances. They can be both scalars or
75	*atoler vectors of length n (in the scalar case pass the addresses of
76	variables where you have placed the tolerance values).
77
78	itoler Switch for atoler and rtoler :
79	itoler=0 : both atoler and rtoler are scalars, the code keeps
80	roughly the local error of y[i] below
81	rtoler*abs(y[i])+atoler.
82	itoler=1 : both rtoler and atoler are vectors, the code keeps
83	the local error of y[i] below
84	rtoler[i]*abs(y[i])+atoler[i].
85
86	solout A pointer to the output function called during integration.
87	If iout >= 1, it is called after every successful step. If iout = 0,
88	pass a pointer equal to NULL. solout must must have the following
89	prototype
90
91	solout (long nr, double xold, double x, double* y, unsigned n, int* irtrn)
92
93	where y is the solution the at nr-th grid point x, xold is the
94	previous grid point and irtrn serves to interrupt the integration
95	(if set to a negative value).
96
97	Continuous output : during the calls to solout, a continuous solution
98	for the interval (xold,x) is available through the function
99
100	contd5(i,s)
101
102	which provides an approximation to the i-th component of the solution
103	at the point s (s must lie in the interval (xold,x)).
104
105	iout Switch for calling solout :
106	iout=0 : no call,
107	iout=1 : solout only used for output,
108	iout=2 : dense output is performed in solout (in this case nrdens
109	must be greater than 0).
110
111	fileout A pointer to the stream used for messages, if you do not want any
112	message, just pass NULL.
113
114	icont An array containing the indexes of components for which dense
115	output is required. If no dense output is required, pass NULL.
116
117	licont The number of cells in icont.
118
119
120	Sophisticated setting of parameters
121	-----------------------------------
122
123	Several parameters have a default value (if set to 0) but, to better
124	adapt the code to your problem, you can specify particular initial
125	values.
126
127	uround The rounding unit, default 2.3E-16 (this default value can be
128	replaced in the code by DBL_EPSILON providing float.h defines it
129	in your system).
130
131	safe Safety factor in the step size prediction, default 0.9.
132
133	fac1 Parameters for step size selection; the new step size is chosen
134	fac2 subject to the restriction fac1 <= hnew/hold <= fac2.
135	Default values are fac1=0.2 and fac2=10.0.
136
137	beta The "beta" for stabilized step size control (see section IV.2 of our
138	book). Larger values for beta ( <= 0.1 ) make the step size control
139	more stable. dopri5 needs a larger beta than Higham & Hall. Negative
140	initial value provoke beta=0; default beta=0.04.
141
142	hmax Maximal step size, default xend-x.
143
144	h Initial step size, default is a guess computed by the function hinit.
145
146	nmax Maximal number of allowed steps, default 100000.
147
148	meth Switch for the choice of the method coefficients; at the moment the
149	only possibility and default value are 1.
150
151	nstiff Test for stiffness is activated when the current step number is a
152	multiple of nstiff. A negative value means no test and the default
153	is 1000.
154
155	nrdens Number of components for which dense outpout is required, default 0.
156	For 0 < nrdens < n, the components have to be specified in icont[0],
157	icont[1], ... icont[nrdens-1]. Note that if nrdens=0 or nrdens=n, no
158	icont is needed, pass NULL.
159
160
161	Memory requirements
162	-------------------
163
164	The function dopri5 allocates dynamically 8*n doubles for the method
165	stages, 5*nrdens doubles for the interpolation if dense output is
166	performed and n unsigned if 0 < nrdens < n.
167
168
169
170	OUTPUT PARAMETERS
171	-----------------
172
173	y numerical solution at x=xRead() (see below).
174
175	dopri5 returns the following values
176
177	1 : computation successful,
178	2 : computation successful interrupted by solout,
179	-1 : input is not consistent,
180	-2 : larger nmax is needed,
181	-3 : step size becomes too small,
182	-4 : the problem is probably stff (interrupted).
183
184
185	Several functions provide access to different values :
186
187	xRead x value for which the solution has been computed (x=xend after
188	successful return).
189
190	hRead Predicted step size of the last accepted step (useful for a
191	subsequent call to dopri5).
192
193	nstepRead Number of used steps.
194	naccptRead Number of accepted steps.
195	nrejctRead Number of rejected steps.
196	nfcnRead Number of function calls.
197	*/
198
199	#include <stdio.h>
200	#include <limits.h>
201
202	typedef void FcnEqDiff(unsigned n, double x, double y, double f);
203	typedef void SolTrait(long nr, double xold, double x, double* y, unsigned n, int* irtrn);
204
205	extern int dopri5(
206	unsigned n, /* dimension of the system <= UINT_MAX-1*/
207	FcnEqDiff fcn, / function computing the value of f(x,y) */
208	double x, /* initial x-value */
209	double* y, /* initial values for y */
210	double xend, /* final x-value (xend-x may be positive or negative) */
211	double* rtoler, /* relative error tolerance */
212	double* atoler, /* absolute error tolerance */
213	int itoler, /* switch for rtoler and atoler */
214	SolTrait solout, / function providing the numerical solution during integration */
215	int iout, /* switch for calling solout */
216	FILE* fileout, /* messages stream */
217	double uround, /* rounding unit */
218	double safe, /* safety factor */
219	double fac1, /* parameters for step size selection */
220	double fac2,
221	double beta, /* for stabilized step size control */
222	double hmax, /* maximal step size */
223	double h, /* initial step size */
224	long nmax, /* maximal number of allowed steps */
225	int meth, /* switch for the choice of the coefficients */
226	long nstiff, /* test for stiffness */
227	unsigned nrdens, /* number of components for which dense outpout is required */
228	unsigned* icont, /* indexes of components for which dense output is required, >= nrdens */
229	unsigned licont /* declared length of icon */
230	);
231
232	extern double contd5(
233	unsigned ii, /* index of desired component */
234	double x /* approximation at x */
235	);
236
237	extern long nfcnRead (void); /* encapsulation of statistical data */
238	extern long nstepRead (void);
239	extern long naccptRead (void);
240	extern long nrejctRead (void);
241	extern double hRead (void);
242	extern double xRead (void);
243