Numerical ODE Integration¶

Solve ordinary differential equations with np_odeint, which wraps Fortran ODEPACK/DLSODE. Supports non-stiff (Adams) and stiff (BDF) methods.

Two calling conventions:

  • Expression mode — write the ODE symbolically: np_odeint([f1, f2, ...], [t, y1, y2, ...], y0, tspan)
  • Function mode — pass a callback: np_odeint(f, y0, tspan) where f(t, y) takes scalar t and ndarray y
In [15]:
load("numerics")$
load("numerics-integrate")$
load("ax-plots")$

1. Exponential Decay¶

The simplest ODE: $\dot{y} = -y$, $y(0) = 1$. Exact solution $y(t) = e^{-t}$.

In [16]:
tspan : np_linspace(0.0, 5.0, 50)$
result : np_odeint([-y], [t, y], [1.0], tspan)$

/* Extract columns and compare with exact solution */
t_col : np_slice(result, all, 0)$
y_col : np_slice(result, all, 1)$
y_exact : np_exp(np_scale(-1.0, t_col))$
max_err : np_max(np_abs(np_sub(y_col, y_exact)))$
print("Max error vs exp(-t):", max_err)$

Max error vs exp(-t): 1.1436221038912953e-7
In [17]:
ax_draw2d(
  color="blue", line_width=2, name="np_odeint",
  lines(np_to_list(t_col), np_to_list(y_col)),
  color="red", dash="dash", name="exp(-t)",
  explicit(exp(-t), t, 0, 5),
  title="Exponential Decay", xlabel="t", ylabel="y",
  grid=true, showlegend=true
)$
No description has been provided for this image

2. Harmonic Oscillator¶

A 2D system: $\dot{x} = v$, $\dot{v} = -x$, with $(x_0, v_0) = (1, 0)$. Exact solution: $x(t) = \cos t$, $v(t) = -\sin t$.

In [18]:
tspan2 : np_linspace(0.0, float(4*%pi), 200)$
res2 : np_odeint([v, -x], [t, x, v], [1.0, 0.0], tspan2)$

t2 : np_slice(res2, all, 0)$
x2 : np_slice(res2, all, 1)$
v2 : np_slice(res2, all, 2)$
In [19]:
ax_draw2d(
  color="blue", line_width=2, name="x(t) = cos t",
  lines(np_to_list(t2), np_to_list(x2)),
  color="orange", line_width=2, name="v(t) = -sin t",
  lines(np_to_list(t2), np_to_list(v2)),
  title="Harmonic Oscillator", xlabel="t", ylabel="",
  grid=true, showlegend=true
)$
No description has been provided for this image
In [20]:
/* Phase portrait: unit circle */
ax_draw2d(
  color="blue", line_width=2,
  lines(np_to_list(x2), np_to_list(v2)),
  title="Phase Portrait (unit circle)", xlabel="x", ylabel="v",
  grid=true
)$
No description has been provided for this image

3. Lotka-Volterra Predator-Prey¶

$$\dot{x} = \alpha x - \beta x y, \qquad \dot{y} = \delta x y - \gamma y$$

Classic ecological oscillations with closed orbits in the phase plane.

In [21]:
/* Parameters: alpha=1.5, beta=1, delta=1, gamma=3 */
tspan3 : np_linspace(0.0, 15.0, 500)$
res3 : np_odeint(
  [1.5*prey - prey*pred, prey*pred - 3*pred],
  [t, prey, pred], [10.0, 5.0], tspan3)$

