Solve ordinary differential equations (initial and boundary value problems). Supports stiff/non-stiff systems, event detection, Hamiltonian/symplectic integration, parameter sweeps, and phase space analysis. Use for any ODE system in physics, engineering, or applied math.
Solve ordinary differential equations (initial and boundary value problems). Supports stiff/non-stiff systems, event detection, Hamiltonian/symplectic integration, parameter sweeps, and phase space analysis. Use for any ODE system in physics, engineering, or applied math.
Solve systems of ordinary differential equations using scipy's battle-tested integrators. Covers initial value problems (IVP), boundary value problems (BVP), event detection, stiff systems, and Hamiltonian dynamics with symplectic integrators.
When to Use
Any initial value problem: dx/dt = f(t, x)
Boundary value problems: solve with constraints at two endpoints
Hamiltonian systems needing energy-conserving integration
Parameter sweeps over ODE systems
Systems with discrete events (bouncing ball, switching dynamics)
Do NOT Use When
Solving PDEs (use pde-solver or fluidsim)
Fitting ODE parameters to data (use physics-fitting, then come back here)
Discovering governing equations from data (use symbolic-regression or sindy)
Core Workflows
1. Basic IVP (Initial Value Problem)
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# Define the system: dy/dt = f(t, y)defharmonic_oscillator(t, y, omega=2.0):
"""Simple harmonic oscillator: x'' + omega^2 * x = 0"""
x, v = y
return [v, -omega**2 * x]
# Solve
t_span = (0, 10)
y0 = [1.0, 0.0] # x(0)=1, v(0)=0
sol = solve_ivp(harmonic_oscillator, t_span, y0,
method='RK45', t_eval=np.linspace(0, 10, 1000),
rtol=1e-10, atol=1e-12)
# Plot
plt.plot(sol.t, sol.y[0], label=)
plt.plot(sol.t, sol.y[], label=)
plt.xlabel()
plt.ylabel()
plt.legend()
plt.grid(, alpha=)
plt.savefig(, dpi=, bbox_inches=)
'x(t)'
1
'v(t)'
'Time [s]'
'State'
True
0.3
'harmonic_oscillator.png'
150
'tight'
2. Stiff Systems
Use method='Radau' or method='BDF' for stiff problems:
defstiff_system(t, y):
"""Van der Pol oscillator (stiff for large mu)"""
mu = 1000# stiffness parameter
x, v = y
return [v, mu * (1 - x**2) * v - x]
sol = solve_ivp(stiff_system, (0, 3000), [2.0, 0.0],
method='Radau', rtol=1e-8, atol=1e-10,
max_step=10.0)
How to detect stiffness:
Explicit solver (RK45) takes extremely many steps or fails
System has widely separated timescales
Jacobian eigenvalues have large negative real parts
3. Event Detection
Find when specific conditions are met during integration:
defprojectile(t, y):
"""Projectile motion with drag"""
x, vx, z, vz = y
g = 9.80665
drag = 0.01# drag coefficient
speed = np.sqrt(vx**2 + vz**2)
return [vx, -drag * speed * vx,
vz, -g - drag * speed * vz]
defhit_ground(t, y):
"""Event: z = 0 (projectile hits ground)"""return y[2]
hit_ground.terminal = True# stop integration
hit_ground.direction = -1# only when z is decreasingdefmax_height(t, y):
"""Event: vz = 0 (apex)"""return y[3]
max_height.direction = -1
sol = solve_ivp(projectile, (0, 100), [0, 50, 0, 50],
events=[hit_ground, max_height],
max_step=0.1, dense_output=True)
print(f"Impact at t = {sol.t_events[0][0]:.3f} s")
print(f"Max height at t = {sol.t_events[1][0]:.3f} s")
4. Hamiltonian Systems (Symplectic Integration)
For energy-conserving systems, use symplectic methods:
