Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

Initial

  • Loading branch information...
commit 6ca6d98414bd8e797f3f8a1d7e7f6de8ab5b65c3 0 parents
@pinard authored
Showing with 797 additions and 0 deletions.
  1. +32 −0 Makefile
  2. +212 −0 README
  3. +209 −0 example.c
  4. +299 −0 rke.c
  5. +45 −0 rke.h
32 Makefile
@@ -0,0 +1,32 @@
+# Ordinary differential equation solver, Runge-Kutta-England technique.
+# Copyright © 1988 Free Software Foundation, Inc.
+# François Pinard <pinard@iro.umontreal.ca>, 1988.
+
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 1, or (at your option)
+# any later version.
+
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software Foundation,
+# Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+
+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
212 README
@@ -0,0 +1,212 @@
+=========================================
+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
+----------------------
+
+A reference on Runge-Kutta integration might be:
+
+ Moursund, David G.; Duris, Charles S.: `Elementary Theory and
+ Application of Numerical Analysis`, McGram-Hill, 1967.
+
+The main part of the algorithm has been taken from an appendix in:
+
+ Pritsker, A. Alan B.: `The GASP IV Simulation Language`, Jonh
+ Wiley & Sons, 1974.
+
+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.
+
+Routines of the package
+=======================
+
+See the reentrency variables in `rke.h`_.
+
+See the individual routines headers in `rke.c`_.
+
+Some example programs
+=====================
+
+Here is a small collection of examples to the usage of rke, given here to
+supplement the documentation given previously.
+
+The file `example.c`_ contains the exact code to solve each of these.
+
+Integration under a normal curve
+--------------------------------
+
+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::
+
+ exp (-0.5 * x * x) / sqrt (2.0 * pi)
+
+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``
+---------------------------------
+
+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
+------------------------------
+
+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.
209 example.c
@@ -0,0 +1,209 @@
+/* Ordinary differential equation solver, Runge-Kutta-England technique.
+ Copyright © 1988 Free Software Foundation, Inc.
+ François Pinard <pinard@iro.umontreal.ca>, 1988.
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 1, or (at your option)
+ any later version.
+
+ This program is distributed in the hope that it will be useful, but
+ WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software Foundation,
+ Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+*/
+
+
+#include "rke.h"
+
+
+/* Check how close we can get back to our initial conditions. */
+
+print_return (back, initial)
+ double back; /* Value from solving forth and back */
+ double initial; /* True value before the whole process */
+{
+ printf (" returning to %12.6lf, got %12.6lf\n", initial, back);
+}
+
+
+
+/* Print statistics about number of steps. */
+
+print_steps (var)
+ rke_variables var;
+{
+ printf (" using %3d accepted and %3d rejected steps\n",
+ var->accepted_steps, var->rejected_steps);
+}
+
+
+/* Integration under a normal curve. */
+
+#include <math.h>
+
+static double example_1_const; /* 1.0 / sqrt (2 * pi) */
+
+
+static int
+problem_function_1 (t, v, d)
+ double t; /* Abcissa point */
+ double v[1]; /* v[0] is the surface under the curve */
+ double d[1]; /* d[0] is the normal curve itself */
+{
+ d[0] = example_1_const * exp (-0.5 * t * t);
+ return 1;
+}
+
+
+static
+example_1 ()
+{
+ rke_variables p;
+ double t;
+ double v[1];
+
+ example_1_const = 1.0 / sqrt (2 * 3.1415926);
+
+ p = rke_init (1, problem_function_1);
+
+ t = -1.0; /* Start at -1.0 */
+ v[0] = 0.0; /* Surface is 0.0 at this point */
+
+ /* Now, simply move to +1.0, and collect the answer. */
+
+ if (rke_solve (p, &t, v, 1.0))
+ printf ("\nProbability = %12.6lf.\n", v[0]);
+ else
+ printf ("\nProbability not computed, error.\n");
+ print_steps (p);
+
+ /* Just undo this, to see if we get back where we started. */
+
+ if (rke_solve (p, &t, v, -1.0))
+ print_return (v[0], 0.0);
+ else
+ printf (" return to start not computed, error.\n");
+ print_steps (p);
+
+ rke_term (p);
+}
+
+
+
+/* Rediscovering cos and sin. */
+
+static int
+problem_function_2 (t, v, d)
+ double t;
+ double v[2]; /* v[0] is cos t */
+ /* v[1] is sin t */
+ double d[2]; /* d[0] is d cos t / dt == - sin t */
+ /* d[1] is d sin t / dt == cos t */
+{
+ d[0] = -v[1];
+ d[1] = v[0];
+ return 1;
+}
+
+
+static
+example_2 ()
+{
+ rke_variables p;
+ double t;
+ double v[2];
+
+ p = rke_init (2, problem_function_2);
+
+ t = 0.0; /* Start where we know the values */
+ v[0] = 1.0; /* cos 0 = 1.0 */
+ v[1] = 0.0; /* sin 0 = 0.0 */
+
+ /* Now, simply move to 1.5, and collect the answer. */
+
+ if (rke_solve (p, &t, v, 1.5))
+ printf ("\ncos (1.5) = %12.6lf.\n", v[0]);
+ else
+ printf ("\ncos (1.5) not computed, error.\n");
+ print_steps (p);
+
+ /* Just undo this, to see if we get back where we started. */
+
+ if (rke_solve (p, &t, v, 0.0))
+ {
+ print_return (v[0], 1.0);
+ print_return (v[1], 0.0);
+ }
+ else
+ printf (" return to start not computed, error.\n");
+ print_steps (p);
+
+ rke_term (p);
+}
+
+
+/* Box slowing by friction in air. */
+
+static int
+problem_function_3 (t, v, d)
+ double t; /* Time so far */
+ double v[2]; /* v[0] is the distance so far */
+ /* v[1] is the current speed */
+ double d[2]; /* d[0] is also the current speed */
+ /* d[1] is the current acceleration */
+{
+ d[0] = v[1];
+ d[1] = -0.01 * v[1] * v[1];
+ return 1;
+}
+
+
+static
+example_3 ()
+{
+ rke_variables p;
+ double t;
+ double v[2];
+
+ p = rke_init (2, problem_function_3);
+
+ t = 0.0; /* Start the clock... */
+ v[0] = 0.0; /* ... with no distance so far */
+ v[1] = 100.0; /* ... but some initial speed */
+
+ /* Now, simply ask the clock to be 5.0, and collect the answer. */
+
+ if (rke_solve (p, &t, v, 5.0))
+ printf ("\nDistance = %12.6lf.\n", v[0]);
+ else
+ printf ("\nDistance not computed, error.\n");
+ print_steps (p);
+
+ /* Just undo this, to see if we get back where we started. */
+
+ if (rke_solve (p, &t, v, 0.0))
+ {
+ print_return (v[0], 0.0);
+ print_return (v[1], 100.0);
+ }
+ else
+ printf (" return to start not computed, error.\n");
+ print_steps (p);
+
+ rke_term (p);
+}
+
+
+/* Main program. */
+
+main ()
+{
+ example_1 ();
+ example_2 ();
+ example_3 ();
+}
299 rke.c
@@ -0,0 +1,299 @@
+/* Ordinary differential equation solver, Runge-Kutta-England technique.
+ Copyright © 1988 Free Software Foundation, Inc.
+ François Pinard <pinard@iro.umontreal.ca>, 1988.
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 1, or (at your option)
+ any later version.
+
+ This program is distributed in the hope that it will be useful, but
+ WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software Foundation,
+ 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. */
+
+rke_variables
+rke_init (number, routine) /* Newly allocated reentrancy block */
+ int number; /* Number of simultaneous equations */
+ int (*routine) (); /* Problem function routine */
+{
+ rke_variables var;
+
+ var = (rke_variables) malloc (sizeof (struct struct_rke_variables));
+
+ var->n_equations = number;
+ var->eval_routine = routine;
+ var->minimum_step = 0.000001;
+ var->maximum_step = 1000000.0;
+ var->current_step = 1.0;
+ var->error_slope = 0.0000001;
+ var->error_biais = 0.00000001;
+ var->accepted_steps = 0;
+ var->rejected_steps = 0;
+
+ return var;
+}
+
+
+
+/* Terminate a set of equations. */
+
+void
+rke_term (var)
+ rke_variables var; /* Reentrency block */
+{
+ free (var);
+}
+
+
+/* Main routine of the module, ODE solver. */
+
+
+
+/* Perform a consistent move of time in the system. */
+
+int /* !0 if success */
+rke_solve (var, time, variables, aimed_time)
+ rke_variables var; /* Reentrency block */
+ double *time; /* Current value of time */
+ double variables[]; /* Current variables */
+ double aimed_time; /* Value of time which is aimed for */
+{
+ double whole_step; /* Signed integration step size */
+ double quarter_step; /* 0.25 * whole_step */
+ double half_step; /* 0.50 * whole_step */
+ double three_quarter_step; /* 0.75 * whole_step */
+ double estimated_error; /* Error as estimated by England method */
+ double allowable_error; /* Maximum error that user permits */
+ int within_tolerance; /* Allowable error has not been passed */
+ int all_errors_small; /* All errors within 2% of tolerances */
+ int length_of_array; /* Length of temporary arrays, is bytes */
+ int k; /* Index in various arrays */
+
+ double *dp, *vt, *v, *d;
+ double *a1, *a2, *a3, *a4, *a5, *a6, *a7;
+
+ /* Allocate the work arrays. */
+
+ length_of_array = var->n_equations * sizeof (double);
+
+ dp = (double *) alloca (length_of_array);
+ vt = (double *) alloca (length_of_array);
+ v = (double *) alloca (length_of_array);
+ d = (double *) alloca (length_of_array);
+ a1 = (double *) alloca (length_of_array);
+ a2 = (double *) alloca (length_of_array);
+ a3 = (double *) alloca (length_of_array);
+ a4 = (double *) alloca (length_of_array);
+ a5 = (double *) alloca (length_of_array);
+ a6 = (double *) alloca (length_of_array);
+ a7 = (double *) alloca (length_of_array);
+
+ /* The integration will continue if a minimum step could bring the
+ system closer to the time that is aimed for, even if we have to
+ overshoot it a little. */
+
+ while (2 * fabs (aimed_time - *time) > var->minimum_step)
+ {
+
+ /* Evaluate initial step size and direction. */
+
+ if ((whole_step = aimed_time - *time) > 0.0)
+ {
+ if (whole_step > var->current_step)
+ whole_step = var->current_step;
+ }
+ else
+ {
+ if (whole_step < - var->current_step)
+ whole_step = - var->current_step;
+ }
+
+ /* Evaluate initial differentials. */
+
+ if (! (*var->eval_routine) (*time, variables, dp))
+ return 0;
+
+ do
+
+ /* Loop integrating at this time point until integration error is
+ within tolerances. In any case, adjust integration step size. */
+
+ {
+ /* Calculate various step sizes. */
+
+ quarter_step = 0.25 * whole_step;
+ half_step = quarter_step + quarter_step;
+ three_quarter_step = half_step + quarter_step;
+
+ /* Perform a partial computation for one step of Runge-Kutta
+ 4th order integration, as far as necessary to chain it to
+ England method for estimating integration errors. */
+
+ for (k = 0; k < var->n_equations; ++k)
+ {
+ a1[k] = half_step * dp[k];
+ v[k] = variables[k]
+ + 0.5*a1[k];
+ }
+
+ if (! (*var->eval_routine) (*time + quarter_step, v, d))
+ return 0;
+
+ for (k = 0; k < var->n_equations; ++k)
+ {
+ a2[k] = half_step * d[k];
+ v[k] = variables[k]
+ + 0.25 * (a1[k] + a2[k]);
+ }
+
+ if (! (*var->eval_routine) (*time + quarter_step, v, d))
+ return 0;
+
+ for (k = 0; k < var->n_equations; ++k)
+ {
+ a3[k] = half_step * d[k];
+ v[k] = variables[k]
+ + (-a2[k] + a3[k] + a3[k]);
+ }
+
+ if (! (*var->eval_routine) (*time + half_step, v, d))
+ return 0;
+
+ for (k = 0; k < var->n_equations; ++k)
+ {
+ a4[k] = half_step * d[k];
+ vt[k] = variables[k]
+ + (a1[k] + 4.0*a3[k] + a4[k]) / 6.0;
+ }
+
+ if (! (*var->eval_routine) (*time + half_step, vt, d))
+ return 0;
+
+ for (k = 0; k < var->n_equations; ++k)
+ {
+ a5[k] = half_step * d[k];
+ v[k] = vt[k]
+ + 0.5*a5[k];
+ }
+
+ if (! (*var->eval_routine) (*time + three_quarter_step, v, d))
+ return 0;
+
+ for (k = 0; k < var->n_equations; ++k)
+ {
+ a6[k] = half_step * d[k];
+ v[k] = vt[k]
+ + 0.25*(a5[k] + a6[k]);
+ }
+
+ if (! (*var->eval_routine) (*time + three_quarter_step, v, d))
+ return 0;
+
+ for (k = 0; k < var->n_equations; ++k)
+ {
+ a7[k] = half_step * d[k];
+ v[k] = variables[k]
+ + (-a1[k] - 96.0*a2[k] + 92.0*a3[k] - 121.0*a4[k]
+ + 144.0*a5[k] + 6.0*a6[k] - 12.0*a7[k]) / 6.0;
+ }
+
+ /* Perform England error analysis on partial integration. */
+
+ if (! (*var->eval_routine) (*time + whole_step, v, d))
+ return 0;
+
+ within_tolerance = 1;
+ all_errors_small = 1;
+
+ for (k = 0; k < var->n_equations; ++k)
+ {
+ estimated_error
+ = fabs ((-a1[k] + 4.0*a3[k] + 17.0*a4[k]
+ - 23.0*a5[k] + 4.0*a7[k] - half_step*d[k])
+ / 90.0);
+ allowable_error = fabs (whole_step)
+ * (var->error_slope*fabs (vt[k]) + var->error_biais);
+ if (estimated_error > allowable_error)
+ {
+ within_tolerance = 0;
+ break;
+ }
+ else if (estimated_error > 0.02 * allowable_error)
+ all_errors_small = 0;
+ }
+ if (within_tolerance)
+ {
+ ++var->accepted_steps;
+
+ /* Complete the Runge-Kutta step and return values. */
+
+ for (k = 0; k < var->n_equations; ++k)
+ v[k] = vt[k] + (-a6[k] + a7[k] + a7[k]);
+
+ if (! (*var->eval_routine) (*time + whole_step, v, d))
+ return 0;
+
+ *time += whole_step;
+
+ for (k = 0; k < var->n_equations; ++k)
+ variables[k] = vt[k]
+ + (a5[k] + 4.0*a7[k] + half_step*d[k]) / 6.0;
+
+ /* Increment step size if desirable. */
+
+ if (all_errors_small
+ && fabs (whole_step) == var->current_step)
+ if (2 * var->current_step > var->maximum_step)
+ var->current_step = var->maximum_step;
+ else
+ var->current_step *= 2;
+ }
+ else
+ {
+
+ ++var->rejected_steps;
+
+ /* Decrement step size if possible. */
+
+ if (fabs (whole_step) > var->minimum_step)
+ {
+ if (var->current_step < 2 * var->minimum_step)
+ var->current_step = var->minimum_step;
+ else
+ var->current_step *= 0.5;
+ if (aimed_time > *time)
+ whole_step = var->current_step;
+ else
+ whole_step = - var->current_step;
+ }
+ else
+ return 0; /* Convergence failed */
+ }
+ }
+
+ while (!within_tolerance);
+ }
+ return 1; /* Convergence succeeded */
+}
45 rke.h
@@ -0,0 +1,45 @@
+/* Ordinary differential equation solver, Runge-Kutta-England technique.
+ Copyright © 1988 Free Software Foundation, Inc.
+ François Pinard <pinard@iro.umontreal.ca>, 1988.
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 1, or (at your option)
+ any later version.
+
+ This program is distributed in the hope that it will be useful, but
+ WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software Foundation,
+ Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+*/
+
+
+
+typedef struct struct_rke_variables
+ {
+
+ /* The following are saved from rke_init call arguments. */
+
+ int n_equations; /* Number of simultaneous equations */
+ int (*eval_routine) (); /* Routine to compute derivatives */
+
+ /* These may be changed by the user between two solve calls. */
+
+ double minimum_step; /* Minimum allowable step size */
+ double maximum_step; /* Maximum allowable step size */
+ double current_step; /* Current integration step size */
+ double error_slope; /* Slope of maximum error per time unit */
+ double error_biais; /* Biais of maximum error per time unit */
+ int accepted_steps; /* Accumulated number of accepted steps */
+ int rejected_steps; /* Accumulated number of rejected steps */
+ }
+*rke_variables;
+
+
+extern rke_variables rke_init (int, int (*) ());
+extern void rke_term (rke_variables);
+extern int rke_solve (rke_variables, double *, double[], double);
Please sign in to comment.
Something went wrong with that request. Please try again.