t3 : np_slice(res3, all, 0)$
prey3 : np_slice(res3, all, 1)$
pred3 : np_slice(res3, all, 2)$
In [22]:
ax_draw2d(
  color="steelblue", line_width=2, name="Prey",
  lines(np_to_list(t3), np_to_list(prey3)),
  color="tomato", line_width=2, name="Predator",
  lines(np_to_list(t3), np_to_list(pred3)),
  title="Lotka-Volterra", xlabel="t", ylabel="Population",
  grid=true, showlegend=true
)$
No description has been provided for this image
In [23]:
/* Phase plane: closed orbits */
ax_draw2d(
  color="purple", line_width=2,
  lines(np_to_list(prey3), np_to_list(pred3)),
  title="Lotka-Volterra Phase Plane",
  xlabel="Prey", ylabel="Predator", grid=true
)$
No description has been provided for this image

4. Stiff System: Van der Pol Oscillator (BDF)¶

$$\dot{x} = v, \qquad \dot{v} = \mu(1 - x^2) v - x$$

With $\mu = 10$ the system is moderately stiff — sharp transitions between slow drift and fast jumps. We use method=bdf to handle the stiffness.

In [24]:
tspan4 : np_linspace(0.0, 40.0, 2000)$
res4 : np_odeint(
  [vv, 10*(1 - xx^2)*vv - xx],
  [t, xx, vv], [2.0, 0.0], tspan4, 1e-8, 1e-8, bdf)$

t4 : np_slice(res4, all, 0)$
x4 : np_slice(res4, all, 1)$
v4 : np_slice(res4, all, 2)$
In [25]:
ax_draw2d(
  color="blue", line_width=1,
  lines(np_to_list(t4), np_to_list(x4)),
  title="Van der Pol (mu=10, BDF)", xlabel="t", ylabel="x",
  grid=true
)$
No description has been provided for this image
In [26]:
/* Relaxation oscillation in the phase plane */
ax_draw2d(
  color="blue", line_width=1,
  lines(np_to_list(x4), np_to_list(v4)),
  title="Van der Pol Phase Plane", xlabel="x", ylabel="v",
  grid=true
)$
No description has been provided for this image

5. Function Mode: Lorenz Attractor¶

Function mode passes a callback f(t, y) instead of symbolic expressions. The callback receives scalar t and a 1D ndarray y, and returns a Maxima list of derivatives.

$$\dot{x} = \sigma(y - x), \quad \dot{y} = x(\rho - z) - y, \quad \dot{z} = xy - \beta z$$

Standard chaos parameters: $\sigma=10$, $\rho=28$, $\beta=8/3$.

In [27]:
lorenz(t_lr, y_lr) := block(
  [lx, ly, lz, sigma_lr, rho_lr, beta_lr],
  sigma_lr : 10,  rho_lr : 28,  beta_lr : 8/3,
  lx : np_ref(y_lr, 0),
  ly : np_ref(y_lr, 1),
  lz : np_ref(y_lr, 2),
  [sigma_lr*(ly - lx),
   lx*(rho_lr - lz) - ly,
   lx*ly - beta_lr*lz])$

tspan5 : np_linspace(0.0, 40.0, 4000)$
res5 : np_odeint(lorenz, [1.0, 1.0, 1.0], tspan5)$

lx5 : np_slice(res5, all, 1)$
ly5 : np_slice(res5, all, 2)$
lz5 : np_slice(res5, all, 3)$
In [28]:
/* x-z projection of the Lorenz attractor */
ax_draw2d(
  color="navy", marker_size=1,
  points(np_to_list(lx5), np_to_list(lz5)),
  title="Lorenz Attractor (x-z projection)",
  xlabel="x", ylabel="z", grid=true
)$
No description has been provided for this image

Summary¶

Example Mode Method Key feature
Exponential decay Expression Adams Simplest ODE
Harmonic oscillator Expression Adams 2D system, phase portrait
Lotka-Volterra Expression Adams Nonlinear oscillations
Van der Pol Expression BDF Stiff relaxation oscillations
Lorenz attractor Function Adams Chaotic 3D system

Expression mode compiles the RHS into a native Lisp closure (coerce-float-fun) — the Maxima evaluator is not in the integration loop. Function mode calls through the Maxima evaluator on every DLSODE step, but can use vectorized ndarray operations and procedural logic.