SciPy numerical toolkit for physics: ODE/PDE solving, FFT analysis, optimization, numerical integration, and sparse linear algebra with real-world examples.
SciPy numerical toolkit for physics: ODE/PDE solving, FFT analysis, optimization, numerical integration, and sparse linear algebra with real-world examples.
SciPy provides a comprehensive suite of numerical algorithms essential for physics
simulation and data analysis. This skill covers ODE solving, PDE discretization,
numerical integration, FFT-based spectral analysis, nonlinear optimization, and
sparse linear algebra.
1. ODE Solving with solve_ivp
SciPy's solve_ivp supports multiple integration methods. Choose the right one:
Method
Best for
Notes
RK45
Non-stiff, smooth solutions
Default, 4th/5th order
RK23
Non-stiff, low accuracy
Cheaper per step
DOP853
Non-stiff, high accuracy
8th order Dormand-Prince
Radau
Stiff problems
Implicit, expensive but robust
BDF
Very stiff (e.g., chemistry)
Backward differentiation
LSODA
Automatic stiff detection
Wraps ODEPACK
1.1 Lorenz Attractor
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
deflorenz_system(t, state, sigma=10.0, rho=28.0, beta=8.0 / 3.0):
"""
Lorenz attractor ODEs.
dx/dt = sigma * (y - x)
dy/dt = x * (rho - z) - y
dz/dt = x * y - beta * z
"""
x, y, z = state
dxdt = sigma * (y - x)
dydt = x * (rho - z) - y
dzdt = x * y - beta * z
return [dxdt, dydt, dzdt]
defsimulate_lorenz(t_span=(0, 50), t_eval=, initial_state=, method=):
initial_state :
initial_state = [, , ]
t_eval :
t_eval = np.linspace(t_span[], t_span[], )
sol = solve_ivp(
lorenz_system,
t_span,
initial_state,
method=method,
t_eval=t_eval,
rtol=,
atol=,
dense_output=,
)
sol
():
fig = plt.figure(figsize=(, ))
ax = fig.add_subplot(, projection=)
ax.plot(sol.y[], sol.y[], sol.y[], lw=, alpha=, color=)
ax.set_xlabel()
ax.set_ylabel()
ax.set_zlabel()
ax.set_title()
plt.tight_layout()
fig
__name__ == :
sol = simulate_lorenz()
()
()
fig = plot_lorenz_attractor(sol)
plt.show()
None
None
"RK45"
"""Simulate the Lorenz attractor and return the solution object."""
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
import matplotlib.pyplot as plt
defparticle_in_box_hamiltonian(N=500, L=1.0, hbar=1.0, m=1.0):
"""
Construct the 1D particle-in-a-box Hamiltonian using finite differences.
H = -hbar^2/(2m) * d^2/dx^2
Discretized on N interior points with Dirichlet boundary conditions.
Returns
-------
H : sparse matrix (CSC)
Hamiltonian matrix.
x : ndarray
Grid points.
"""
dx = L / (N + 1)
x = np.linspace(dx, L - dx, N)
diag_main = 2.0 * np.ones(N)
diag_off = -1.0 * np.ones(N - 1)
kinetic = sp.diags([diag_off, diag_main, diag_off], [-1, 0, 1], format="csc")
H = (hbar ** 2 / (2 * m * dx ** 2)) * kinetic
return H, x
defsolve_particle_in_box(N=500, L=1.0, n_eigvals=6):
"""
Solve for the lowest n_eigvals eigenstates of the particle in a box.
Analytical energies: E_n = n^2 * pi^2 * hbar^2 / (2 m L^2)
"""
H, x = particle_in_box_hamiltonian(N=N, L=L)
# Use eigsh for symmetric matrices — much faster than full diagonalization
eigenvalues, eigenvectors = spla.eigsh(H, k=n_eigvals, which="SM")
# Sort by energy
idx = np.argsort(eigenvalues)
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
# Normalize eigenvectors
dx = L / (N + 1)
for i inrange(n_eigvals):
norm = np.sqrt(np.trapz(eigenvectors[:, i] ** 2, dx=dx))
eigenvectors[:, i] /= norm
# Analytical energies (hbar=m=1)
n_vals = np.arange(1, n_eigvals + 1)
E_exact = n_vals ** 2 * np.pi ** 2 / (2 * L ** 2)
print("\nParticle-in-a-box energy levels (hbar=m=1):")
print(f"{'n':>4}{'Numerical':>14}{'Analytical':>14}{'Rel. Error':>12}")
for i, (En, Ea) inenumerate(zip(eigenvalues, E_exact)):
print(f"{i+1:>4}{En:>14.6f}{Ea:>14.6f}{abs(En-Ea)/Ea:>12.2e}")
# Plot wavefunctions
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
for i inrange(n_eigvals):
psi = eigenvectors[:, i]
# Fix sign conventionif psi[N // 4] < 0:
psi = -psi
axes[0].plot(x, psi + eigenvalues[i], label=f"n={i+1}", lw=1.5)
axes[0].set_xlabel("x / L")
axes[0].set_ylabel("Energy + ψ(x)")
axes[0].set_title("Particle-in-a-Box Wavefunctions")
axes[0].legend(loc="upper left", fontsize=8)
axes[1].scatter(n_vals, eigenvalues, label="Numerical", zorder=5)
axes[1].plot(n_vals, E_exact, "r--", label="Analytical", lw=1.5)
axes[1].set_xlabel("Quantum number n")
axes[1].set_ylabel("Energy")
axes[1].set_title("Energy Levels Comparison")
axes[1].legend()
plt.tight_layout()
plt.show()
return eigenvalues, eigenvectors, x
if __name__ == "__main__":
solve_particle_in_box(N=800, n_eigvals=8)
7. Complete Example A — Damped Pendulum: RK45 vs Radau
This comprehensive example wraps everything together: ODE solving with two stiff
and non-stiff methods, phase-space analysis, and Poincaré section extraction.
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
defpendulum_ode(t, y, b=0.25, c=5.0):
"""
Simple pendulum ODE:
theta'' + b * theta' + c * sin(theta) = 0
State: y = [theta, omega]
"""
theta, omega = y
dydt = [omega, -b * omega - c * np.sin(theta)]
return dydt
defpendulum_event(t, y, **kwargs):
"""Event: detect zero-crossings of theta (Poincaré section)."""return y[0]
pendulum_event.terminal = False
pendulum_event.direction = 1defrun_pendulum_comparison(theta0=np.pi - 0.1, omega0=0.0, t_end=30.0):
"""
Solve the pendulum ODE with RK45 and Radau, compare trajectories.
"""
y0 = [theta0, omega0]
t_eval = np.linspace(0, t_end, 3000)
solutions = {}
for method in ["RK45", "Radau"]:
sol = solve_ivp(
pendulum_ode,
(0, t_end),
y0,
method=method,
t_eval=t_eval,
rtol=1e-9,
atol=1e-11,
events=pendulum_event,
)
solutions[method] = sol
print(f"{method}: nfev={sol.nfev}, success={sol.success}, steps={sol.t.size}")
# Compute energy (should be conserved for b=0)defenergy(theta, omega, c=5.0):
return0.5 * omega ** 2 - c * np.cos(theta)
fig, axes = plt.subplots(2, 2, figsize=(13, 8))
colors = {"RK45": "steelblue", "Radau": "darkorange"}
for method, sol in solutions.items():
theta = sol.y[0]
omega = sol.y[1]
E = energy(theta, omega)
axes[0, 0].plot(sol.t, theta, label=method, color=colors[method], lw=1)
axes[0, 1].plot(theta, omega, color=colors[method], label=method, lw=0.8, alpha=0.8)
axes[1, 0].plot(sol.t, E - E[0], color=colors[method], label=method, lw=1)
# Difference between methods
theta_rk45 = solutions["RK45"].y[0]
theta_radau = solutions["Radau"].y[0]
diff = np.abs(theta_rk45 - theta_radau)
axes[1, 1].semilogy(solutions["RK45"].t, diff + 1e-16, color="purple", lw=1)
axes[1, 1].set_title("RK45 vs Radau |Δθ|")
axes[1, 1].set_xlabel("Time (s)")
axes[1, 1].set_ylabel("|Δθ| (rad)")
axes[0, 0].set_title("Pendulum Angle")
axes[0, 0].set_xlabel("Time (s)")
axes[0, 0].set_ylabel("θ (rad)")
axes[0, 0].legend()
axes[0, 1].set_title("Phase Portrait")
axes[0, 1].set_xlabel("θ (rad)")
axes[0, 1].set_ylabel("ω (rad/s)")
axes[0, 1].legend()
axes[1, 0].set_title("Energy Drift (b=0.25, damped)")
axes[1, 0].set_xlabel("Time (s)")
axes[1, 0].set_ylabel("ΔE")
axes[1, 0].legend()
plt.suptitle("Damped Pendulum: RK45 vs Radau", fontsize=13, fontweight="bold")
plt.tight_layout()
plt.show()
if __name__ == "__main__":
run_pendulum_comparison()
8. Complete Example B — Frequency Analysis Pipeline
