ode45-cash-karp
v1.1.0
Published
Integrate a system of Ordinary Differential Equations using the Fifth Order Adaptive Cash-Karp method
Downloads
28
Readme
ode45-cash-karp
Integrate a system of Ordinary Differential Equations using the Fifth Order Adaptive Cash-Karp method
Introduction
This module integrates a system of ordinary differential equations of the form
where is a vector of length . Given time step , the Cash-Karp method uses a fifth order Runge-Kutta scheme with a fourth order embedded estimator in order to control the error. In other words, the same intermediate values used in calculating the fifth order update can be used to calculate a fourth order estimate. The difference yields an error estimate, and the error estimate controls the timestep .
Install
$ npm install ode45-cash-karp
Example
Van der Pol oscillator
var ode45 = require('ode45-cash-karp')
// The derivative function for a Van der Pol oscillator:
var vanderpol = function(dydt, y, t) {
dydt[0] = y[1]
dydt[1] = 4 * (1-y[0]*y[0])*y[1] - y[0]
}
// Initialize:
var y0 = [2,0],
t0 = 0,
dt0 = 1e-3,
integrator = ode45( y0, vanderpol, t0, dt0 )
// Integrate up to tmax:
var tmax = 10, t = [], y = []
while( integrator.step( tmax ) ) {
// Store the solution at this timestep:
t.push( integrator.t )
y.push( integrator.y )
}
Common patterns
A single adaptive step:
integrator.step()
integrator.y // current state
integrator.t // current time
integrator.dt // newly adapted timestep
A single adaptive step, returning false if reached:
integrator.step( 1.5 ) // returns true if t < 1.5
integrator.step( 1.5 ) // returns false if t = 1.5
Ten adaptive timesteps in sequence:
integrator.steps( 10 )
Ten timesteps, halting early and returning false if is reached::
integrator.steps( 10, 1.5 )
Take any number of timesteps until :
integrator.steps( Infinity, 1.5 )
API
require('ode45-cash-karp')( y0, deriv, t0, dt0 [, options] )
Arguments:
y0
: an array or typed array containing initial conditions. This vector is updated in-place with each integrator step.deriv
: a function that calculates the derivative. Format isfunction( dydt, y, t )
. Inputs are current statey
and current timet
, output is the calculated derivativedydt
.t0
: initial time .dt0
: initial time step .options
: an optional associative array of options. Valid parameters are:tol
: The target error level to be acheived. Default is: .maxIncreaseFactor
: The maximum factor by which to increase the timestep if the error tolerance is met. Default value is 10. This limit is applied at the end of a successful timestep.maxDecreaseFactor
: The maximum factor by which to decrease the timestep if the error tolerance is not met. Default value is 10. This limit is applied on each trial step until the error tolerance is acheived.dtMinMag
: The minimum allowed magnitude of . If limit is exceeded during adaptation, a warning is printed to the console and the timestep completes with clipped to the the prescribed magnitude. Ifundefined
, this limit is ignored. Default value isundefined
.dtMaxMag
: The maximum allowed magnitude of . This limit is applied at the beginning of each step. If a timestep larger than this magnitude is requested, the timestep is executed with clipped to the prescribed magnitude. Ifundefined
, this limit is ignored. Default value isundefined
.errorScaleFunction
: The function used to compute a normalizing factor for the error in a given dimension. See below for details.errorReduceFunction
: The reduce operation by which errors in each dimension are combined into a single error metric. See below for details.errorPostFunction
: An operation applied to the total error. For example, if using the norm this would be a square root. See below for details.verbose
: Log covergence warnings. Default istrue
.
Returns:
Initialized integrator object.
Properties:
n
: dimension ofy0
.y
: current state; a reference to input arrayy0
.deriv
: function that calculates the derivative.t
: current time, incremented bydt
on each time step.dt
: current time step . Initialized from inputdt0
. May be changed, but will be overwritten with each adaptive step in order to acheive the prescribed error bound.- all options are copied to properties on the integrator object and may be changed at any time.
Methods:
.step( [tLimit] )
: takes a single step of the integrator and stores the result in-place in they
property. Returns true iftLimit
was not provided or ift
has not reached the limit, otherwise returns false, meaningt
has reachedtLimit
..steps( n, [tLimit] )
: takesn
steps of the integrator, storing the result in-place in they
property. Exits early iftLimit
is reached. Returns true iftLimit
was not provided or ift
has not reached the limit, otherwise returns false, meaningt
has reachedtLimit
. Note that, for example,.steps( Infinity, 10 )
is valid and will take whatever number of step is required to reach .
Error Estimation
Ideally, there would be no choices in error computation since this library would implement the best possible choices, but I've left this configurable.
errorScaleFunction: function( i, dt, y, dydt )
This function receives the dimension number i
, the current timestep dt
, the current state y
, and the derivative calculated at the beginning of the step, dydt
. It must return a normalization factor by which the error in the given dimension is normalized. It is executed once at the beginning of each timestep and not for subsequent trial steps. By default, it is:
function errorScaleFunction( i, dt, y, dydt ) {
return Math.abs(y) + Math.abs(dt * dydt) + 1e-32
}
errorReduceFunction: function( i, accumulatedError, errorEstimate )
This function performs a reduce operation on the per-dimension error. accumulatedError
is initially zero. The function must add the error estimate in dimension i
and return a new error estimate. By default, the error reduce function simply returns the maximum error:
function errorReduceFunction( i, accumulatedError, errorEstimate ) {
return Math.max( accumulatedError, Math.abs(errorEstimate))
}
errorPostFunction: function( accumulatedError, errorEstimate )
This function applies a mapping to the total reduced error resulting from errorReduceFunction
. For the norm, this would just be Math.sqrt
; for the norm, this is simply a no-op:
function errorPostFunction( accumulatedError ) {
return accumulatedError
}
Credits
(c) 2015 Ricky Reusser. MIT License