Skip to content

# pinard/RKE

### Subversion checkout URL

You can clone with
or
.

# Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also compare across forks.

# Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also compare across forks.
...
• 3 commits
• 5 files changed
• 0 commit comments
• 1 contributor
Commits on Jun 03, 2008
 pinard .gitignore: New file aa2318c pinard README: Reformulated from reST to Org cc7073b pinard rke.h, rke.c: Do not declare alloca nor malloca, include i… …nstead 7400ca3
Showing with 92 additions and 209 deletions.
1. +2 −0  .gitignore
2. +0 −2  Makefile
3. +89 −193 README
4. +0 −12 rke.c
5. +1 −2  rke.h
2  .gitignore
 @@ -0,0 +1,2 @@ +*.o +example
2  Makefile
 @@ -20,13 +20,11 @@ CC = gcc CFLAGS = -g LDFLAGS = - example: example.o rke.o $(CC)$(LDFLAGS) example.o rke.o -o example -lm example.o: rke.h example.c rke.o: rke.h rke.c - clean: rm -f *~ *.o example
282 README
 @@ -1,212 +1,108 @@ -========================================= -README for the Runge-Kutta-England solver -========================================= - - :Author: François Pinard - :Email: pinard@iro.umontreal.ca - :Date: 1988-02 - :Copyright: © 1988 Free Software Foundation, Inc. - :Copying: See file COPYING for copying conditions. - -.. contents:: -.. sectnum:: - -.. _example.c: /rke/example.c -.. _rke.c: /rke/rke.c -.. _rke.h: /rke/rke.h - -To install, put yourself in the "rke" directory and type make. -To use, link "rke.o" with your program. File "rke.texinfo" contains -the documentation in GNU documentation format. Send any bug reports, -comments or criticisms to me at mailto:pinard@iro.umontreal.ca. - -General presentation -==================== - -This section presents briefly the theory of the problem. It also explains -the main characteristics of the solver. - -The user provides the problem to be solved as a problem function that -computes the derivatives, given estimates of the independent and all -dependent variables. Even if step size adjustment and error control -are somewhat automatic, the user may tune parameters of this automatic -behaviour if necessary. - -Definition of the problem -------------------------- - -This module approximates a numerical solution to problems that may be -expressed as systems of ordinary differential equations of the first -order. A simple trick may be used to work on higher order systems, -see the third example in example.c_. Using this module is easy -and requires no special knowledge, besides the intuitive notion of a -derivative. - -The system to be solved looks like:: - - d1(t) = f1(t, v1(t), v2(t), ..., vn(t)) - d2(t) = f2(t, v1(t), v2(t), ..., vn(t)) - ... - dn(t) = fn(t, v1(t), v2(t), ..., vn(t)) - -with the special notations:: - - d1(t) meaning the value of [d v1(t) / dt] at time t - d2(t) meaning the value of [d v2(t) / dt] at time t - ... - dn(t) meaning the value of [d vn(t) / dt] at time t - -and given a consistent initial state:: - - t, v1(t), v2(t), ..., vn(t). - -By consistent, I mean that the values v1, v2, ..., vn are all taken -at the same time t. Note that even if the variable t is often -used to represent the concept of time, this is absolutely artificial and -this variable t may be used for anything useful, of course. - -The function defining the problem ---------------------------------- - -The user should provide a single function which, given t and an -array of values for v1, v2, ..., vn, computes f1, f2, ..., fn -and returns an array of values for the derivatives d1, d2, ..., -dn. This function defines the problem to solve. See examples in -example.c_. - -Note that the problem function, as implemented, should not modify -anything in its argument array of values v1, v2, ..., vn. It should -not depend for its computations on any previous result stored into -its result array d1, d2, ..., dn; this array should be considered -undefined at the beginning of the problem function. Be careful: in C, -arrays are counted from 0, not 1. - -Note also that the problem function should normally return a non-zero -value. A zero value would indicate that some derivative could not be -computed and would cause the solver routine to return an error. - -The solver routine, and error control -------------------------------------- - -The solver routine will use this problem function in the process of -estimating other consistent sets of values of the functions v1(t), -v2(t), ..., vn(t) for any chosen value of t. The routine -receives a consistent set of values and a goal time, it will update -in place the consistent set of values into another consistent set of -values, in such a way that the updated t value will be very near the -goal time. The updated t will not necessary equal the goal time, but -should not be further than half the minimum step size. - -The routine implements a variable step size that adjusts automatically -during the computations, so to limit accumulated errors. If the estimated -error for any variable goes over some prescribed limit, the step size is -shortened and the computation attempted again. If all the estimated errors -are comfortably low, the step size is lengthened so to increase the -efficiency of the computations. The user supplied function is called 9 -times per successful step and 7 times per rejected step. - -The number of steps required augments linearly with the "distance" t -has to move and linearly with the inverse of the step size. More steep -the functions, more steps needed. - -There are a minimum and a maximum allowable step sizes, which the user -may change at any time using the pointer into the reentrency variables. -The minimum step size should not be zero. If the minimum or the maximum -step size is too small, the CPU consumption may increase drastically, -especially if the module gets desperate, numerical problems may also -develop due to the limited representation of real numbers. If the -minimum step size is too big, the module might declare itself incapable -of limiting the integration errors of fluctuating functions. If the -maximum step size is too big, some sharp phenomenon may go unnoticed by -the module. - -The prescribed maximum error is the sum of a maximum relative error and -a maximum absolute error, and is represented by the same linear function -for each of the integration variables. The slope of the linear function -tells the maximum relative error per unit of time. The bias of the -linear function tells the maximum absolute error per unit of time. If -the two coefficients of the linear function are too low, the solver will -tend to use a lot of CPU. Those coefficients may be changed at any time -between two solve calls using the pointer into the reentrancy variables. - -Reentrency considerations -------------------------- - -The initialization routine allocates a new block of variables for a -problem and returns a pointer to the block. The termination routine -deallocates the block. The solver routine receives this pointer and -any consistent set of values t, v1(t), v2(t), ..., vn(t). It is -possible and even easy to work simultaneously with different systems -and/or problems, provided that each system keeps for itself the pointer -to its block of variables and its consistent set of values. - -References and history ----------------------- +#+TITLE: RKE solver +#+SETUPFILE: ~/fp/notes/setup/pinard.org -A reference on Runge-Kutta integration might be: + /The Runge-Kutta-England tool is a differential equation solver written in C, which I used a few times in continuous simulation contexts. This is a fourth order Runge-Kutta solver combined with an heuristic for implementing an adaptative step size./ - Moursund, David G.; Duris, Charles S.: Elementary Theory and - Application of Numerical Analysis, McGram-Hill, 1967. +* Installation -The main part of the algorithm has been taken from an appendix in: +This is an old project of mine, ported to C in 1988-02. To install, fetch a copy and call *make*, like this: + +#+BEGIN_SRC sh +git clone git://github.com/pinard/RKE.git +cd RKE +make +#+END_SRC + +The reentrency variables are declared in [[https://github.com/pinard/RKE/blob/master/rke.h][rke.h]], and the routines are found within [[https://github.com/pinard/RKE/blob/master/rke.c][rke.c]]. To use, link /rke.o/ with your program, or modify the [[https://github.com/pinard/RKE/blob/master/Makefile][Makefile]] at will. Send any bug reports, comments or criticisms to me at mailto:pinard@iro.umontreal.ca. + +* General presentation + +This section presents briefly the theory of the problem. It also explains the main characteristics of the solver. + +The user provides the problem to be solved as a problem function that computes the derivatives, given estimates of the independent and all dependent variables. Even if step size adjustment and error control are somewhat automatic, the user may tune parameters of this automatic behaviour if necessary. + +** Definition of the problem + +This module approximates a numerical solution to problems that may be expressed as systems of ordinary differential equations of the first order. A simple trick may be used to work on higher order systems, see the [[Box%20slowing%20by%20friction%20in%20air][third example]]. Using this module is easy and requires no special knowledge, besides the intuitive notion of a derivative. + +The system to be solved looks like: + +#+BEGIN_VERSE +$d_1(t) = f1(t, v1(t), v2(t), …, vn(t))$ +$d2(t) = f2(t, v1(t), v2(t), …, vn(t))$ +… +$dn(t) = fn(t, v1(t), v2(t), …, vn(t))$ +#+END_VERSE + +with the special notations: + +#+BEGIN_EXAMPLE +d_1(t) meaning the value of [d v1(t) / dt] at time t +d2(t) meaning the value of [d v2(t) / dt] at time t +… +dn(t) meaning the value of [d vn(t) / dt] at time t +#+END_EXAMPLE + +and given a consistent initial state $t, v1(t), v2(t), …, vn(t)$. - Pritsker, A. Alan B.: The GASP IV Simulation Language, Jonh - Wiley & Sons, 1974. +By consistent, I mean that the values $v1, v2, …, vn$ are all taken at the same time $t$. Note that even if the variable $t$ is often used to represent the concept of time, this is absolutely artificial and this variable $t$ may be used for anything useful, of course. -I translated it from FORTRAN to Cyber Pascal in 1982, making the module -usable for a few applications in pure continuous simulation. In 1988, I -translated the Pascal module into C, making slight improvements and trying -to conform to GNU standards. +** The function defining the problem -Routines of the package -======================= +The user should provide a single function which, given $t$ and an array of values for $v1, v2, …, vn$, computes $f1, f2, …, fn$ and returns an array of values for the derivatives $d_1, d2, …, dn$. This function defines the problem to solve. See [[Some%20example%20programs][example programs]] below. -See the reentrency variables in rke.h_. +Note that the problem function, as implemented, should not modify anything in its argument array of values $v1, v2, …, vn$. It should not depend for its computations on any previous result stored into its result array $d_1, d2, …, dn$; this array should be considered undefined at the beginning of the problem function. Be careful: in C, arrays are counted from 0, not 1. -See the individual routines headers in rke.c_. +Note also that the problem function should normally return a non-zero value. A zero value would indicate that some derivative could not be computed and would cause the solver routine to return an error. + +** The solver routine, and error control + +The solver routine will use this problem function in the process of estimating other consistent sets of values of the functions $v1(t), v2(t), …, vn(t)$ for any chosen value of $t$. The routine receives a consistent set of values and a goal time, it will update in place the consistent set of values into another consistent set of values, in such a way that the updated $t$ value will be very near the goal time. The updated t will not necessary equal the goal time, but should not be further than half the minimum step size. + +The routine implements a variable step size that adjusts automatically during the computations, so to limit accumulated errors. If the estimated error for any variable goes over some prescribed limit, the step size is shortened and the computation attempted again. If all the estimated errors are comfortably low, the step size is lengthened so to increase the efficiency of the computations. The user supplied function is called 9 times per successful step and 7 times per rejected step. + +The number of steps required augments linearly with the "distance" $t$ has to move and linearly with the inverse of the step size. More steep the functions, more steps needed. + +There are a minimum and a maximum allowable step sizes, which the user may change at any time using the pointer into the reentrency variables. The minimum step size should not be zero. If the minimum or the maximum step size is too small, the CPU consumption may increase drastically, especially if the module gets desperate, numerical problems may also develop due to the limited representation of real numbers. If the minimum step size is too big, the module might declare itself incapable of limiting the integration errors of fluctuating functions. If the maximum step size is too big, some sharp phenomenon may go unnoticed by the module. + +The prescribed maximum error is the sum of a maximum relative error and a maximum absolute error, and is represented by the same linear function for each of the integration variables. The slope of the linear function tells the maximum relative error per unit of time. The bias of the linear function tells the maximum absolute error per unit of time. If the two coefficients of the linear function are too low, the solver will tend to use a lot of CPU. Those coefficients may be changed at any time between two solve calls using the pointer into the reentrancy variables. + +** Reentrency considerations + +The initialization routine allocates a new block of variables for a problem and returns a pointer to the block. The termination routine deallocates the block. The solver routine receives this pointer and any consistent set of values $t, v1(t), v2(t), …, vn(t)$. It is possible and even easy to work simultaneously with different systems and/or problems, provided that each system keeps for itself the pointer to its block of variables and its consistent set of values. + +** References and history + +A reference on Runge-Kutta integration might be: + +#+BEGIN_QUOTE +Moursund, David G.; Duris, Charles S.: /Elementary Theory and Application of Numerical Analysis/, McGram-Hill, 1967. +#+END_QUOTE + +The main part of the algorithm has been taken from an appendix in: -Some example programs -===================== +#+BEGIN_QUOTE +Pritsker, A. Alan B.: /The GASP IV Simulation Language/, Jonh Wiley & Sons, 1974. +#+END_QUOTE -Here is a small collection of examples to the usage of rke, given here to -supplement the documentation given previously. +I translated RKE from FORTRAN to Cyber Pascal in 1982, making it usable for a few applications in pure continuous simulation. In 1988, I translated the [[http://pinard.progiciels-bpi.ca/bonjour/rke.html][Pascal module]] into C, making slight improvements and trying to conform to GNU standards. For a while, it has been part of my [[http://fp-etc.progiciels-bpi.ca/index.html][FP etc.]] project, but now has [[https://github.com/pinard/RKE/][its own home]] on GitHub. -The file example.c_ contains the exact code to solve each of these. +* Some example programs -Integration under a normal curve --------------------------------- +File [[https://github.com/pinard/RKE/blob/master/example.c][example.c]] contains a solution examples for three small problems illustrating the usage of this RKE solver. These examples are explained below. -Simple integration is a particular case of solving ordinary differential -equations. Let us compute the probability that a normally distributed -variable will be no further than one standard deviation from its mean, -that is, let us integrate:: +** Integration under a normal curve - exp (-0.5 * x * x) / sqrt (2.0 * pi) +Simple integration is a particular case of solving ordinary differential equations. Let us compute the probability that a normally distributed variable will be no further than one standard deviation from its mean, that is, let us integrate $e^{x^2\over 2}\over\sqrt{2\pi}$ between $-1$ and $+1$. This is incredibly simple. Note that the problem function does not even use the $v[0]$ variable. -between -1.0 and 1.0. This is incredibly simple. Note that the problem -function does not even use the v[0] variable. +** Rediscovering $cos$ and $sin$ -Rediscovering cos and sin ---------------------------------- +Let us now seek two functions ($cos$ and $sin$), the first being the derivative of the second, the second being the negative of the derivative of the first. The first function has value 1 at point 0, the second function has value 0 at point 0. We want the value of the first function at 1.5 (then: $cos 1.5$). This is quite simple. Note that the problem function does not even use the $t$ variable. Isn't that beautiful that such a straight definition is able to sweep a circle? π is hidden somewhere, there… -Let us now seek two functions (cos and sin), the first being the -derivative of the second, the second being the negative of the derivative -of the first. The first function has value 1.0 at point 0.0, the second -function has value 0.0 at point 0.0. We want the value of the first -function at 1.5 (then: cos 1.5). This is quite simple. Note that the -problem function does not even use the t variable. Isn't that beautiful -that such a straight definition is able to sweep a circle? π is hidden -somewhere, there... +** Box slowing by friction in air -Box slowing by friction in air ------------------------------- +As an example of ordinary differential equations of level higher than one, let have some box of unitary mass moving in the air at speed 100, to which a friction force is applied, equal to $-0.01$ times the square of the speed. We want to know the run distance after 5 units of time. The two functions to solve would be the total distance and the speed, for which derivatives would be the speed and the acceleration. The fact that the speed appears in two places link the functions and increase the level of the differential equations, and this is the main trick. Note that it is not even required to use $t$ in the problem function. -As an example of ordinary differential equations of level higher than -one, let have some box of unitary mass moving in the air at speed 100.0, -to which a friction force is applied, equal to -0.01 times the square of -the speed. We want to know the run distance after 5.0 units of time. -The two functions to solve would be the total distance and the speed, -for which derivatives would be the speed and the acceleration. The fact -that the speed appears in two places link the functions and increase the -level of the differential equations, and this is the main trick. Note -that it is not even required to use t in the problem function. +# Local Variables: +# mode: org +# End:
12 rke.c
 @@ -17,19 +17,11 @@ Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ - - /* Global declarations, initialization and termination. */ - #include "rke.h" - extern double fabs (); -extern char *malloc (); -extern char *alloca (); - - /* Initialize a new system of equations. */ @@ -55,8 +47,6 @@ rke_init (number, routine) /* Newly allocated reentrancy block */ return var; } - - /* Terminate a set of equations. */ void @@ -69,8 +59,6 @@ rke_term (var) /* Main routine of the module, ODE solver. */ - - /* Perform a consistent move of time in the system. */ int /* !0 if success */
3  rke.h
 @@ -17,7 +17,7 @@ Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ - +#include typedef struct struct_rke_variables { @@ -39,7 +39,6 @@ typedef struct struct_rke_variables } *rke_variables; - extern rke_variables rke_init (int, int (*) ()); extern void rke_term (rke_variables); extern int rke_solve (rke_variables, double *, double[], double);

### No commit comments for this range

Something went wrong with that request. Please try again.