import sympy as sp
# Define symbols with assumptions
m, g, x, y = sp.symbols('m, g, x, y', real=True, positive=True)
theta = sp.symbols('theta', real=True)5.1 Introduction
SymPy provides powerful tools for symbolic vector algebra and matrix manipulations. This section covers defining variables, vectors, computing norms, and working with rotation matrices, following the best practices for mechanics.
Variables and Expressions
It is good practice to state everything that is known about symbolic variables using “assumptions”. For example, in mechanics, mass and gravity are often real and positive.
Vectors
Vectors are defined using sp.Matrix.
# Define a symbolic vector with known magnitude and direction
# Example: Weight vector pointing down in Y
WW = m * g * sp.Matrix([0, -1, 0])
WW\(\displaystyle \left[\begin{matrix}0\\- g m\\0\end{matrix}\right]\)
You can also define vectors with unknown components:
F_Ax, F_Ay = sp.symbols('F_Ax, F_Ay', real=True)
FFA = sp.Matrix([F_Ax, F_Ay, 0])
FFA\(\displaystyle \left[\begin{matrix}F_{Ax}\\F_{Ay}\\0\end{matrix}\right]\)
Vector Operations
Standard operations like addition, subtraction, cross products, and dot products are straightforward.
Addition
# Sum of forces
RR = WW + FFA
RR\(\displaystyle \left[\begin{matrix}F_{Ax}\\F_{Ay} - g m\\0\end{matrix}\right]\)
Cross Product
The cross product \(\mathbb{A} \times \mathbb{B}\) is computed using .cross():
# Moment MM = rr x FFA
rr = sp.Matrix([1, 0, 0]) # Position vector
MM = rr.cross(FFA)
MM\(\displaystyle \left[\begin{matrix}0\\0\\F_{Ay}\end{matrix}\right]\)
Dot Product
The dot product \(\mathbb{A} \cdot \mathbb{B}\) is computed using .dot():
# Power P = F . v
vv = sp.Matrix([1, 1, 0])
P = FFA.dot(vv)
P\(\displaystyle F_{Ax} + F_{Ay}\)
Norm and Normalization
The norm (magnitude) \(||\mathbb{A}||\) is computed using .norm():
# Magnitude of force FFA
FFA_mag = FFA.norm()
FFA_mag\(\displaystyle \sqrt{F_{Ax}^{2} + F_{Ay}^{2}}\)
To get the unit vector (direction), use .normalized():
# Direction of FFA
ee_A = FFA.normalized()
ee_A\(\displaystyle \left[\begin{matrix}\frac{F_{Ax}}{\sqrt{F_{Ax}^{2} + F_{Ay}^{2}}}\\\frac{F_{Ay}}{\sqrt{F_{Ax}^{2} + F_{Ay}^{2}}}\\0\end{matrix}\right]\)
Matrices
SymPy has its own Matrix class for symbolic linear algebra.
M = sp.Matrix([[1, x], [y, 1]])
print(f"Matrix M:\n{M}")
# Determinant
det = M.det()
print(f"Determinant: {det}")
# Inverse
inv = M.inv()
print(f"Inverse:\n{inv}")Matrix M:
Matrix([[1, x], [y, 1]])
Determinant: -x*y + 1
Inverse:
Matrix([[-1/(x*y - 1), x/(x*y - 1)], [y/(x*y - 1), -1/(x*y - 1)]])Rotation Matrices
SymPy offers helper functions for rotation matrices: sp.rot_ccw_axis1 (X-axis), sp.rot_ccw_axis2 (Y-axis), and sp.rot_ccw_axis3 (Z-axis).
# Rotation matrix around Z-axis
Rz = sp.rot_ccw_axis3(theta)
Rz\(\displaystyle \left[\begin{matrix}\cos{\left(\theta \right)} & - \sin{\left(\theta \right)} & 0\\\sin{\left(\theta \right)} & \cos{\left(\theta \right)} & 0\\0 & 0 & 1\end{matrix}\right]\)
We can rotate a vector by matrix multiplication using the @ operator. Rotational matrices are applied on the left side.
# Rotate FFA around Z-axis
FFA_rotated = Rz @ FFA
FFA_rotated\(\displaystyle \left[\begin{matrix}F_{Ax} \cos{\left(\theta \right)} - F_{Ay} \sin{\left(\theta \right)}\\F_{Ax} \sin{\left(\theta \right)} + F_{Ay} \cos{\left(\theta \right)}\\0\end{matrix}\right]\)
Projection
The vector projection of \(\mathbb{A}\) along \(\mathbb{B}\) can be computed using .project().
AA = 20 * sp.Matrix([sp.cos(sp.pi/2), sp.sin(sp.pi/2), 0])
BB = sp.Matrix([sp.sqrt(2), sp.sqrt(3), 0])
# Vector projection of A along B
AA_proj_BB = AA.project(BB)
AA_proj_BB\(\displaystyle \left[\begin{matrix}4 \sqrt{6}\\12\\0\end{matrix}\right]\)
Substitution
You can substitute values or other variables into expressions.
# Define an expression
expr = x**2 + y**2
# Substitute x = 2, y = 3
value = expr.subs({x: 2, y: 3})
print(f"Value at x=2, y=3: {value}")Value at x=2, y=3: 13Solving Equations
SymPy can solve algebraic equations symbolically.
# Solve x^2 - 4 = 0
solution = sp.solve(x**2 - 4, x)
print(f"Solutions for x: {solution}")Solutions for x: [2]Solving Systems of Equations
SymPy’s sp.solve can handle systems of linear and non-linear equations. It is best practice to explicitly define equations using sp.Eq(LHS, RHS). If you pass an expression without sp.Eq, SymPy assumes it equals zero.
import sympy as sp
# Define symbols
x, y, F = sp.symbols('x y F', real=True)
# Define equations
# Eq 1: x + y + F = 0
# Eq 2: x - 2*y = 5
eq1 = sp.Eq(x + y + F, 0)
eq2 = sp.Eq(x - 2*y, 5)
# Solve for x and y
sol = sp.solve([eq1, eq2], [x, y])
x_sol = sol[x]
x_sol\(\displaystyle \frac{5}{3} - \frac{2 F}{3}\)
Differentiation
# Derivative of sin(x) * e^x
f = sp.sin(x) * sp.exp(x)
df_dx = sp.diff(f, x)
print(f"Derivative: {df_dx}")Derivative: exp(x)*sin(x) + exp(x)*cos(x)Integration
# Indefinite integral of cos(x)
integral = sp.integrate(sp.cos(x), x)
print(f"Indefinite integral: {integral}")
# Definite integral of x^2 from 0 to 1
def_integral = sp.integrate(x**2, (x, 0, 1))
print(f"Definite integral: {def_integral}")Indefinite integral: sin(x)
Definite integral: 1/3Differential Equations
Use sp.dsolve to find analytical solutions to ordinary differential equations (ODEs).
# Define function and derivative
t = sp.symbols('t', real=True)
x = sp.Function('x')(t)
k = sp.symbols('k', positive=True)
v0 = sp.symbols('v0', real=True)
# Define ODE: x'' = -k^2 * x (Harmonic Oscillator)
ode = sp.Eq(x.diff(t, 2), -k**2 * x)
# Initial conditions: x(0) = 0, x'(0) = v0
ics = {x.subs(t, 0): 0, x.diff(t).subs(t, 0): v0}
# Solve
sol_ode = sp.dsolve(ode, ics=ics)
sol_ode\(\displaystyle x{\left(t \right)} = \frac{v_{0} \sin{\left(k t \right)}}{k}\)
Converting to NumPy (Lambdify)
For numerical evaluation of symbolic expressions (e.g., for plotting or iterative solvers), use lambdify.
import numpy as np
# Convert symbolic solution for x to a Python function
# The first argument is a list of variables the function will accept
x_func = sp.lambdify([F], x_sol, 'numpy')
# Evaluate for a specific value
x_val = x_func(10)
x_val-4.999999999999999When lambdifying a vector (Matrix), the result is often a list of lists or a 2D array. To get a standard 1D NumPy array for vector operations, use .flatten().
# Define a symbolic vector
m, g = sp.symbols('m g', real=True)
WW = m * g * sp.Matrix([0, -1, 0])
# Lambdify
WW_func = sp.lambdify([m, g], WW, 'numpy')
# Evaluate and flatten
WW_np = WW_func(20, 9.81).flatten()
WW_nparray([ 0. , -196.2, 0. ])