Computer algebra systems: SymPy, SageMath, and Mathematica for research
A skill for using computer algebra systems (CAS) in mathematical research. Covers symbolic differentiation, integration, equation solving, series expansion, linear algebra, and polynomial arithmetic using SymPy, SageMath, and Mathematica, with practical workflows for research mathematics.
from sympy import (
symbols, expand, factor, simplify, cancel, apart,
sin, cos, exp, log, sqrt, pi, oo, I,
Rational, Eq, solve, solveset, S
)
x, y, z, t, n, k = symbols("x y z t n k")
a, b, c = symbols("a b c", real=True)
# Expression manipulation
expr = (x + 1) ** 3
expanded = expand(expr) # x**3 + 3*x**2 + 3*x + 1
factored = factor(expanded) # (x + 1)**3
# Trigonometric simplification
from sympy import trigsimp
trig_expr = sin(x)**2 + cos(x)**2
simplified = trigsimp(trig_expr) # 1
# Partial fraction decomposition
rational = (x**2 + 2*x + 3) / ((x + 1) * (x + 2) * (x + 3))
partial = apart(rational, x)
# 3/(2*(x + 3)) - 2/(x + 2) + 1/(2*(x + 1))
from sympy import diff, integrate, limit, series, Sum, Product
# Differentiation
f = x**3 * exp(-x) * sin(x)
f_prime = diff(f, x)
f_double_prime = diff(f, x, 2)
# Integration
# Definite integral
area = integrate(exp(-x**2), (x, -oo, oo)) # sqrt(pi)
# Indefinite integral
antideriv = integrate(x * sin(x), x) # -x*cos(x) + sin(x)
# Limits
lim_result = limit(sin(x) / x, x, 0) # 1
lim_inf = limit((1 + 1/n)**n, n, oo) # E (Euler's number)
# Taylor series
taylor = series(exp(x) * cos(x), x, 0, n=6)
# 1 + x - x**3/3 - x**4/6 + ...
# Summation
harmonic = Sum(1/k, (k, 1, n))
partial_sum = harmonic.doit() # harmonic(n) -- returns harmonic number
geometric = Sum(x**k, (k, 0, oo))
closed_form = geometric.doit() # Piecewise(1/(1 - x), Abs(x) < 1)
# Algebraic equations
solutions = solve(x**3 - 6*x**2 + 11*x - 6, x) # [1, 2, 3]
# System of equations
system_sol = solve([
2*x + 3*y - 7,
x - y + 1
], [x, y]) # {x: 4/5, y: 9/5}
# Differential equations
from sympy import Function, dsolve, Derivative
f = Function("f")
# f''(x) + f(x) = 0 (simple harmonic oscillator)
ode = Eq(f(x).diff(x, 2) + f(x), 0)
general_solution = dsolve(ode, f(x))
# f(x) = C1*sin(x) + C2*cos(x)
# With initial conditions
particular = dsolve(ode, f(x), ics={f(0): 1, f(x).diff(x).subs(x, 0): 0})
# f(x) = cos(x)
from sympy import Matrix, eye, zeros, det, Rational
# Define a symbolic matrix
A = Matrix([
[1, 2, 3],
[4, 5, 6],
[7, 8, 10]
])
# Basic operations
print(f"Determinant: {det(A)}") # -3
print(f"Inverse:\n{A.inv()}")
print(f"Eigenvalues: {A.eigenvals()}")
print(f"Rank: {A.rank()}")
# Characteristic polynomial
lam = symbols("lambda")
char_poly = (A - lam * eye(3)).det()
char_poly = expand(char_poly)
# Jordan normal form
P, J = A.jordan_form()
# Null space and column space
null = A.nullspace()
col_space = A.columnspace()
# Symbolic matrix with parameters
M = Matrix([
[a, b],
[c, a]
])
eigenvals = M.eigenvals() # {a - sqrt(b*c): 1, a + sqrt(b*c): 1}
# SageMath syntax (Python-based, but with enhanced number theory)
# Run in SageMath environment or via sage -python
"""
# Prime factorization
factor(2024) # 2^3 * 11 * 23
# Modular arithmetic
R = IntegerModRing(17)
R(3)^(-1) # multiplicative inverse of 3 mod 17
# Elliptic curves
E = EllipticCurve(QQ, [-1, 0])
E.rank()
E.torsion_subgroup()
E.gens()
# Polynomial rings
R.<x,y> = PolynomialRing(QQ)
I = R.ideal(x^2 + y^2 - 1, x - y)
I.groebner_basis() # [y^2 - 1/2, x - y]
# Group theory
G = SymmetricGroup(4)
G.order() # 24
G.center()
G.normal_subgroups()
"""
"""
# SageMath combinatorics
Partitions(10).cardinality() # 42
# Graph theory
G = graphs.PetersenGraph()
G.chromatic_number() # 3
G.is_vertex_transitive() # True
G.automorphism_group().order() # 120
# Posets and lattices
P = posets.BooleanLattice(3)
P.is_lattice()
P.mobius_function(P.bottom(), P.top())
"""
(* Symbolic integration *)
Integrate[x^n * Exp[-x], {x, 0, Infinity}, Assumptions -> n > -1]
(* Result: Gamma[1 + n] *)
(* Solve a PDE *)
DSolve[D[u[x, t], t] == k * D[u[x, t], {x, 2}], u[x, t], {x, t}]
(* Asymptotic expansion *)
Series[Gamma[n + 1], {n, Infinity, 3}]
(* Minimize with constraints *)
NMinimize[{x^2 + y^2, x + y >= 1}, {x, y}]
(* Compute a sum in closed form *)
Sum[1/k^2, {k, 1, Infinity}] (* Pi^2/6 *)
Common CAS workflow in mathematical research:
from sympy import simplify, Abs
def verify_identity(lhs, rhs):
"""Verify a proposed mathematical identity symbolically."""
diff = simplify(lhs - rhs)
if diff == 0:
return "VERIFIED: identity holds symbolically"
else:
return f"NOT VERIFIED: difference = {diff}"
# Example: verify Cauchy-Schwarz for 2D
a1, a2, b1, b2 = symbols("a1 a2 b1 b2", real=True)
lhs = (a1*b1 + a2*b2)**2
rhs = (a1**2 + a2**2) * (b1**2 + b2**2)
diff = expand(rhs - lhs)
# (a1*b2 - a2*b1)**2 >= 0, confirming Cauchy-Schwarz