Python for Engineering Mathematics

Python with NumPy, SymPy, SciPy and Matplotlib is the free standard toolkit for engineering computation. Learn how symbolic (exact) and numerical (floating-point) computing differ, core NumPy arrays and vectorisation, linear algebra, SymPy calculus and ODE solving, floating-point pitfalls, and a clean workflow for engineering scripts.

Hand calculation teaches understanding; computers deliver answers at scale. Python lets you do both: SymPy manipulates formulas exactly like a textbook, while NumPy crunches millions of numbers per second for simulations and data.

Knowing which tool to use — and knowing floating-point numbers are approximations — separates reliable engineering code from code that quietly gives wrong answers.

An analogy

SymPy is the mathematician with pencil and paper — exact but slow; NumPy is the calculator with a thousand hands — fast but rounding every result to about 16 digits.

The scientific Python stack

Each library has a clear job:

Vectorise, don't loop

A Python for-loop over a million elements is roughly 50–100× slower than one NumPy call, because NumPy runs compiled C loops. Write formulas on whole arrays: stress = F / A for arrays F and A. Use broadcasting to combine arrays of different shapes without copies.

Floating-point realities

0.1 + 0.2 == 0.3 is False in every language using IEEE-754 doubles. Compare with a tolerance (np.isclose). Subtracting nearly equal numbers loses digits (catastrophic cancellation) — which is why a finite-difference step that is too small gives worse derivatives, as the plot shows.

Reliable engineering code

Keep units consistent (or use a units library like pint), test against analytic cases, use solve rather than explicit inverses, prefer library routines (tested and optimised) over hand-written ones, fix random seeds for reproducibility, and put calculations in functions with docstrings. SymPy's lambdify turns derived formulas into fast NumPy functions — derive symbolically, evaluate numerically.

Formulas

Worked example: Solving a linear system

Given: A = [[4,−1,0],[−1,4,−1],[0,−1,4]], b = [2,4,10]

Find: x.

Answer: x = [1, 2, 3]

Worked example: Symbolic derivative

Given: f = x³e^(−x)

Find: f′ and its zeros.

Answer: Maximum at x = 3

Worked example: Floating-point comparison

Given: 0.1 + 0.2

Find: Why is it not equal to 0.3?

Answer: Use np.isclose(a, b) or abs(a − b) < tol

Symbolic and numeric side by side

import numpy as np
import sympy as sp

# --- Symbolic: SymPy does exact algebra and calculus ---
x = sp.symbols("x")
f = x**3 * sp.exp(-x)
print("f'(x)       =", sp.simplify(sp.diff(f, x)))
print("int_0^oo f  =", sp.integrate(f, (x, 0, sp.oo)))
y = sp.Function("y")
ode = sp.dsolve(y(x).diff(x, 2) + 4*y(x), y(x), ics={y(0): 1, y(x).diff(x).subs(x, 0): 0})
print("ODE         :", ode)

# --- Numerical: NumPy does fast arrays and linear algebra ---
A = np.array([[4.0, -1, 0], [-1, 4, -1], [0, -1, 4]])
b = np.array([2.0, 4, 10])
print("solve Ax=b  =", np.linalg.solve(A, b))
print("eigenvalues =", np.round(np.linalg.eigvalsh(A), 4))
t = np.linspace(0, 1, 5)
print("vectorised  =", np.round(np.sin(np.pi * t), 4))
# f'(x)       = x**2*(3 - x)*exp(-x)
# int_0^oo f  = 6
# ODE         : Eq(y(x), cos(2*x))
# solve Ax=b  = [1. 2. 3.]
# eigenvalues = [2.5858 4.     5.4142]
# vectorised  = [0.     0.7071 1.     0.7071 0.    ]

SymPy returns exact objects (6, cos(2x)); NumPy returns floating-point arrays. The integral 6 = Γ(4) matches the special-functions lesson — a good habit is to cross-check tools against known results.

Common mistakes

Applications

Key takeaways

Going further