Skip to main content Skills Marketplace Découvrez et explorez les compétences IA créées par la communauté.
Installer avec Codex ou Claude Copiez ce prompt, collez-le dans Codex, Claude ou un autre assistant, puis laissez-le vérifier la page du skill et l'installer pour vous.
Afficher les détails du prompt Une commande directe contourne le prompt de vérification. Examinez la source avant de l'exécuter.
npx skills add https://github.com/xjtulyc/awesome-rosetta-skills --skill sympy-symbolicLa commande reste sur une seule ligne. Faites défiler horizontalement pour la vérifier avant de la copier.
Vous préférez une copie locale ? Téléchargez les fichiers actuellement disponibles dans SkillsMP.
Télécharger Zip Téléchargement... Métiers associés SOC
Basé sur la classification professionnelle SOC
name sympy-symbolic description SymPy symbolic computation for physics: algebra, calculus, linear algebra, ODE solving, Lagrangian mechanics, quantum physics, and code generation.
tags ["sympy","symbolic-computation","physics","mechanics","quantum","code-generation"] version 1.0.0 authors [{"name":"awesome-rosetta-skills contributors","github":"@awesome-rosetta-skills"}] license MIT platforms ["claude-code","codex","gemini-cli","cursor"] dependencies ["sympy>=1.12","numpy>=1.24.0","matplotlib>=3.7.0","scipy>=1.11.0"] last_updated 2026-03-17
SymPy Symbolic Computation for Physics
SymPy is a pure-Python computer algebra system. This skill demonstrates how to use
SymPy for symbolic physics: manipulating equations, solving ODEs analytically,
deriving equations of motion via Lagrangian mechanics, performing quantum mechanics
calculations, and converting symbolic expressions to fast numerical code.
1. Symbolic Algebra Fundamentals
import sympy as sp
x, y, z = sp.symbols("x y z" , real=True )
t = sp.Symbol("t" , positive=True )
n = sp.Symbol("n" , integer=True , positive=True )
epsilon, omega, phi = sp.symbols("epsilon omega phi" , positive=True )
expr = (x + y) ** 4
print ("Expanded:" , sp.expand(expr))
factored = sp.factor(x**4 - 1 )
print ("Factored:" , factored)
trig_expr = sp.sin(x) ** 2 + sp.cos(x) ** 2
print ("Trig simplified:" , sp.trigsimp(trig_expr))
f = sp.exp(sp.I * x)
series = sp.series(f, x, 0 , 6 )
print ("Euler series:" , series)
quadratic = sp.Eq(x**2 + 3 *x + 2 , 0 )
solutions = sp.solve(quadratic, x)
print ("Quadratic roots:" , solutions)
eq1 = sp.Eq(2 *x + y, 5 )
eq2 = sp.Eq(x - y, 1
print
"System solution:"
2. Symbolic Calculus import sympy as sp
x, t, a, b = sp.symbols("x t a b" , real=True )
f = sp.sin(x**2 ) * sp.exp(-x)
df = sp.diff(f, x)
d2f = sp.diff(f, x, 2 )
print ("f(x) =" , f)
print ("f'(x) =" , sp.simplify(df))
print ("f''(x) =" , sp.simplify(d2f))
phi = sp.Function("phi" )(x, t)
wave_eq = sp.diff(phi, t, 2 ) - sp.diff(phi, x, 2 )
print ("\nWave equation LHS:" , wave_eq)
I1 = sp.integrate(sp.exp(-x**2 ), (x, -sp.oo, sp.oo))
print ("\nGaussian integral:" , I1)
I2 = sp.integrate(sp.sin(x) * sp.cos(x), x)
print ("Trig integral:" , sp.simplify(I2))
I3 = sp.integrate(sp.exp(-a * x**2 ), (x, 0 , sp.oo))
print ("Parameterized Gaussian:" , sp.simplify(I3))
lim1 = sp.limit(sp.sin(x) / x, x, 0 )
print ("\nlim_{x->0} sin(x)/x =" , lim1)
lim2 = sp.limit((1 + 1 /x)**x, x, sp.oo)
print ("lim_{x->inf} (1+1/x)^x =" , lim2)
def fourier_coefficient (f, L, n, var ):
"""Compute the nth Fourier coefficient of f on [-L, L]."""
a_n = (1 / L) * sp.integrate(f * sp.cos(n * sp.pi * var / L), (var, -L, L))
b_n = (1 / L) * sp.integrate(f * sp.sin(n * sp.pi * var / L), (var, -L, L))
return sp.simplify(a_n), sp.simplify(b_n)
n = sp.Symbol("n" , integer=True , positive=True )
L = sp.pi
a_n, b_n = fourier_coefficient(x, L, n, x)
print ("\nFourier coefficients of f(x)=x on [-pi,pi]:" )
print (" a_n =" , a_n)
print (" b_n =" , b_n)
3. Symbolic Linear Algebra import sympy as sp
a, b, c, d = sp.symbols("a b c d" )
M = sp.Matrix([[a, b], [c, d]])
print ("Matrix M:\n" , M)
print ("Determinant:" , M.det())
print ("Inverse:\n" , sp.simplify(M.inv()))
eigenvals = M.eigenvals()
eigenvects = M.eigenvects()
print ("\nEigenvalues:" , eigenvals)
for eigenval, mult, vecs in eigenvects:
print (f" λ = {eigenval} , multiplicity = {mult} , eigenvector = {vecs[0 ].T} " )
sx = sp.Matrix([[0 , 1 ], [1 , 0 ]])
sy = sp.Matrix([[0 , -sp.I], [sp.I, 0 ]])
sz = sp.Matrix([[1 , 0 ], [0 , -1 ]])
print ("\nPauli matrix commutation relations:" )
print ("[sx, sy] = 2i*sz:" , sp.simplify(sx * sy - sy * sx) == 2 * sp.I * sz)
print ("[sy, sz] = 2i*sx:" , sp.simplify(sy * sz - sz * sy) == 2 * sp.I * sx)
print ("[sz, sx] = 2i*sy:" , sp.simplify(sz * sx - sx * sz) == 2 * sp.I * sy)
theta = sp.Symbol("theta" , real=True )
R = sp.Matrix([
[sp.cos(theta), -sp.sin(theta)],
[sp.sin(theta), sp.cos(theta)],
])
print ("\nRotation matrix R @ R.T:" )
print (sp.simplify(R * R.T))
print ("det(R):" , sp.simplify(R.det()))
A = sp.Matrix([[2 , 1 , -1 ], [-3 , -1 , 2 ], [-2 , 1 , 2 ]])
b_vec = sp.Matrix([8 , -11 , -3 ])
solution = A.solve(b_vec)
print ("\nLinear system solution:" , solution.T)
4. Solving Differential Equations Symbolically import sympy as sp
t = sp.Symbol("t" , positive=True )
x = sp.Function("x" )
omega = sp.Symbol("omega" , positive=True )
ode_sho = sp.Eq(x(t).diff(t, 2 ) + omega**2 * x(t), 0 )
sol_sho = sp.dsolve(ode_sho, x(t))
print ("SHO solution:" , sol_sho)
C1, C2 = sp.symbols("C1 C2" )
ics = {x(0 ): 1 , x(t).diff(t).subs(t, 0 ): 0 }
sol_ics = sp.dsolve(ode_sho, x(t), ics=ics)
print ("SHO with ICs:" , sol_ics)
gamma, omega0 = sp.symbols("gamma omega_0" , positive=True )
ode_damped = sp.Eq(x(t).diff(t, 2 ) + 2 *gamma*x(t).diff(t) + omega0**2 *x(t), 0 )
sol_damped = sp.dsolve(ode_damped, x(t))
print ("\nDamped oscillator solution:" , sp.simplify(sol_damped.rhs))
F0, omega_d = sp.symbols("F_0 omega_d" , positive=True )
ode_driven = sp.Eq(
x(t).diff(t, 2 ) + 2 *gamma*x(t).diff(t) + omega0**2 *x(t),
F0 * sp.cos(omega_d * t)
)
sol_driven = sp.dsolve(ode_driven, x(t))
print ("\nDriven oscillator particular solution:" )
print (sp.simplify(sol_driven.rhs))
x_var = sp.Symbol("x" , real=True )
alpha, k = sp.symbols("alpha k" , positive=True )
u = sp.exp(-alpha * k**2 * t) * sp.sin(k * x_var)
residual = sp.diff(u, t) - alpha * sp.diff(u, x_var, 2 )
print ("\nHeat equation residual (should be 0):" , sp.simplify(residual))
5. Lagrangian Mechanics — Double Pendulum Derive the equations of motion for a double pendulum using the Lagrangian formalism.
import sympy as sp
from sympy.physics.mechanics import dynamicsymbols, LagrangesMethod, Particle, Point, ReferenceFrame
def derive_double_pendulum ():
"""
Derive equations of motion for a planar double pendulum via Lagrangian mechanics.
q1 = angle of rod 1 from vertical
q2 = angle of rod 2 from vertical
"""
t = sp.Symbol("t" )
m1, m2, l1, l2, g = sp.symbols("m_1 m_2 l_1 l_2 g" , positive=True )
q1 = dynamicsymbols("q_1" )
q2 = dynamicsymbols("q_2" )
q1d = dynamicsymbols("q_1" , 1 )
q2d = dynamicsymbols("q_2" , 1 )
x1 = l1 * sp.sin(q1)
y1 = -l1 * sp.cos(q1)
x2 = x1 + l2 * sp.sin(q2)
y2 = y1 - l2 * sp.cos(q2)
x1d = sp.diff(x1, t)
y1d = sp.diff(y1, t)
x2d = sp.diff(x2, t)
y2d = sp.diff(y2, t)
T = sp.Rational(1 , 2 ) * m1 * (x1d**2 + y1d**2 ) + \
sp.Rational(1 , 2 ) * m2 * (x2d**2 + y2d**2 )
V = m1 * g * y1 + m2 * g * y2
L = T - V
L = sp.simplify(L)
print ("Lagrangian L = T - V:" )
print (sp.simplify(L))
def euler_lagrange (L, q, qd, t ):
"""Compute d/dt(dL/dqd) - dL/dq = 0."""
dL_dqd = sp.diff(L, qd)
d_dt_dL_dqd = sp.diff(dL_dqd, t)
dL_dq = sp.diff(L, q)
return sp.simplify(d_dt_dL_dqd - dL_dq)
eq1 = euler_lagrange(L, q1, q1d, t)
eq2 = euler_lagrange(L, q2, q2d, t)
print ("\nEuler-Lagrange equation 1 (q1):" )
print (sp.simplify(eq1))
print ("\nEuler-Lagrange equation 2 (q2):" )
print (sp.simplify(eq2))
def linearize (expr, coords ):
result = expr
for q in coords:
result = result.subs(sp.sin(q), q).subs(sp.cos(q), 1 )
return sp.simplify(result)
eq1_lin = linearize(eq1, [q1, q2])
eq2_lin = linearize(eq2, [q1, q2])
print ("\nLinearized EOM (small angles):" )
print ("EOM1:" , eq1_lin)
print ("EOM2:" , eq2_lin)
return eq1, eq2, eq1_lin, eq2_lin
if __name__ == "__main__" :
eqs = derive_double_pendulum()
6. Quantum Harmonic Oscillator import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
def quantum_harmonic_oscillator ():
"""
Symbolic solution of the quantum harmonic oscillator.
H psi_n = E_n psi_n
E_n = hbar * omega * (n + 1/2)
psi_n(x) = N_n * H_n(xi) * exp(-xi^2/2) where xi = sqrt(m*omega/hbar)*x
"""
x, xi = sp.symbols("x xi" , real=True )
n = sp.Symbol("n" , nonnegative=True , integer=True )
hbar, omega, m = sp.symbols("hbar omega m" , positive=True )
E_n = hbar * omega * (n + sp.Rational(1 , 2 ))
print ("Energy eigenvalues:" )
for k in range (6 ):
print (f" E_{k} = {E_n.subs(n, k)} * hbar*omega" )
def hermite (n_val, xi_sym ):
"""Generate Hermite polynomial H_n(xi) via recurrence."""
if n_val == 0 :
return sp.Integer(1 )
elif n_val == 1 :
return 2 * xi_sym
else :
H_prev2 = sp.Integer(1 )
H_prev1 = 2 * xi_sym
for k in range (2 , n_val + 1 ):
H_curr = 2 * xi_sym * H_prev1 - 2 * (k - 1 ) * H_prev2
H_prev2 = H_prev1
H_prev1 = H_curr
return sp.expand(H_curr)
def norm_constant (n_val ):
return 1 / sp.sqrt(2 **n_val * sp.factorial(n_val) * sp.sqrt(sp.pi))
print ("\nWavefunctions psi_n(xi):" )
wavefunctions = []
for k in range (5 ):
Hn = hermite(k, xi)
Nn = norm_constant(k)
psi = Nn * Hn * sp.exp(-xi**2 / 2 )
psi_simplified = sp.simplify(psi)
wavefunctions.append(psi_simplified)
print (f" psi_{k} (xi) = {psi_simplified} " )
inner_00 = sp.integrate(wavefunctions[0 ]**2 , (xi, -sp.oo, sp.oo))
inner_01 = sp.integrate(wavefunctions[0 ] * wavefunctions[1 ], (xi, -sp.oo, sp.oo))
inner_11 = sp.integrate(wavefunctions[1 ]**2 , (xi, -sp.oo, sp.oo))
print (f"\nOrthonormality check:" )
print (f" <0|0> = {sp.simplify(inner_00)} " )
print (f" <0|1> = {sp.simplify(inner_01)} " )
print (f" <1|1> = {sp.simplify(inner_11)} " )
psi_funcs = [sp.lambdify(xi, psi, "numpy" ) for psi in wavefunctions]
xi_vals = np.linspace(-4 , 4 , 500 )
fig, ax = plt.subplots(figsize=(10 , 6 ))
colors = ["steelblue" , "darkorange" , "green" , "red" , "purple" ]
for k, (psi_np, color) in enumerate (zip (psi_funcs, colors)):
psi_vals = psi_np(xi_vals)
ax.plot(xi_vals, psi_vals + k + 0.5 , label=f"ψ_{k} (ξ)" , color=color, lw=1.8 )
ax.axhline(k + 0.5 , color=color, ls="--" , alpha=0.3 )
V = xi_vals**2 / 2
ax.plot(xi_vals, V / 5 , "k-" , lw=1.5 , alpha=0.3 , label="V(ξ)/5" )
ax.set_xlim(-4 , 4 )
ax.set_ylim(-0.5 , 5.5 )
ax.set_xlabel("ξ = x/x₀" )
ax.set_ylabel("ψₙ(ξ) + Eₙ/ℏω" )
ax.set_title("Quantum Harmonic Oscillator Wavefunctions" )
ax.legend(loc="upper right" )
plt.tight_layout()
plt.show()
return wavefunctions
if __name__ == "__main__" :
quantum_harmonic_oscillator()
7. Code Generation: lambdify and CSE import sympy as sp
import numpy as np
import time
def code_generation_demo ():
"""
Demonstrate converting symbolic expressions to fast numerical functions.
lambdify: converts SymPy expression to numpy/scipy callable.
cse: common subexpression elimination for optimized computation.
"""
x, y, z = sp.symbols("x y z" , real=True )
expr = (sp.sin(x) * sp.cos(y) + sp.exp(-x**2 - y**2 )) / (1 + sp.sqrt(x**2 + y**2 + z**2 ))
f_np = sp.lambdify([x, y, z], expr, modules="numpy" )
x_arr = np.random.randn(100_000 )
y_arr = np.random.randn(100_000 )
z_arr = np.random.randn(100_000 )
t0 = time.perf_counter()
result = f_np(x_arr, y_arr, z_arr)
print (f"lambdify evaluation on 100k points: {time.perf_counter()-t0:.4 f} s" )
print (f" Result sample: {result[:5 ]} " )
expr_list = [
sp.sin(x)**2 + sp.cos(x)*sp.cos(y),
sp.sin(x)**2 * sp.exp(-y),
sp.cos(x)**2 - sp.sin(x)*sp.sin(y),
]
replacements, reduced_exprs = sp.cse(expr_list)
print ("\nCSE replacements:" )
for sym, val in replacements:
print (f" {sym} = {val} " )
print ("Reduced expressions:" )
for e in reduced_exprs:
print (f" {e} " )
from sympy.printing.pycode import pycode
print ("\nGenerated Python code snippet:" )
for sym, val in replacements:
print (f" {sym} = {pycode(val)} " )
omega, t_sym = sp.symbols("omega t" , real=True )
gamma_sym = sp.Symbol("gamma" , positive=True )
transfer_fn = omega**2 / (omega**2 - t_sym**2 + 2 *sp.I*gamma_sym*t_sym)
print ("\nLaTeX representation of transfer function:" )
print (sp.latex(transfer_fn))
from scipy.integrate import quad
integrand_sym = sp.exp(-x**2 ) * sp.cos(x)
integrand_np = sp.lambdify(x, integrand_sym, modules="numpy" )
result_int, err = quad(integrand_np, -np.inf, np.inf)
exact = np.sqrt(np.pi) * np.exp(-0.25 )
print (f"\nIntegral of exp(-x^2)*cos(x): {result_int:.8 f} (exact: {exact:.8 f} )" )
if __name__ == "__main__" :
code_generation_demo()
8. Maxwell Equations Symbolic Manipulation import sympy as sp
from sympy.vector import CoordSys3D, curl, divergence, gradient
def maxwell_equations_demo ():
"""
Demonstrate Maxwell equations using SymPy's vector calculus.
In Gaussian units:
div E = 4*pi*rho
div B = 0
curl E = -1/c * dB/dt
curl B = 4*pi/c * J + 1/c * dE/dt
"""
N = CoordSys3D("N" )
x, y, z, t = sp.symbols("x y z t" , real=True )
epsilon0, mu0, c = sp.symbols("epsilon_0 mu_0 c" , positive=True )
omega_k, k_val = sp.symbols("omega k" , positive=True )
E0 = sp.Symbol("E_0" , real=True )
E_x = E0 * sp.cos(k_val * z - omega_k * t)
E_field = E_x * N.i + 0 * N.j + 0 * N.k
B_y = E0 / c * sp.cos(k_val * z - omega_k * t)
B_field = 0 * N.i + B_y * N.j + 0 * N.k
div_E = divergence(E_field)
print ("div(E) =" , sp.simplify(div_E))
div_B = divergence(B_field)
print ("div(B) =" , sp.simplify(div_B))
curl_E = curl(E_field)
dB_dt = (E0 / c) * omega_k * sp.sin(k_val * z - omega_k * t) * N.j
faraday_lhs = curl_E
faraday_rhs = -1 /c * dB_dt
print ("\nFaraday's law verification:" )
print (" curl(E) =" , faraday_lhs)
print (" -1/c * dB/dt =" , sp.simplify(faraday_rhs))
substituted = sp.simplify(
sp.nsimplify(faraday_lhs.dot(N.j) - faraday_rhs.dot(N.j)).subs(k_val, omega_k / c)
)
print (f" Difference (k=omega/c): {substituted} " )
print ("\nDispersion relation for electromagnetic waves in vacuum:" )
print (" k = omega/c => v_phase = omega/k = c" )
return E_field, B_field
if __name__ == "__main__" :
maxwell_equations_demo()
9. Complete Example A — Double Pendulum: Symbolic → Numerical import sympy as sp
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
def double_pendulum_symbolic_to_numerical ():
"""
(a) Derive double pendulum EOM symbolically
(b) Lambdify to fast numerical functions
(c) Simulate and plot
"""
t = sp.Symbol("t" )
g_sym, l_sym, m_sym = sp.symbols("g l m" , positive=True )
q1 = sp.Function("q1" )(t)
q2 = sp.Function("q2" )(t)
dq1 = q1.diff(t)
dq2 = q2.diff(t)
x1 = l_sym * sp.sin(q1)
y1 = -l_sym * sp.cos(q1)
x2 = x1 + l_sym * sp.sin(q2)
y2 = y1 - l_sym * sp.cos(q2)
T = sp.Rational(1 , 2 ) * m_sym * (
(x1.diff(t))**2 + (y1.diff(t))**2 +
(x2.diff(t))**2 + (y2.diff(t))**2
)
T = sp.trigsimp(sp.expand(T))
V = m_sym * g_sym * (y1 + y2)
L = sp.expand(T - V)
dL_dq1 = L.diff(q1)
dL_ddq1 = L.diff(dq1)
EL1 = sp.expand(dL_ddq1.diff(t) - dL_dq1)
dL_dq2 = L.diff(q2)
dL_ddq2 = L.diff(dq2)
EL2 = sp.expand(dL_ddq2.diff(t) - dL_dq2)
ddq1, ddq2 = sp.symbols("ddq1 ddq2" )
EL1_sub = EL1.subs({q1.diff(t, 2 ): ddq1, q2.diff(t, 2 ): ddq2})
EL2_sub = EL2.subs({q1.diff(t, 2 ): ddq1, q2.diff(t, 2 ): ddq2})
sol = sp.solve([EL1_sub, EL2_sub], [ddq1, ddq2])
ddq1_expr = sol[ddq1]
ddq2_expr = sol[ddq2]
Q1, W1, Q2, W2 = sp.symbols("Q1 W1 Q2 W2" , real=True )
subs_dict = {q1: Q1, dq1: W1, q2: Q2, dq2: W2}
ddq1_num = ddq1_expr.subs(subs_dict)
ddq2_num = ddq2_expr.subs(subs_dict)
g_val, l_val, m_val = 9.81 , 1.0 , 1.0
param_subs = {g_sym: g_val, l_sym: l_val, m_sym: m_val}
ddq1_np = sp.lambdify([Q1, W1, Q2, W2], ddq1_num.subs(param_subs), modules="numpy" )
ddq2_np = sp.lambdify([Q1, W1, Q2, W2], ddq2_num.subs(param_subs), modules="numpy" )
def ode_system (t_val, state ):
q1v, w1v, q2v, w2v = state
return [
w1v,
ddq1_np(q1v, w1v, q2v, w2v),
w2v,
ddq2_np(q1v, w1v, q2v, w2v),
]
y0 = [np.pi / 2 , 0.0 , np.pi / 4 , 0.0 ]
t_span = (0 , 20 )
t_eval = np.linspace(*t_span, 5000 )
sol_num = solve_ivp(ode_system, t_span, y0, method="DOP853" ,
t_eval=t_eval, rtol=1e-9 , atol=1e-11 )
q1_t = sol_num.y[0 ]
q2_t = sol_num.y[2 ]
x1_t = l_val * np.sin(q1_t)
y1_t = -l_val * np.cos(q1_t)
x2_t = x1_t + l_val * np.sin(q2_t)
y2_t = y1_t - l_val * np.cos(q2_t)
fig, axes = plt.subplots(1 , 2 , figsize=(13 , 5 ))
axes[0 ].plot(sol_num.t, q1_t, label="q₁(t)" , lw=1 )
axes[0 ].plot(sol_num.t, q2_t, label="q₂(t)" , lw=1 , ls="--" )
axes[0 ].set_xlabel("Time (s)" )
axes[0 ].set_ylabel("Angle (rad)" )
axes[0 ].legend()
axes[0 ].set_title("Double Pendulum Angles" )
axes[1 ].plot(x2_t, y2_t, lw=0.4 , alpha=0.7 , color="steelblue" )
axes[1 ].set_aspect("equal" )
axes[1 ].set_xlabel("x₂ (m)" )
axes[1 ].set_ylabel("y₂ (m)" )
axes[1 ].set_title("Trajectory of Bob 2" )
plt.tight_layout()
plt.show()
return sol_num
if __name__ == "__main__" :
double_pendulum_symbolic_to_numerical()
10. Complete Example B — Symbolic ODE Solve then Plot import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
def symbolic_then_numerical_ode ():
"""
(b) Solve a 2D system of ODEs symbolically, then compare with numerical solution.
System (predator-prey linearization about fixed point):
dx/dt = a*x - b*x*y → linearized: dx/dt = a*x - b*y
dy/dt = -c*y + d*x*y → linearized: dy/dt = d*x - c*y
"""
t = sp.Symbol("t" , positive=True )
x = sp.Function("x" )
y = sp.Function("y" )
a, b, c, d = sp.symbols("a b c d" , positive=True )
eq1 = sp.Eq(x(t).diff(t), a * x(t) - b * y(t))
eq2 = sp.Eq(y(t).diff(t), d * x(t) - c * y(t))
print ("Linearized Lotka-Volterra system:" )
print (" " , eq1)
print (" " , eq2)
system = [eq1, eq2]
sol = sp.dsolve(system, [x(t), y(t)])
print ("\nSymbolic solution:" )
for s in sol:
print (" " , s)
params = {a: 0.5 , b: 0.1 , c: 0.5 , d: 0.1 }
x_sol_sym = sol[0 ].rhs.subs(params)
y_sol_sym = sol[1 ].rhs.subs(params)
C1, C2, C3, C4 = sp.symbols("C1 C2 C3 C4" )
x0_val, y0_val = 10.0 , 5.0
ic_eqs = [
sp.Eq(x_sol_sym.subs(t, 0 ), x0_val),
sp.Eq(y_sol_sym.subs(t, 0 ), y0_val),
]
try :
const_sol = sp.solve(ic_eqs, [C1, C2, C3, C4])
x_particular = x_sol_sym.subs(const_sol)
y_particular = y_sol_sym.subs(const_sol)
x_np = sp.lambdify(t, x_particular, modules="numpy" )
y_np = sp.lambdify(t, y_particular, modules="numpy" )
t_arr = np.linspace(0 , 20 , 1000 )
x_sym_vals = x_np(t_arr)
y_sym_vals = y_np(t_arr)
has_symbolic = True
except Exception as e:
print (f"Symbolic IC application failed: {e} " )
has_symbolic = False
def lotka_volterra_linear (t_val, state ):
xv, yv = state
a_n, b_n, c_n, d_n = 0.5 , 0.1 , 0.5 , 0.1
return [a_n * xv - b_n * yv, d_n * xv - c_n * yv]
t_arr = np.linspace(0 , 20 , 1000 )
sol_num = solve_ivp(lotka_volterra_linear, (0 , 20 ), [x0_val, y0_val],
t_eval=t_arr, method="RK45" , rtol=1e-10 )
fig, axes = plt.subplots(1 , 2 , figsize=(13 , 5 ))
if has_symbolic:
axes[0 ].plot(t_arr, x_sym_vals.real, "r--" , lw=1.5 , label="x(t) symbolic" )
axes[0 ].plot(t_arr, y_sym_vals.real, "b--" , lw=1.5 , label="y(t) symbolic" )
axes[0 ].plot(sol_num.t, sol_num.y[0 ], "r-" , lw=1 , alpha=0.7 , label="x(t) numerical" )
axes[0 ].plot(sol_num.t, sol_num.y[1 ], "b-" , lw=1 , alpha=0.7 , label="y(t) numerical" )
axes[0 ].set_xlabel("Time" )
axes[0 ].set_ylabel("Population" )
axes[0 ].legend()
axes[0 ].set_title("Linearized Lotka-Volterra" )
axes[1 ].plot(sol_num.y[0 ], sol_num.y[1 ], color="steelblue" , lw=1.2 )
axes[1 ].scatter([x0_val], [y0_val], color="red" , zorder=5 , label="Initial state" )
axes[1 ].set_xlabel("x (prey)" )
axes[1 ].set_ylabel("y (predator)" )
axes[1 ].legend()
axes[1 ].set_title("Phase Portrait" )
plt.tight_layout()
plt.show()
if __name__ == "__main__" :
symbolic_then_numerical_ode()
Quick Reference Task SymPy API Define symbol with assumptions sp.Symbol("x", real=True, positive=True)Expand expression sp.expand(expr)Factor polynomial sp.factor(expr)Simplify (general) sp.simplify(expr)Trig simplification sp.trigsimp(expr)Differentiate sp.diff(expr, var, n)Integrate (indefinite) sp.integrate(expr, var)Integrate (definite) sp.integrate(expr, (var, a, b))Power series sp.series(expr, var, x0, n)Solve algebraic equation sp.solve(eq, var)Solve ODE sp.dsolve(ode, func, ics={})Matrix operations sp.Matrix([[...]]); M.det(), M.inv(), M.eigenvects()Convert to numpy sp.lambdify([vars], expr, modules="numpy")Common subexpression elimination sp.cse([expr1, expr2])LaTeX output sp.latex(expr)Pretty print sp.pprint(expr, use_unicode=True)
Performance Tips
Use sp.lambdify with modules="numpy" for vectorized numerical evaluation.
Apply sp.cse before lambdifying multi-output expressions to avoid redundant computation.
Set assumptions on symbols (e.g., real=True, positive=True) — SymPy simplifies more aggressively.
Use sp.nsimplify to convert floating-point results back to exact rational/algebraic form.
For large matrix operations, consider sp.MatrixSymbol for abstract manipulation.
