Numerical solution of differential equations

Cadabra has functionality to determine numerical solutions to ordinary differential equations. The equations themselves are written as standard Cadabra expressions, and the result is (one or more) interpolating function objects. These can then be used for further analysis, e.g. plotting. To see how this works, it is best to look at some concrete examples.

Damped coupled harmonic oscillators

The system we look at here contains two functions $x(\tau)$ and $y(\tau)$, which satisfy two differential equations. In order to enter these into Cadabra, we first have to declare these two functions and the variable on which they depend.

\tau::Coordinate; \dot{#}::Derivative(\tau); {x,y}::Depends(\tau);

$\displaystyle{}\text{Property Coordinate attached to }\tau.$

$\displaystyle{}\text{Property Derivative attached to }\backslash\texttt{dot}\{\#\}.$

$\displaystyle{}\text{Property Depends attached to }\left[x, y\right].$

The differential equations are given below. There is one parameter $\gamma$ which describes the damping rate.

ho:= [ \dot{x} = y, \dot{y} = -x + \gamma y ];

$\displaystyle{}\left[\dot{x} = y, \dot{y} = -\,x +\gamma y\right]$

Note how, by virtue of the property that $x$ and $y$ depend on $\tau$, we do not have to write that dependence explicitly anymore (so you do not have to write $x(\tau)$ ). Let us specialise to a particular value of the damping,

substitute(ho, $\gamma->-0.2$);

$\displaystyle{}\left[\dot{x} = y, \dot{y} = -\,x -0.2\,y\right]$

We can now solve this system with a call to ndsolve which sets the initial values of $x$ and $y$, and the range of $\tau$ for which we want to solve for the functions.

sol = ndsolve(ho, {$x$: 0, $y$: 1}, ($\tau$, 0, 10));

$\displaystyle{}\left[x = \square{} (\tau), y = \square{} (\tau)\right]$

The result is a list of two equations, one for $x$ and one for $y$, expressing them both in terms of 'anonymous functions'. These are interpolating functions which can be evaluated numerically. If we want to plot the result, do the following:

from cdb.graphics.plot import * plot(sol, title="Damped coupled oscillators");

Note how we did not have to specify the range of the variable $\tau$ for plotting, as the plot command has automatically figured out the range for which the interpolating functions are defined (if you want, you can of course specify the range by hand).

Lotka-Volterra equations

The following example is a little bit more complicated. It shows how to numerically solve a system of two ODEs, and then adds a slider so you can interact with the value of the constant which appears in the equations.

from cdb.interact.slider import * from cdb.graphics.plot import * fig=None

t::Coordinate. {x, y}::Depends(t). \dot{#}::Derivative(t). K = 0.45

The equations themselves describe a very simple predator-prey system, in which $x$ are the number of prey and $y$ are the number of predators.

lotka_volterra := [ \dot{x} = x*(1-y), \dot{y} = - @(K) y*(1-x) ];

$\displaystyle{}\left[\dot{x} = x \left(1 -\,y\right), \dot{y} = -0.69\,y \left(1 -\,x\right)\right]$

We solve these numerically with some initial conditions, and then feed the solutions of ndsolve straight into plot to visualise them.

solution = ndsolve( lotka_volterra, {$x$: 0.6, $y$: 0.2}, ($t$, 0, 30) ) fig = plot(solution, {$t$: (0,30)}, title=f"Lotka-Volterra with K={K}", fig=fig, yrange=[0, 5]);

By adding a slider for the variable $K$, we can create a dynamical plot which updates automatically when the slider is changed:

slider("K", K, 0.1, 1.5, 0.05);

If you drag the slider above, the numerical solution to the ODEs will be re-computed and then re-plotted.

Stopping conditions

Sometimes you do not know the range of the independent variable, but you only know that you want to stop the numerical integration when a particular condition is satisfied (e.g. when the function reaches a particular value). You can pass conditions like these into ndsolve using the stop parameter.

from cdb.graphics.plot import plot

t::Coordinate; x::Depends(t); \dot{#}::Derivative(t);

$\displaystyle{}\text{Property Coordinate attached to }t.$

$\displaystyle{}\text{Property Depends attached to }x.$

$\displaystyle{}\text{Property Derivative attached to }\backslash\texttt{dot}\{\#\}.$

eq:= \dot{x} = x/5+1/10;

$\displaystyle{}\dot{x} = \frac{1}{5}\,x +\frac{1}{10}\,$

We want to evaluate this until $x$ gets larger than $2$, so we do

sol2 = ndsolve(eq, {$x$: 0}, ($t$, 0, 100), stop=$x > 2$);

$\displaystyle{}\left[x = \square{} (t)\right]$

As above, we can now plot this by simply passing the solution to the plot function, which will figure out the maximal range of $t$ for which the numerical solution has been obtained.

plot(sol2);