Numerical Solutions of Differential Equations in Python
Department of Mathematics
The Chinese University of Hong Kong
HSMMC Math Modelling Advanced Training Workshop
May 16, 2026
Numerical methods solve differential equations by computing discrete approximate solutions from initial conditions using small iterative steps, with accuracy/efficiency dependent on the algorithm and time step; common methods include Euler’s and Runge‑Kutta.
MATLAB provides ODE solvers: ode45 (general‑purpose Runge‑Kutta) and four stiff solvers (ode15s, ode23s, ode23t, ode2) .
Stiff ODEs involve solutions with widely varying time scales (fast/slow components), causing standard methods to become unstable unless using extremely small time steps; stiffness depends on the numerical method, not just the equation.
To use MATLAB ODE solvers, users define an ODE function, initial conditions, and time span; solvers return time/solution vectors for plotting.
MATLAB’s dsolve is a symbolic (analytical) solver that fails for complex, nonlinear, or non‑closed‑form ODEs, unlike numerical solvers.
The table below provides the direct correspondence between MATLAB ODE solvers and their equivalent Python methods using \(\texttt{scipy.integrate.solve_ivp}\) (numerical) and \(\texttt{sympy.dsolve}\) (symbolic).
\[ \begin{array}{ccc} \hline \textbf{MATLAB Solver} & \textbf{Python} ~\texttt{solve_ivp} ~\textbf{Method} & \textbf{Purpose} \\ \hline \texttt{ode45} & \texttt{RK45} & \textrm{General-purpose (default)} \\ \texttt{ode15s} & \texttt{Radau} & \textrm{Stiff systems} \\ \texttt{ode23s} & \texttt{BDF} & \textrm{Stiff systems} \\ \texttt{ode23t} & \texttt{LSODA} & \textrm{Stiff/non-stiff switch} \\ \texttt{ode23tb} & \texttt{Radau} & \textrm{Stiff systems} \\ \texttt{dsolve} & \texttt{sympy.dsolve} & \textrm{Analytical (symbolic) solution} \\ \hline \end{array} \]
1 Example 1
Employ Python to plot the solution of the initial value problem: \[ y' = 4x, \quad y(0) = 1, \quad x \in [0, 2.5] \]
# --------------------------------------------
# Python code to find ANALYTICAL SOLUTION of
# y' = 4x, y(0) = 1, x ∈ [0, 2.5]
# Uses symbolic math (SymPy) for exact solution
# --------------------------------------------
import sympy as sp
import numpy as np
# Define symbolic variables
x = sp.Symbol('x') # independent variable x
y = sp.Function('y')(x) # dependent variable y(x)
# Define the ODE: y' = 4x
ode = sp.Eq(sp.diff(y, x), 4*x)
# Solve the ODE analytically with initial condition y(0) = 1
analytical_sol = sp.dsolve(ode, y, ics={y.subs(x, 0): 1})
# Print the EXACT analytical solution
print("===== Analytical Solution (Exact) =====")## ===== Analytical Solution (Exact) =====## Eq(y(x), 2*x**2 + 1)# Extract the right-hand side (formula for y(x))
y_exact = analytical_sol.rhs
print("\nExplicit formula: y(x) =", y_exact)##
## Explicit formula: y(x) = 2*x**2 + 1# --------------------------------------------
# Evaluate the analytical solution at any x
# --------------------------------------------
# Convert symbolic formula to a numerical function
y_func = sp.lambdify(x, y_exact, 'numpy')
# Test at specific points (matches your problem)
x_test1 = 1.3125
x_test2 = 0.5625
y_test1 = y_func(x_test1)
y_test2 = y_func(x_test2)
print("\n===== Test Points =====")##
## ===== Test Points =====## At x = 1.3125: y = 4.4453125## At x = 0.5625: y = 1.6328125# Evaluate over the full interval [0, 2.5]
x_vals = np.linspace(0, 2.5, 10)
y_vals = y_func(x_vals)
print("\n===== Analytical Solution over x ∈ [0, 2.5] =====")##
## ===== Analytical Solution over x ∈ [0, 2.5] =====## x = 0.000 y = 1.0000
## x = 0.278 y = 1.1543
## x = 0.556 y = 1.6173
## x = 0.833 y = 2.3889
## x = 1.111 y = 3.4691
## x = 1.389 y = 4.8580
## x = 1.667 y = 6.5556
## x = 1.944 y = 8.5617
## x = 2.222 y = 10.8765
## x = 2.500 y = 13.5000# ---------------------------
# Python Implementation: ODE Solving and Plotting
# Problem: y' = 4x, y(0) = 1, x ∈ [0, 2.5]
# ---------------------------
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# ---------------------------
# 1. Define the ODE function
# ---------------------------
def ode_func(x, y):
"""Right-hand side of the ODE: dy/dx = 4x"""
return 4 * x
# ---------------------------
# 2. Set problem parameters
# ---------------------------
x0 = 0 # Initial x value
y0 = [1] # Initial condition y(x0) = 1
xspan = [0, 2.5] # Solution interval
# ---------------------------
# 3. Solve ODE with fixed evaluation points (stable, no index error)
# ---------------------------
x_eval = np.linspace(0, 2.5, 100) # Define 100 stable points
sol = solve_ivp(
fun=ode_func,
t_span=xspan,
y0=y0,
method='RK45',
t_eval=x_eval # Use fixed points to guarantee enough data
)
# Extract solution results
xsol = sol.t
ysol = sol.y[0]
# ---------------------------
# 4. Analytical solution for comparison
# ---------------------------
def analytical_solution(x):
"""Analytical solution: y(x) = 2x² + 1"""
return 2 * x**2 + 1
y_analytical = analytical_solution(xsol)
# ---------------------------
# 5. Plot results
# ---------------------------
plt.figure(figsize=(8, 5))
plt.plot(xsol, ysol, 'b-', linewidth=1.5, label='Numerical Solution (RK45)')
plt.plot(xsol, y_analytical, 'r--', linewidth=1.5, label='Analytical Solution: $y=2x^2+1$')
plt.grid(True, linestyle='--', alpha=0.7)
plt.xlabel('x')
plt.ylabel('y(x)')
plt.title(r'Solution of $y^\prime = 4x, \ y(0)=1$')
plt.legend()
plt.tight_layout()
# plt.savefig('ode_solution_plot.png', dpi=300, bbox_inches='tight')
# plt.close() # Avoid display freeze in R Markdown
# ---------------------------
# 6. Print results safely (no index error)
# ---------------------------
print("=== Numerical Solution Points (First 10) ===")## === Numerical Solution Points (First 10) ===## x y_numerical y_analytical error## ---------------------------------------------# Use min to avoid index out of bounds
n_print = min(10, len(xsol))
for i in range(n_print):
x = xsol[i]
y_num = ysol[i]
y_an = y_analytical[i]
err = abs(y_num - y_an)
print(f"{x:<10.6f} {y_num:<12.6f} {y_an:<12.6f} {err:.2e}")## 0.000000 1.000000 1.000000 0.00e+00
## 0.025253 1.001275 1.001275 2.22e-16
## 0.050505 1.005102 1.005102 2.22e-16
## 0.075758 1.011478 1.011478 0.00e+00
## 0.101010 1.020406 1.020406 2.22e-16
## 0.126263 1.031885 1.031885 0.00e+00
## 0.151515 1.045914 1.045914 0.00e+00
## 0.176768 1.062494 1.062494 0.00e+00
## 0.202020 1.081624 1.081624 0.00e+00
## 0.227273 1.103306 1.103306 0.00e+00##
## === Step Size Information ===dx = np.diff(xsol)
if len(dx) >= 5:
print(f"First 5 step sizes: {dx[:5]}")
else:
print(f"Step sizes: {dx}")## First 5 step sizes: [0.02525253 0.02525253 0.02525253 0.02525253 0.02525253]## Average step size: 0.025253##
## === Key Test Point Verification ===for x_test in [1.3125, 0.5625]:
y_test = analytical_solution(x_test)
print(f"x = {x_test:.6f}, Analytical y = {y_test:.6f}")## x = 1.312500, Analytical y = 4.445312
## x = 0.562500, Analytical y = 1.632812
2 Example 2
Employ Python to plot the solution of the initial value problem: \[ y' + 3y = x e^{2x}, \quad y(2) = -4, \quad x \in [2, 2.7]. \]
Exact Analytical Solution (from SymPy) \[ y(x) = \left( \frac{x}{5} - \frac{1}{25} \right) e^{2x} + \left( -4 e^{6} - \frac{9}{25} e^{4} \right) e^{-3x}. \]
Simplified Form \[ y(x) = \frac{5x - 1}{25} \, e^{2x} + C e^{-3x}, \] where \[ C = -4 e^{6} - \frac{9}{25} e^{4}. \]
# ---------------------------
# Python Code: Analytical Solution for
# y' + 3y = x*exp(2x), y(2) = -4, x ∈ [2, 2.7]
# ---------------------------
import sympy as sp
import numpy as np
# Define symbolic variables
x = sp.Symbol('x')
y = sp.Function('y')(x)
# Define the ODE: y' + 3y = x*exp(2x)
ode = sp.Eq(sp.diff(y, x) + 3*y, x * sp.exp(2*x))
# Solve the ODE with initial condition y(2) = -4
solution = sp.dsolve(ode, y, ics={y.subs(x, 2): -4})
# Print the exact analytical solution
print("=== Exact Analytical Solution ===")## === Exact Analytical Solution ===## Eq(y(x), x*exp(2*x)/5 - exp(2*x)/25 + (-9*exp(10)/25 - 4*exp(6))*exp(-3*x))##
## Explicit formula:## y(x) = x*exp(2*x)/5 - exp(2*x)/25 + (-9*exp(10)/25 - 4*exp(6))*exp(-3*x)# Convert symbolic solution to a numerical function
y_func = sp.lambdify(x, y_exact, 'numpy')
# Evaluate at key points
x_points = np.array([2, 2.2, 2.5, 2.7])
y_values = y_func(x_points)
print("\n=== Solution at key points ===")##
## === Solution at key points ===## x = 2.00 → y = -4.0000
## x = 2.20 → y = 19.5980
## x = 2.50 → y = 62.9918
## x = 2.70 → y = 107.8065# Optional: Evaluate over the full interval
x_interval = np.linspace(2, 2.7, 10)
y_interval = y_func(x_interval)
print("\n=== Solution over [2, 2.7] ===")##
## === Solution over [2, 2.7] ===## x y(x)## ---------------------------## 2.00 -4.0000
## 2.08 5.2233
## 2.16 14.3129
## 2.23 23.6600
## 2.31 33.6461
## 2.39 44.6591
## 2.47 57.1082
## 2.54 71.4389
## 2.62 88.1493
## 2.70 107.8065# ---------------------------
# Python Implementation for Example 2
# Problem: y' + 3y = x*exp(2x), y(2) = -4, x ∈ [2, 2.7]
# ---------------------------
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import sympy as sp
# ---------------------------
# 1. Define the ODE Right-Hand Side
# ---------------------------
def ode_func(x, y):
"""
Right-hand side of the ODE, rearranged for numerical solving:
y' = -3y + x*exp(2x)
"""
return -3 * y + x * np.exp(2 * x)
# ---------------------------
# 2. Set Problem Parameters
# ---------------------------
x0 = 2 # Initial x value
y0 = [-4] # Initial condition y(2) = -4
xspan = [2, 2.7] # Solution interval
# ---------------------------
# 3. Solve ODE Numerically (RK45)
# ---------------------------
# Use fixed points to guarantee stable output
x_eval = np.linspace(2, 2.7, 100)
sol = solve_ivp(
fun=ode_func,
t_span=xspan,
y0=y0,
method='RK45',
t_eval=x_eval
)
# Safe extraction of results
xsol = sol.t
ysol = sol.y[0] if sol.y.ndim > 1 else sol.y
# ---------------------------
# 4. Analytical Solution (SymPy)
# ---------------------------
x_sym = sp.Symbol('x')
y_sym = sp.Function('y')(x_sym)
ode_sym = sp.Eq(sp.diff(y_sym, x_sym) + 3 * y_sym, x_sym * sp.exp(2 * x_sym))
solution_sym = sp.dsolve(ode_sym, y_sym, ics={y_sym.subs(x_sym, 2): -4})
# Convert symbolic solution to a numerical function
y_analytical_func = sp.lambdify(x_sym, solution_sym.rhs, 'numpy')
y_analytical = y_analytical_func(xsol)
# ---------------------------
# 5. Plot Solutions
# ---------------------------
plt.figure(figsize=(8, 5))
plt.plot(xsol, ysol, 'b-', linewidth=1.5, label='Numerical Solution (RK45)')
plt.plot(xsol, y_analytical, 'r--', linewidth=1.5, label='Analytical Solution')
plt.grid(True, linestyle='--', alpha=0.7)
plt.xlabel('x')
plt.ylabel('y(x)')
plt.title(r'Solution of $y^\prime + 3y = x e^{2x}, \ y(2)=-4$')
plt.legend()
# plt.tight_layout()
# plt.close() # Prevents R Markdown display freeze
# ---------------------------
# 6. Safe Verification (FIXED: NO INDEX ERROR)
# ---------------------------
print("=== Exact Analytical Solution ===")## === Exact Analytical Solution ===## Eq(y(x), x*exp(2*x)/5 - exp(2*x)/25 + (-9*exp(10)/25 - 4*exp(6))*exp(-3*x))##
## === Verification at Key Points ===test_points = [2.0, 2.2, 2.5, 2.7]
# Use precomputed solution for safety (NO SOLVE INSIDE LOOP)
for x_test in test_points:
# Get analytical value
y_an = y_analytical_func(x_test)
# Interpolate numerical solution (safe, no index risk)
y_num = np.interp(x_test, xsol, ysol)
print(f"x = {x_test:.2f}: Numerical y = {y_num:.4f}, Analytical y = {y_an:.4f}")## x = 2.00: Numerical y = -4.0000, Analytical y = -4.0000
## x = 2.20: Numerical y = 19.5986, Analytical y = 19.5980
## x = 2.50: Numerical y = 62.9424, Analytical y = 62.9918
## x = 2.70: Numerical y = 107.8438, Analytical y = 107.8065
3 Example 3 System of First-Order ODEs
Employ Python to solve and plot the solution of the initial value problem: \[ \begin{cases} \dfrac{dx}{dt} = y^2 - x \\[6pt] \dfrac{dy}{dt} = 0.5x^2 - y \end{cases} \quad \begin{bmatrix} x(0) \\ y(0) \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}, \quad t \in [0, 3]. \]
# ---------------------------
# Python Implementation for Example 2.7
# System of ODEs:
# dx/dt = y^2 - x
# dy/dt = 0.5*x^2 - y
# Initial conditions: x(0)=1, y(0)=0, t ∈ [0, 3]
# ---------------------------
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# ---------------------------
# 1. Define the ODE System
# ---------------------------
def ode_system(t, state):
"""
Defines the system of first-order ODEs.
Input:
t: independent variable (time)
state: list [x, y] containing the current values of the state variables
Output:
dstate_dt: list [dx/dt, dy/dt] containing the derivatives
"""
x, y = state
dxdt = y**2 - x
dydt = 0.5 * x**2 - y
return [dxdt, dydt]
# ---------------------------
# 2. Set Problem Parameters
# ---------------------------
t0 = 0 # Initial time
u0 = [1, 0] # Initial state [x(0), y(0)]
tspan = [t0, 3] # Time interval [0, 3]
# ---------------------------
# 3. Solve the System Numerically (RK45, equivalent to MATLAB ode45)
# ---------------------------
# Define fixed evaluation points for stable plotting
t_eval = np.linspace(0, 3, 100)
sol = solve_ivp(
fun=ode_system,
t_span=tspan,
y0=u0,
method='RK45', # Adaptive Runge-Kutta method (matches MATLAB ode45)
t_eval=t_eval # Use fixed points for guaranteed output
)
# Extract results
tsol = sol.t # Time points
xsol = sol.y[0] # Solution x(t)
ysol = sol.y[1] # Solution y(t)
# ---------------------------
# 4. Plot the Solutions
# ---------------------------
plt.figure(figsize=(8, 5))
plt.plot(tsol, xsol, 'b-', linewidth=1.5, label='x(t)')
plt.plot(tsol, ysol, 'r-', linewidth=1.5, label='y(t)')
plt.grid(True, linestyle='--', alpha=0.7)
plt.xlabel('t')
plt.ylabel('x(t), y(t)')
plt.title(r'Solutions of the System: $\frac{dx}{dt}=y^2-x$, $\frac{dy}{dt}=0.5x^2-y$')
plt.legend()
plt.tight_layout()
# Save figure for LaTeX/R Markdown
# plt.savefig('ode_system_plot.png', dpi=300, bbox_inches='tight')
# plt.close() # Prevents display freeze in R Markdown
# ---------------------------
# 5. Print Verification Information
# ---------------------------
print("=== System ODE Solution (RK45) ===")## === System ODE Solution (RK45) ===## Initial conditions: x(0)=1, y(0)=0## Time interval: [0, 3]##
## === Values at Key Time Points ===test_times = [0, 1, 2, 3]
for t in test_times:
# Interpolate to get values at the test time
x_val = np.interp(t, tsol, xsol)
y_val = np.interp(t, tsol, ysol)
print(f"t = {t:.1f}: x(t) = {x_val:.4f}, y(t) = {y_val:.4f}")## t = 0.0: x(t) = 1.0000, y(t) = 0.0000
## t = 1.0: x(t) = 0.3756, y(t) = 0.1175
## t = 2.0: x(t) = 0.1428, y(t) = 0.0601
## t = 3.0: x(t) = 0.0535, y(t) = 0.0245
4 Example 4 System of First-Order ODEs with Fixed Time Step
Employ Python to solve and plot the solution of the initial value problem: \[ \begin{cases} \dfrac{dy_2}{dt} + 3y_2 = 6y_1 \\[6pt] \dfrac{dy_1}{dt} + 5y_1 = \sin(2t) \end{cases} \quad \begin{bmatrix} y_1(0) \\ y_2(0) \end{bmatrix} = \begin{bmatrix} 5 \\ -2 \end{bmatrix}, \quad t \in [0, 0.5]. \] Rearranged into standard form: \[ \begin{cases} \dfrac{dy_1}{dt} = -5y_1 + \sin(2t) \\[6pt] \dfrac{dy_2}{dt} = 6y_1 - 3y_2 \end{cases} \]
# ---------------------------
# Python Implementation for Example 2.8
# System of ODEs:
# dy1/dt = -5*y1 + sin(2*t)
# dy2/dt = 6*y1 - 3*y2
# Initial conditions: y1(0)=5, y2(0)=-2, t ∈ [0, 0.5]
# Fixed time step: 0.001
# ---------------------------
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# ---------------------------
# 1. Define the ODE System
# ---------------------------
def ode_system(t, state):
"""
Defines the system of first-order ODEs in standard form.
Input:
t: independent variable (time)
state: list [y1, y2] containing the current values of the state variables
Output:
dstate_dt: list [dy1/dt, dy2/dt] containing the derivatives
"""
y1, y2 = state
dy1dt = -5 * y1 + np.sin(2 * t)
dy2dt = 6 * y1 - 3 * y2
return [dy1dt, dy2dt]
# ---------------------------
# 2. Set Problem Parameters
# ---------------------------
t0 = 0 # Initial time
y0 = [5, -2] # Initial state [y1(0), y2(0)]
tspan = [t0, 0.5] # Time interval [0, 0.5]
dt = 0.001 # Fixed time step (as specified in the problem)
# ---------------------------
# 3. Solve the System Numerically (RK45, equivalent to MATLAB ode45)
# ---------------------------
# Create time points with fixed step size 0.001 (matches MATLAB's tspan)
t_eval = np.arange(t0, 0.5 + dt, dt)
sol = solve_ivp(
fun=ode_system,
t_span=tspan,
y0=y0,
method='RK45', # Adaptive Runge-Kutta method (matches MATLAB ode45)
t_eval=t_eval # Use fixed time step points for guaranteed output
)
# Extract results
tsol = sol.t # Time points (with step size 0.001)
y1sol = sol.y[0] # Solution y1(t)
y2sol = sol.y[1] # Solution y2(t)
# ---------------------------
# 4. Plot the Solutions
# ---------------------------
plt.figure(figsize=(8, 5))
plt.plot(tsol, y1sol, 'b-', linewidth=1.5, label='y1(t)')
plt.plot(tsol, y2sol, 'r-', linewidth=1.5, label='y2(t)')
plt.grid(True, linestyle='--', alpha=0.7)
plt.xlabel('t')
plt.ylabel('y1(t), y2(t)')
plt.title(r'Solutions of the System: $\frac{dy_1}{dt}=-5y_1+\sin(2t)$, $\frac{dy_2}{dt}=6y_1-3y_2$')
plt.legend()
plt.tight_layout()
# Save figure for LaTeX/R Markdown
# plt.savefig('ode_system_plot_2.png', dpi=300, bbox_inches='tight')
# plt.close() # Prevents display freeze in R Markdown
# ---------------------------
# 5. Print Verification Information
# ---------------------------
print("=== System ODE Solution (RK45) ===")## === System ODE Solution (RK45) ===## Initial conditions: y1(0)=5, y2(0)=-2## Time interval: [0, 0.5]## Fixed time step: 0.001##
## === Values at Key Time Points ===test_times = [0, 0.1, 0.2, 0.3, 0.4, 0.5]
for t in test_times:
# Interpolate to get values at the test time (safe, no index risk)
y1_val = np.interp(t, tsol, y1sol)
y2_val = np.interp(t, tsol, y2sol)
print(f"t = {t:.2f}: y1(t) = {y1_val:.4f}, y2(t) = {y2_val:.4f}")## t = 0.00: y1(t) = 5.0000, y2(t) = -2.0000
## t = 0.10: y1(t) = 3.0408, y2(t) = 0.5352
## t = 0.20: y1(t) = 1.8685, y2(t) = 1.6270
## t = 0.30: y1(t) = 1.1716, y2(t) = 1.9683
## t = 0.40: y1(t) = 0.7616, y2(t) = 1.9447
## t = 0.50: y1(t) = 0.5240, y2(t) = 1.7648
5 Example 5 3-Variable System of ODEs with Fixed Time Step
Employ Python to solve and plot the solution of the initial value problem: \[ \begin{cases} \dfrac{dx}{dt} = y - z \\[6pt] \dfrac{dy}{dt} = 0.45z - e^{-t} \\[6pt] \dfrac{dz}{dt} = -0.25xy + t^2 \end{cases} \quad \begin{bmatrix} x(0) \\ y(0) \\ z(0) \end{bmatrix} = \begin{bmatrix} -2 \\ -5 \\ 7 \end{bmatrix}, \quad t \in [0, 4]. \]
# ---------------------------
# Python Implementation for Example 2.9
# 3-Variable System of ODEs:
# dx/dt = y - z
# dy/dt = 0.45*z - exp(-t)
# dz/dt = -0.25*x*y + t^2
# Initial conditions: x(0)=-2, y(0)=-5, z(0)=7
# Time interval: t ∈ [0, 4]
# Fixed time step: 0.05
# ---------------------------
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# ---------------------------
# 1. Define the ODE System
# ---------------------------
def ode_system(t, state):
"""
Defines the system of three first-order ODEs.
Input:
t: independent variable (time)
state: list [x, y, z] containing the current values of the state variables
Output:
dstate_dt: list [dx/dt, dy/dt, dz/dt] containing the derivatives
"""
x, y, z = state
dxdt = y - z
dydt = 0.45 * z - np.exp(-t)
dzdt = -0.25 * x * y + t**2
return [dxdt, dydt, dzdt]
# ---------------------------
# 2. Set Problem Parameters
# ---------------------------
t0 = 0 # Initial time
u0 = [-2, -5, 7] # Initial state [x(0), y(0), z(0)]
tspan = [t0, 4] # Time interval [0, 4]
dt = 0.05 # Fixed time step (as specified in the problem)
# ---------------------------
# 3. Solve the System Numerically (RK45, equivalent to MATLAB ode45)
# ---------------------------
# Create time points with fixed step size 0.05 (matches MATLAB's tspan)
t_eval = np.arange(t0, 4 + dt, dt)
sol = solve_ivp(
fun=ode_system,
t_span=tspan,
y0=u0,
method='RK45', # Adaptive Runge-Kutta method (matches MATLAB ode45)
t_eval=t_eval # Use fixed time step points for guaranteed output
)
# Extract results
tsol = sol.t # Time points (with step size 0.05)
xsol = sol.y[0] # Solution x(t)
ysol = sol.y[1] # Solution y(t)
zsol = sol.y[2] # Solution z(t)
# ---------------------------
# 4. Plot the Solutions
# ---------------------------
plt.figure(figsize=(8, 5))
plt.plot(tsol, xsol, 'b-', linewidth=1.5, label='x(t)')
plt.plot(tsol, ysol, 'r-', linewidth=1.5, label='y(t)')
plt.plot(tsol, zsol, 'y-', linewidth=1.5, label='z(t)')
plt.grid(True, linestyle='--', alpha=0.7)
plt.xlabel('t')
plt.ylabel('x(t), y(t), z(t)')
plt.title(r'Solutions of the 3-Variable ODE System')
plt.legend()
plt.tight_layout()
# Save figure for LaTeX/R Markdown
# plt.savefig('ode_system_plot_3.png', dpi=300, bbox_inches='tight')
# plt.close() # Prevents display freeze in R Markdown
# ---------------------------
# 5. Print Verification Information
# ---------------------------
print("=== 3-Variable System ODE Solution (RK45) ===")## === 3-Variable System ODE Solution (RK45) ===## Initial conditions: x(0)=-2, y(0)=-5, z(0)=7## Time interval: [0, 4]## Fixed time step: 0.05##
## === Values at Key Time Points ===test_times = [0, 1, 2, 3, 4]
for t in test_times:
# Interpolate to get values at the test time (safe, no index risk)
x_val = np.interp(t, tsol, xsol)
y_val = np.interp(t, tsol, ysol)
z_val = np.interp(t, tsol, zsol)
print(f"t = {t:.1f}: x(t) = {x_val:.4f}, y(t) = {y_val:.4f}, z(t) = {z_val:.4f}")## t = 0.0: x(t) = -2.0000, y(t) = -5.0000, z(t) = 7.0000
## t = 1.0: x(t) = -10.3561, y(t) = -3.7487, z(t) = 0.3960
## t = 2.0: x(t) = -10.4528, y(t) = -5.9611, z(t) = -9.7007
## t = 3.0: x(t) = -4.2575, y(t) = -12.9089, z(t) = -19.7099
## t = 4.0: x(t) = -2.4123, y(t) = -21.5900, z(t) = -17.1283
6 Example 6 Stiff 3-Variable ODE System (Robertson Reaction Kinetics)
Employ Python to solve and plot the solution of the initial value problem: \[ \begin{cases} \dfrac{dy_1}{dt} = -k_1 y_1 + k_3 y_2 y_3 \\[6pt] \dfrac{dy_2}{dt} = k_1 y_1 - k_2 y_2^2 - k_3 y_2 y_3 \\[6pt] \dfrac{dy_3}{dt} = k_2 y_2^2 \end{cases} \quad \begin{bmatrix} y_1(0) \\ y_2(0) \\ y_3(0) \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}, \] with rate constants \(k_1 = 0.04\), \(k_2 = 10^4\), \(k_3 = 3 \times 10^7\), over the interval \(t \in [0, 4 \times 10^6]\).
# ---------------------------
# Python Implementation for Example 2.13
# Stiff 3-Variable ODE System (Robertson Reaction Kinetics)
# dy1/dt = -k1*y1 + k3*y2*y3
# dy2/dt = k1*y1 - k2*y2^2 - k3*y2*y3
# dy3/dt = k2*y2^2
# Initial conditions: y1(0)=1, y2(0)=0, y3(0)=0
# Rate constants: k1=0.04, k2=1e4, k3=3e7
# Time interval: t ∈ [0, 4e6]
# ---------------------------
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# ---------------------------
# 1. Define Parameters and ODE System
# ---------------------------
k1 = 0.04
k2 = 1e4
k3 = 3e7
def ode_system(t, state):
"""
Defines the stiff system of three first-order ODEs (Robertson kinetics).
Input:
t: independent variable (time)
state: list [y1, y2, y3] containing the current values of the state variables
Output:
dstate_dt: list [dy1/dt, dy2/dt, dy3/dt] containing the derivatives
"""
y1, y2, y3 = state
dy1dt = -k1 * y1 + k3 * y2 * y3
dy2dt = k1 * y1 - k2 * y2**2 - k3 * y2 * y3
dy3dt = k2 * y2**2
return [dy1dt, dy2dt, dy3dt]
# ---------------------------
# 2. Set Problem Parameters
# ---------------------------
t0 = 0
u0 = [1, 0, 0]
tspan = [t0, 4e6]
# ---------------------------
# 3. Solve ODE (FIXED: No t_eval error)
# ---------------------------
# SOLUTION: Use linear space instead of logspace (avoids 0 → 1 jump)
t_eval = np.linspace(0, 4e6, 200)
sol = solve_ivp(
fun=ode_system,
t_span=tspan,
y0=u0,
method='Radau', # Stiff solver (like MATLAB ode15s)
t_eval=t_eval
)
# Extract results safely
tsol = sol.t
y1sol = sol.y[0]
y2sol = sol.y[1]
y3sol = sol.y[2]
# ---------------------------
# 4. Plot Solutions
# ---------------------------
plt.figure(figsize=(8, 10))
plt.subplot(311)
plt.plot(tsol, y1sol, 'b-', linewidth=1.5)
plt.grid(True, alpha=0.7)
plt.ylabel('$y_1(t)$')
plt.title('Stiff Robertson ODE System Solution')
plt.subplot(312)
plt.plot(tsol, y2sol, 'r-', linewidth=1.5)
plt.grid(True, alpha=0.7)
plt.ylabel('$y_2(t)$')
plt.subplot(313)
plt.plot(tsol, y3sol, 'g-', linewidth=1.5)
plt.grid(True, alpha=0.7)
plt.xlabel('t')
plt.ylabel('$y_3(t)$')
# plt.tight_layout()
# plt.close() # Critical for R Markdown
# ---------------------------
# 5. Print Verification (SAFE: No index / out-of-bounds error)
# ---------------------------
print("=== Stiff 3-Variable ODE System Solution ===")## === Stiff 3-Variable ODE System Solution ===## Initial conditions: y1(0)=1, y2(0)=0, y3(0)=0## Time interval: [0, 4e6]##
## === Key Time Points ===test_times = [0, 1e6, 2e6, 3e6, 4e6]
for t in test_times:
y1_val = np.interp(t, tsol, y1sol)
y2_val = np.interp(t, tsol, y2sol)
y3_val = np.interp(t, tsol, y3sol)
print(f"t = {t:.1e} | y1 = {y1_val:.6f} | y2 = {y2_val:.6e} | y3 = {y3_val:.6f}")## t = 0.0e+00 | y1 = 1.000000 | y2 = 0.000000e+00 | y3 = 0.000000
## t = 1.0e+06 | y1 = 0.996242 | y2 = 3.535549e-07 | y3 = 0.003757
## t = 2.0e+06 | y1 = 0.995268 | y2 = 2.804462e-07 | y3 = 0.004731
## t = 3.0e+06 | y1 = 0.994586 | y2 = 2.449147e-07 | y3 = 0.005414
## t = 4.0e+06 | y1 = 0.994042 | y2 = 2.224892e-07 | y3 = 0.005958