This is a port to JavaScript (actually, TypeScript) of E. Hairer and G. Wanner's implementation of the Gragg-Bulirsch-Stoer method of integrating systems of differential equations. The original code is written in idiomatic Fortran; this code tries to present an idiomatic JavaScript interface while preserving all the of the virtues of the original code, including its speed, configurability, and compact memory footprint.

(We'll write the usage examples in plain JavaScript. The usage from TypeScript is very similar.) ##### One first-order equation

The simplest possible example would be y′ = y, with y(0) = 1: we expect the solution y(x) = exp(x). First we create a solver object, telling how many independent variables there are in the system (in this case just one).

`var odex = ;var s = 1;`

To represent the differential equation, we write a routine that computes y′ given y at the point x. For this example it's very simple:

`var { return y;}`

Since we asked for one independent variable, `y`

is an array of length 1.
We return an array of the same size.

We can solve the equation by supplying the initial data and the start and endpoints. Let's find y(1):

`s; // final x value // { y: [ 2.7182817799042955 ], // outcome: 0, // nStep: 7, // xEnd: 1, // nAccept: 7, // nReject: 0, // nEval: 75 } `

Not bad: the answer `y[1]`

is close to *e*. It would be closer if we requested
more precision. When you create a new `Solver`

object, it is
equipped with a number of properties you can change to control the integration.
You can change a property and re-run the solution:

`sabsoluteTolerance = srelativeTolerance = 1e-10;sy// [ 2.7182818284562535 ] Math - 27182818284562535// 2.7915447731174936e-12 `

You can supply a callback function that runs during the integration to supply intermediate
points of the integration as it proceeds. The callback function is an optional
parameter to `solve`

, which receives the step number, x0, x1 and y(x1). (x0
and x1 represent the interval covered in this integration step).

`sy// 1 0 0 [ 1 ] // 2 0 0.0001 [ 1.0001000050001667 ] // 3 0.0001 0.0007841772783189289 [ 1.000784484825706 ] // 4 0.0007841772783189289 0.004832938716978181 [ 1.0048446362021166 ] // 5 0.004832938716978181 0.01913478583589434 [ 1.0193190291261103 ] // 6 0.01913478583589434 0.0937117110731088 [ 1.0982430889534374 ] // 7 0.0937117110731088 0.2862232977213724 [ 1.3313897183518322 ] // 8 0.2862232977213724 0.7103628434248046 [ 2.034729412908106 ] // 9 0.7103628434248046 1 [ 2.7182818284562535 ] // [ 2.7182818284562535 ] `

You will observe that `odex`

has chosen its own grid points for evaluation.
Adaptive step size is one of the nicest features of this library: you don't
have to worry about it too much.

However, you will often want to sample the data at points of your own choosing.
When you request `denseOutput`

in the `Solver`

parameters, the function you
supply to solve receives a fifth argument which is a closure which you can call to obtain
very accurate y values in the interval [x0, x1]. You call this closure with
the index (within the y vector) of the component you want to evaluate, and the
x value in [x0, x1] where you want to find that y value. One common use case
for this is to obtain otuput at evenly spaced points. To this end, we supply a
canned callback `grid`

which you can use for this:

`sdenseOutput = true; // request interpolation closure in solution callback s;// 0 [ 1 ] // 0.2 [ 1.2214027470178732 ] // 0.4 [ 1.4918240050068732 ] // 0.6 [ 1.8221161568592519 ] // 0.8 [ 2.2255378426172316 ] // 1 [ 2.7182804587510203 ] // [ 2.7182804587510203 ] `

To see how you could use the dense output feature yourself, take a look at
the source to grid.
##### A system of two first order equations
Note that in all these examples, `y`

is a vector: this software is designed to
solve systems. Let's work with the Lotka-Volterra predator-prey system.
The system is:

`dx/dt = a x - b x ydy/dt = c x y - d y`

For odex, we rename *t* to *x*, and then *x* and *y* become `y[0]`

and `y[1]`

.
We write a function LV which binds the constants of the population system
`a`

, `b`

, `c`

, `d`

and returns a function suitable for the integrator.
To represent this system we can write:

`var { return { return a * y0 - b * y0 * y1 c * y0 * y1 - d * y1 ; };};`

Then we can solve it. It's the same as the previous examples, but this time
we need a solver created to handle two independent variables and must supply
initial data for both of them. To find the state of the rabbits and wolves
at time 6, if the state at time zero is {y_{0} = 1, y_{1}
= 1}:

`s = 2;sy// [ 1.6542774481418214, 0.3252864486771545 ] `

To see more of this system of equations in action, you can visit a demo page which allows you to vary the initial conditions with the mouse.

You can integrate second order ordinary differential equations by making a simple transformation to a system of first order equations. Consider Airy's equation: y″ − x y = 0:

In ODEX, we could write y_{0} for y and y_{1} for y′,
so that y″ = y′_{1} and rewrite the system like this:
y′_{0} = y_{1};
y′_{1} − x y_{0} = 0 to get:

`var { return y1 x * y0;}`

There's also a demo page for this equation too.

You might also enjoy a demo of the Lorenz attractor or Van der Pol equation!

This project comes with a mocha test suite. The suite contains other examples of second-order equations which have been translated to systems of first order equations you may examine.