Cadabra manual: Ndsolve

Damped harmonic oscillator

This is the harmonic oscillator with friction example of the Boost ODEint manual.

from cdb.graphics.plot import *

\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].$

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

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

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

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

p=None out=None

for gamma in np.linspace(0, -0.3, 50): sol = ndsolve($[ \dot{x} = y, \dot{y} = -x + @(gamma) y ]$, {$x$:0, $y$:1}, ($\tau$, 0, 10)) p = plot(sol, {$\tau$: (0,10)}, fig = p) out = display(p, out)

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.73\,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]);

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

You can pass conditions into ndsolve, which determine when you want the numerical evolution to stop.

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}\,$

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

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

Cadabra's plot will automatically try to figure out the domain of the function you are trying to plot. If you plot an interpolating function (which is what ndsolve returns), then the domain can be deduced and you do not need to pass it into plot:

plot(sol2);