Class RungeKutta

Object
Model
ODE
RungeKutta
All Implemented Interfaces:
DModel, CLOProvider
Direct Known Subclasses:
ATBT.ODE, LV.ODE, SIR.ODE
public class RungeKutta extends ODE
Implementation of Runge-Kutta method with adaptive step size for the numerical integration of systems of differential equations. Adapted from numerical recipes in C.
Author:
Christoph Hauert

  • Field Details

    • autoDt

      boolean autoDt
      true if the adaptive step sizes should be used. If false the traditional Euler's method with fixed step size is used.
      See Also:

    • doEuler

      boolean doEuler
      true to resort to the traditional Euler's method with fixed step size. This happens if autoDt==false or in multi-species modules with mixed frequency and density dynamics.
      See Also:

    • SAFETY

      private static final double SAFETY
      The safety margin for adjusting time steps. Magic number from Numerical Recipes in C.
      See Also:

    • PGROW

      private static final double PGROW
      The exponent for increasing the time step if error outside margin. Magic number from Numerical Recipes in C.
      See Also:

    • PSHRNK

      private static final double PSHRNK
      The exponent for decreasing the time step if error outside margin. Magic number from Numerical Recipes in C.
      See Also:

    • ERRCON

      private static final double ERRCON
      The maximum error in a single step. The value ERRCON equals (5/SAFETY)^(1/PGROW). Magic number from Numerical Recipes in C.
      See Also:

    • ACCURACY

      private static final double ACCURACY
      The accuracy required for the integration. Magic number from Numerical Recipes in C.
      See Also:

    • yerr

      private double[] yerr
      Helper variables and temporary storage for errors, states and fitness when calculating derivatives for different steps.

    • ytmp

      private double[] ytmp
      Temporary state vector.

    • ftmp

      private double[] ftmp
      Temporary derivative vector.

    • A2

      private static final double A2
      Cash-Karp coefficient for the second stage time.
      See Also:

    • A3

      private static final double A3
      Cash-Karp coefficient for the third stage time.
      See Also:

    • A4

      private static final double A4
      Cash-Karp coefficient for the fourth stage time.
      See Also:

    • A5

      private static final double A5
      Cash-Karp coefficient for the fifth stage time.
      See Also:

    • A6

      private static final double A6
      Cash-Karp coefficient for the sixth stage time.
      See Also:

    • B21

      private static final double B21
      Cash-Karp coefficient for the second stage derivative.
      See Also:

    • B31

      private static final double B31
      Cash-Karp coefficient for the third stage derivative.
      See Also:

    • B32

      private static final double B32
      Cash-Karp coefficient for the third stage derivative.
      See Also:

    • B41

      private static final double B41
      Cash-Karp coefficient for the fourth stage derivative.
      See Also:

    • B42

      private static final double B42
      Cash-Karp coefficient for the fourth stage derivative.
      See Also:

    • B43

      private static final double B43
      Cash-Karp coefficient for the fourth stage derivative.
      See Also:

    • B51

      private static final double B51
      Cash-Karp coefficient for the fifth stage derivative.
      See Also:

    • B52

      private static final double B52
      Cash-Karp coefficient for the fifth stage derivative.
      See Also:

    • B53

      private static final double B53
      Cash-Karp coefficient for the fifth stage derivative.
      See Also:

    • B54

      private static final double B54
      Cash-Karp coefficient for the fifth stage derivative.
      See Also:

    • B61

      private static final double B61
      Cash-Karp coefficient for the sixth stage derivative.
      See Also:

    • B62

      private static final double B62
      Cash-Karp coefficient for the sixth stage derivative.
      See Also:

    • B63

      private static final double B63
      Cash-Karp coefficient for the sixth stage derivative.
      See Also:

    • B64

      private static final double B64
      Cash-Karp coefficient for the sixth stage derivative.
      See Also:

    • B65

      private static final double B65
      Cash-Karp coefficient for the sixth stage derivative.
      See Also:

    • C1

      private static final double C1
      Cash-Karp fifth-order output coefficient.
      See Also:

    • C3

      private static final double C3
      Cash-Karp fifth-order output coefficient.
      See Also:

    • C4

      private static final double C4
      Cash-Karp fifth-order output coefficient.
      See Also:

    • C6

      private static final double C6
      Cash-Karp fifth-order output coefficient.
      See Also:

    • DC5

      private static final double DC5
      Cash-Karp error estimate coefficient.
      See Also:

    • DC1

      private static final double DC1
      Cash-Karp error estimate coefficient.
      See Also:

    • DC3

      private static final double DC3
      Cash-Karp error estimate coefficient.
      See Also:

    • DC4

      private static final double DC4
      Cash-Karp error estimate coefficient.
      See Also:

    • DC6

      private static final double DC6
      Cash-Karp error estimate coefficient.
      See Also:

    • ak2

      private double[] ak2
      Temporary variables for intermediate results required to implement the fifth-order Cash-Karp Runge-Kutta method.

    • ak3

      private double[] ak3
      Third stage derivative vector.

    • ak4

      private double[] ak4
      Fourth stage derivative vector.

    • ak5

      private double[] ak5
      Fifth stage derivative vector.

    • ak6

      private double[] ak6
      Sixth stage derivative vector.

  • Constructor Details

    • RungeKutta

      public RungeKutta(EvoLudo engine)
      Constructs a new model for the numerical integration of the system of ordinary differential equations representing the dynamics specified by the Module module using the EvoLudo pacemaker engine to control the numerical evaluations. The integrator implements the fifth order Runge-Kutta method with adaptive step size.
      Parameters:
      engine - the pacemaker for running the model

  • Method Details

    • unload

      public void unload()
      Description copied from class: Model
      Milestone: Unload this model and free resources (if applicable).
      Overrides:
      unload in class ODE
      See Also:

    • reset

      public void reset()
      Description copied from class: Model
      Milestone: Reset this model
      Overrides:
      reset in class ODE
      See Also:

    • setAutoDt

      public void setAutoDt(boolean autoDt)
      Sets whether adaptive step sizes should be used. This is required for this implementation of the fifth order Cash-Karp Runge-Kutta method. If set to false the integrator defaults to the Euler method with fixed increments.
      Parameters:
      autoDt - the flag to indicate whether adaptive step sizes should be used.
      See Also:

    • getAutoDt

      public boolean getAutoDt()
      Gets whether adaptive step sizes are used. If not the integrator defaults to the Euler method with fixed increments.
      Returns:
      true if using adaptive step sizes

    • check

      public boolean check()
      Description copied from class: Model
      Check consistency of parameters and adjust if necessary (and possible). All issues and modifications should be reported through logger. Some parameters can be adjusted while the model remains active or even while running, whereas others require a reset. An example of the former category is in general simple adjustments of payoffs, while an example of the latter category is a change of the population structure.
      Overrides:
      check in class ODE
      Returns:
      true if reset required

    • deStep

      protected double deStep(double step)
      Attempts a numerical integration step of size step with step size adjustment. Starts with an Euler step but if the result yields:
      1. negative densities/frequencies
      2. frequencies exceeding 1.0
      the step size is reduced such that densities/frequencies are set to 0.0 or 1.0, respectively, and the effective step size is adjusted and recorded in ODE.dtTaken.

      Implementation Notes:

      Fifth-order Runge-Kutta step with monitoring of local truncation error to ensure accuracy and adjust step size. Input are the dependent variable vector y[1..n] and its derivative dydx[1..n] at the starting value of the independent variable x. Also input are the step size to be attempted htry, the required accuracy eps, and the vector yscal[1..n] against which the error is scaled. On output, y and x are replaced by their new values, hdid is the step size that was actually accomplished, and hnext is the estimated next step size. derivs is the user-supplied routine that computes the right-hand side derivatives.

      If you desire constant fractional errors, plug a pointer to y into the pointer to yscal calling slot (no need to copy the values into a different array). If you desire constant absolute errors relative to some maximum values, set the elements of yscal equal to those maximum values. A useful “trick” for getting constant fractional errors except “very” near zero crossings is to set yscal[i] equal to |y[i]| + |h × dydx[i]|. (The routine odeint, below, does this.)

      Copied from Numerical Recipes in C, chapter 16.2, p.718f

      Overrides:
      deStep in class ODE
      Parameters:
      step - the time step to attempt
      Returns:
      squared distance between this state and previous one, (yt-yout)2.

    • rkck

      private boolean rkck(double h)
      Given values for n variables y[1..n] and their derivatives dydx[1..n] known at x, use the fifth-order Cash-Karp Runge-Kutta method to advance the solution over an interval h and return the incremented variables as yout[1..n]. Also return an estimate of the local truncation error in yout using the embedded fourth-order method. The user supplies the routine derivs(x,y,dydx), which returns derivatives dydx at x.

      Copied from Numerical Recipes in C, chapter 16.2, p.719

      Implementation Notes:

      For replicator dynamics all intermediate steps must remain normalized, i.e. sum ytmp = 1, otherwise reject. Currently this requires that in multi-species modules either all or none are frequency based.
      Parameters:
      h - the step to try to take
      Returns:
      true if step successful