How to Solve ODEs Numerically in Python with SciPy’s odeint

1 Numerical Solutions

1.1 odeint

Python solves ordinary differential equations (ODEs) numerically by converting higher‑order equations to first‑order systems and applying Runge–Kutta methods. The scipy.integrate module provides the odeint function, whose basic call signature is:

<code>odeint(func, y0, t, ...)</code>

where func defines the differential equation (or system) as a function of the state vector and the independent variable, y0 is the sequence of initial conditions, and t is the sequence of time points (the first element must be the initial time). The return value sol contains the numerical solution evaluated at each time point; for a system with multiple equations, sol is a matrix whose columns correspond to the solutions of each variable.

1.2 Example 1

Solve the ODE dy/dx = -2*y + x**2 + 2*x with initial condition y(1)=2 over the interval x∈[1,10] with step 0.5.

<code>from scipy.integrate import odeint  # import integration function
from numpy import arange               # generate arithmetic sequence

dy = lambda y, x: -2*y + x**2 + 2*x   # define the ODE
x = arange(1, 10.5, 0.5)               # time series from 1 to 10, step 0.5
sol = odeint(dy, 2, x)                # compute numerical solution
print("x={}\nCorresponding solution y={}".format(x, sol.T))  # display results</code>

1.3 Example 2

Convert a second‑order ODE to a first‑order system and solve it numerically.

<code>from scipy.integrate import odeint
from sympy.abc import t
import numpy as np

def func(y, x):
    y1, y2 = y
    return np.array([y2, -2*y1 - 2*y2])

x = np.arange(0, 10, 0.1)          # create time points
sol1 = odeint(func, [0.0, 1.0], x)  # compute numerical solution</code>

1.4 Example 3 – Lorenz Chaotic Model

The Lorenz model, a classic deterministic system of three first‑order ODEs, exhibits chaotic behavior. Its equations are:

<code>dx/dt = sigma*(y - x)
dy/dt = rho*x - y - x*z
dz/dt = x*y - beta*z</code>

Typical parameter values are sigma = 10 , rho = 28 , and beta = 8/3 .

<code>from scipy.integrate import odeint
import numpy as np

def lorenz(w, t):
    sigma = 10; rho = 28; beta = 8/3
    x, y, z = w
    return np.array([sigma*(y - x), rho*x - y - x*z, x*y - beta*z])

t = np.arange(0, 50, 0.01)               # time points
sol1 = odeint(lorenz, [0.0, 1.0, 0.0], t)  # first initial condition
sol2 = odeint(lorenz, [0.0, 1.0001, 0.0], t)  # slightly perturbed initial condition</code>

The two trajectories diverge rapidly, illustrating the system’s sensitivity to initial conditions (the “butterfly effect”).

Reference

Python数学实验与建模 / 司守奎, 孙玺菁, 科学出版社