Beyond Simple Algebra
While simplifying polynomials is useful, SymPy’s real power for scientists lies in its ability to solve complex calculus problems and differential equations—tasks that usually require expensive software like Mathematica or Maple.
Taylor Series Expansions
In physics and engineering, we often approximate complex functions using Taylor series. SymPy can generate these to any arbitrary order.
Interactive Lab
import sympy as sp
x = sp.symbols('x')
# Taylor series of sin(x) around 0 up to 10th order
series = sp.series(sp.sin(x), x, 0, 10)
print(f"Taylor series of sin(x):\n{series}")
# Remove the 'O' term for numerical evaluation
func_approx = series.removeO()
print(f"\nPolynomial version: {func_approx}")
Expected output
Taylor series of sin(x):
x - x**3/6 + x**5/120 - x**7/5040 + x**9/362880 + O(x**10)
Polynomial version: x**9/362880 - x**7/5040 + x**5/120 - x**3/6 + x
1
2import sympy as sp
3x = sp.symbols('x')
4
5
6series = sp.series(sp.sin(x), x, 0, 10)
7print(f"Taylor series of sin(x):\n{series}")
8
9
10func_approx = series.removeO()
11print(f"\nPolynomial version: {func_approx}")
12
Solving Ordinary Differential Equations (ODEs)
SymPy can solve many classes of ODEs analytically. This is extremely useful for verifying numerical solvers or finding exact solutions for simple physical models.
Example: The Harmonic Oscillator
The equation is .
Interactive Lab
import sympy as sp
x = sp.symbols('x')
f = sp.Function('f')
omega = sp.symbols('omega', positive=True)
# Define the ODE
diffeq = sp.Eq(f(x).diff(x, x) + omega**2 * f(x), 0)
print(f"ODE: {diffeq}")
# Solve the ODE
solution = sp.dsolve(diffeq, f(x))
print(f"General Solution: {solution}")
Expected output
ODE: Eq(omega**2*f(x) + Derivative(f(x), (x, 2)), 0)
General Solution: Eq(f(x), C1*sin(omega*x) + C2*cos(omega*x))
1
2import sympy as sp
3x = sp.symbols('x')
4f = sp.Function('f')
5omega = sp.symbols('omega', positive=True)
6
7
8diffeq = sp.Eq(f(x).diff(x, x) + omega**2 * f(x), 0)
9print(f"ODE: {diffeq}")
10
11
12solution = sp.dsolve(diffeq, f(x))
13print(f"General Solution: {solution}")
14
Multivariable Calculus
SymPy handles gradients, Jacobians, and Hessians with ease.
Interactive Lab
import sympy as sp
x, y, z = sp.symbols('x y z')
f = x**2 + y**2 + z**2
# Partial derivatives
df_dx = sp.diff(f, x)
print(f"df/dx = {df_dx}")
# Gradient vector
grad = [sp.diff(f, var) for var in (x, y, z)]
print(f"Gradient: {grad}")
Expected output
df/dx = 2*x
Gradient: [2*x, 2*y, 2*z]
1
2import sympy as sp
3x, y, z = sp.symbols('x y z')
4
5f = x**2 + y**2 + z**2
6
7
8df_dx = sp.diff(f, x)
9print(f"df/dx = {df_dx}")
10
11
12grad = [sp.diff(f, var) for var in (x, y, z)]
13print(f"Gradient: {grad}")
14
Symbolic Linear Algebra
Calculating the eigenvalues or inverse of a matrix with symbols is a common requirement in theoretical research.
Interactive Lab
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)
# Determinant
print("\nDeterminant:", M.det())
# Inverse
print("\nInverse:\n", M.inv())
# Eigenvalues (returns a dict of value: multiplicity)
print("\nEigenvalues:", M.eigenvals())
Expected output
Matrix M:
Matrix([[a, b], [c, d]])
Determinant: a*d - b*c
Inverse:
Matrix([[d/(a*d - b*c), -b/(a*d - b*c)], [-c/(a*d - b*c), a/(a*d - b*c)]])
Eigenvalues: {a/2 + d/2 - sqrt(a**2 - 2*a*d + 4*b*c + d**2)/2: 1, a/2 + d/2 + sqrt(a**2 - 2*a*d + 4*b*c + d**2)/2: 1}
1
2import sympy as sp
3a, b, c, d = sp.symbols('a b c d')
4
5M = sp.Matrix([[a, b], [c, d]])
6print("Matrix M:\n", M)
7
8
9print("\nDeterminant:", M.det())
10
11
12print("\nInverse:\n", M.inv())
13
14
15print("\nEigenvalues:", M.eigenvals())
16
Global Optimization (Symbolic)
If a function is simple enough, we can find its global minimum by solving where the derivative is zero.
Interactive Lab
import sympy as sp
x = sp.symbols('x')
f = x**4 - 4*x**2 + 3
# Find critical points
critical_points = sp.solve(sp.diff(f, x), x)
print(f"Critical points: {critical_points}")
# Determine which is a minimum using second derivative
for point in critical_points:
second_deriv = sp.diff(f, x, 2).subs(x, point)
nature = "minimum" if second_deriv > 0 else "maximum"
print(f"x = {point} is a {nature} (f''={second_deriv})")
Expected output
Critical points: [0, -sqrt(2), sqrt(2)]
x = 0 is a maximum (f''=-8)
x = -sqrt(2) is a minimum (f''=16)
x = sqrt(2) is a minimum (f''=16)
1
2import sympy as sp
3x = sp.symbols('x')
4
5f = x**4 - 4*x**2 + 3
6
7
8critical_points = sp.solve(sp.diff(f, x), x)
9print(f"Critical points: {critical_points}")
10
11
12for point in critical_points:
13 second_deriv = sp.diff(f, x, 2).subs(x, point)
14 nature = "minimum" if second_deriv > 0 else "maximum"
15 print(f"x = {point} is a {nature} (f''={second_deriv})")
16
In the next section, we will integrate these symbolic results into numerical workflows using the SciPy stack.