Appendix D — NumPy / SciPy / SymPy Cheatsheet

This appendix is a task-organized quick reference for the four libraries used throughout the book: SymPy (symbolic math — exact derivatives, integrals, limits, series), NumPy (numerical arrays — Riemann sums, finite differences), SciPy (numerical integration, ODEs, optimization), and Matplotlib (every plot). It follows the conventions in the continuity tracker §4: SymPy for symbolic work, NumPy for numerical work, SciPy for integration/ODEs/optimization, Matplotlib for plotting, type hints on functions, and results shown in # → comments so you can verify without running anything.

Install instructions live in Appendix C — they are not repeated here. Each recipe lists the chapters that lean on it most heavily.

1. Setup and Imports

The header that opens essentially every notebook in this book.

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
from scipy import integrate, optimize, linalg

# SymPy pretty-printing in Jupyter:
sp.init_printing()

# Reproducible randomness (used in Monte Carlo / data-science tracks):
rng = np.random.default_rng(42)

Computational Note — Keep symbolic and numeric worlds separate. SymPy objects are exact (sp.Rational(1, 3) is one-third forever); NumPy arrays are floating point (1/3 → 0.333...). Convert from symbolic to numeric only at the end, with sp.lambdify (Recipe 2.10).

2. SymPy — Symbolic Calculus

2.1 Declaring symbols and expressions

Used in: every symbolic chapter (Ch. 3, 6, 7, 12).

x, y, t = sp.symbols('x y t', real=True)
n = sp.symbols('n', integer=True, positive=True)
f = x**3 + sp.sin(x)        # an expression, not a Python function

2.2 Limits — sp.limit

Used in: Ch. 3 (the limit), Ch. 4 (continuity), Ch. 17 (improper integrals), Ch. 20–22 (sequences/series).

sp.limit(sp.sin(x)/x, x, 0)            # → 1
sp.limit((1 + 1/x)**x, x, sp.oo)       # → E   (Euler's number e)
sp.limit(1/x, x, 0, '+')               # → oo  (one-sided, from the right)
sp.limit(1/x, x, 0, '-')               # → -oo (one-sided, from the left)

The one-sided forms compute $\lim_{x \to 0^+}$ and $\lim_{x \to 0^-}$.

2.3 Derivatives — sp.diff

Used in: Ch. 6 (the derivative), Ch. 7 (rules), Ch. 8 (implicit/related rates).

sp.diff(x**3 + sp.sin(x), x)        # → 3*x**2 + cos(x)
sp.diff(sp.exp(x**2), x)            # → 2*x*exp(x**2)   (chain rule)
sp.diff(x**4, x, 3)                 # → 24*x   (third derivative)
sp.diff(x**4, x, x, x)              # → 24*x   (equivalent form)

This computes $\frac{d}{dx}(x^3 + \sin x) = 3x^2 + \cos x$ and the third derivative $\frac{d^3}{dx^3}x^4 = 24x$.

2.4 Partial and mixed derivatives

Used in: Ch. 29 (several variables), Ch. 30 (gradient/chain rule).

g = sp.sin(x*y) + x**2 * y
sp.diff(g, x)            # → 2*x*y + y*cos(x*y)        (∂g/∂x)
sp.diff(g, y)            # → x**2 + x*cos(x*y)         (∂g/∂y)
sp.diff(g, x, y)         # → 2*x - x*y*sin(x*y) + cos(x*y)  (mixed ∂²g/∂x∂y)

These are $\frac{\partial g}{\partial x}$, $\frac{\partial g}{\partial y}$, and the mixed partial $\frac{\partial^2 g}{\partial x\,\partial y}$.

2.5 Indefinite and definite integrals — sp.integrate

Used in: Ch. 12 (antiderivatives), Ch. 13 (definite integral), Ch. 14 (FTC), Ch. 15–16 (techniques).

sp.integrate(x**2, x)               # → x**3/3      (indefinite; SymPy omits +C)
sp.integrate(sp.cos(x), x)          # → sin(x)
sp.integrate(x**2, (x, 0, 1))       # → 1/3         (definite ∫₀¹ x² dx)
sp.integrate(sp.exp(-x**2), (x, -sp.oo, sp.oo))  # → sqrt(pi)  (Gaussian)

The last computes $\int_{-\infty}^{\infty} e^{-x^2}\,dx = \sqrt{\pi}$ — the basis of the normal-curve anchor example.

Common Pitfall — SymPy returns the indefinite integral without the +C. Add the constant yourself when reporting an antiderivative.

2.6 Series and Taylor expansions — sp.series

Used in: Ch. 23 (power/Taylor series), Ch. 24 (applications, Euler's formula).

sp.series(sp.exp(x), x, 0, 5)
# → 1 + x + x**2/2 + x**3/6 + x**4/24 + O(x**5)

sp.series(sp.sin(x), x, 0, 8)
# → x - x**3/6 + x**5/120 - x**7/5040 + O(x**8)

# Drop the O() term to get a usable polynomial:
sp.series(sp.cos(x), x, 0, 7).removeO()
# → x**6/720 - x**4/24 + x**2/2 + 1

This is the Taylor expansion about $x = 0$: $e^x = 1 + x + \frac{x^2}{2} + \frac{x^3}{6} + \cdots$.

2.7 Solving equations — sp.solve and sp.dsolve

Used in: Ch. 10–11 (optimization, Newton's method), Ch. 19 (differential equations).

sp.solve(x**2 - 2, x)               # → [-sqrt(2), sqrt(2)]
sp.solve([x + y - 1, x - y - 3], [x, y])   # → {x: 2, y: -1}

# Symbolic ODE: y' + y = 0
y_ = sp.Function('y')
sp.dsolve(sp.Eq(y_(t).diff(t) + y_(t), 0), y_(t))
# → Eq(y(t), C1*exp(-t))

2.8 Simplify, expand, factor

Used in: throughout, especially Ch. 7 and Ch. 15 (cleaning up by-hand results).

sp.simplify(sp.sin(x)**2 + sp.cos(x)**2)   # → 1
sp.expand((x + 1)**3)                       # → x**3 + 3*x**2 + 3*x + 1
sp.factor(x**3 + 3*x**2 + 3*x + 1)          # → (x + 1)**3
sp.trigsimp(2*sp.sin(x)*sp.cos(x))          # → sin(2*x)

2.9 Substitution and exact evaluation

Used in: every chapter that checks a value.

expr = x**2 + 1
expr.subs(x, 3)                 # → 10
expr.subs(x, sp.pi)             # → 1 + pi**2   (kept exact)
float(expr.subs(x, sp.pi))      # → 10.869604401089358
sp.N(sp.pi, 20)                 # → 3.1415926535897932385  (20 digits)

2.10 SymPy → NumPy bridge — sp.lambdify

Used in: every chapter that plots a symbolic result (Ch. 5, 9, 13, 23, 30).

f_expr = x**2 * sp.sin(x)
f = sp.lambdify(x, f_expr, 'numpy')   # fast NumPy-callable function
xs = np.linspace(0, np.pi, 5)
f(xs)
# → array([0.        , 0.43622... , 2.46740... , 3.92568... , 0.        ])

The Key Insight — lambdify is the seam between the exact and the numerical halves of every workflow: differentiate or integrate symbolically, then lambdify to evaluate and plot over a NumPy grid.

2.11 Matrices, gradient, and Jacobian

Used in: Ch. 30 (gradient), Ch. 31 (multivariable optimization), Ch. 33 (Jacobians, change of variables).

F = sp.Matrix([x**2 * y, sp.sin(x + y)])
v = sp.Matrix([x, y])

F.jacobian(v)
# → Matrix([[2*x*y,        x**2     ],
#           [cos(x + y),   cos(x + y)]])

phi = x**2 + y**2
grad = sp.Matrix([phi]).jacobian(v).T     # gradient ∇φ as a column
# → Matrix([[2*x], [2*y]])

H = sp.hessian(phi, (x, y))               # Hessian (second-derivative matrix)
# → Matrix([[2, 0], [0, 2]])

The Jacobian $J$ has entries $J_{ij} = \partial F_i / \partial x_j$; the gradient is $\nabla \phi = (\partial\phi/\partial x,\ \partial\phi/\partial y)$.

3. NumPy — Numerical Arrays

3.1 Building grids — linspace, arange

Used in: every plotting and Riemann-sum recipe.

np.linspace(0, 1, 5)        # → array([0.  , 0.25, 0.5 , 0.75, 1.  ])  (includes endpoint)
np.arange(0, 1, 0.25)       # → array([0.  , 0.25, 0.5 , 0.75])       (excludes endpoint)

3.2 Vectorized functions

Used in: everywhere a function is evaluated on a grid.

x = np.linspace(0, np.pi, 4)
np.sin(x)        # → array([0.        , 0.8660254, 0.8660254, 0.        ])
np.exp([0, 1])   # → array([1.        , 2.71828183])

3.3 Finite-difference derivatives — np.gradient

Used in: Ch. 5 (rates of change), Ch. 6 (the derivative numerically), Ch. 11 (linear approximation).

x = np.linspace(0, np.pi, 100)
y = np.sin(x)
dydx = np.gradient(y, x)        # central differences; same length as x
# dydx ≈ cos(x); e.g. dydx[0] ≈ 1.0, dydx[-1] ≈ -1.0

np.gradient approximates $\frac{dy}{dx}$ using central differences in the interior and one-sided differences at the ends.

3.4 Riemann sums — np.sum

Used in: Ch. 13 (the definite integral — the central numerical recipe of Part III).

def riemann(f, a: float, b: float, n: int, kind: str = "mid") -> float:
    """Riemann sum of f on [a, b] with n subintervals."""
    dx = (b - a) / n
    if kind == "left":
        x = a + dx * np.arange(n)
    elif kind == "right":
        x = a + dx * np.arange(1, n + 1)
    else:  # midpoint
        x = a + dx * (np.arange(n) + 0.5)
    return np.sum(f(x)) * dx

riemann(lambda x: x**2, 0, 1, 1000, "mid")   # → 0.33333325   (true value 1/3)

This approximates $\int_a^b f(x)\,dx \approx \sum_{i} f(x_i)\,\Delta x$.

3.5 Accumulated area — np.cumsum

Used in: Ch. 14 (FTC — the integral as an accumulation function $F(x) = \int_a^x f$).

x = np.linspace(0, np.pi, 1000)
dx = x[1] - x[0]
F = np.cumsum(np.sin(x)) * dx       # running integral of sin from 0 to x
# F[-1] ≈ 2.0   (since ∫₀^π sin x dx = 2); F ≈ 1 - cos(x)

3.6 Curve fitting — np.polyfit

Used in: Ch. 2 (functions/models), Ch. 9 (applications), Data Science track.

xs = np.array([0, 1, 2, 3])
ys = np.array([1, 3, 7, 13])              # roughly x² + x + 1
coeffs = np.polyfit(xs, ys, 2)            # fit a degree-2 polynomial
# → array([1., 1., 1.])   (a=1, b=1, c=1)
p = np.poly1d(coeffs)
p(4)                                       # → 21.0

polyfit returns coefficients highest-degree-first: here $p(x) = x^2 + x + 1$.

4. SciPy — Numerical Integration, ODEs, Optimization

4.1 1-D integration — integrate.quad

Used in: Ch. 13, 16, 17 (especially improper integrals), Ch. 18 (applications).

val, err = integrate.quad(lambda x: np.exp(-x**2), 0, np.inf)
# val → 0.8862269254527579   (= √π / 2),  err → ~7e-09

quad returns the value and an error estimate; it handles infinite limits.

4.2 Double and triple integrals — dblquad, tplquad

Used in: Ch. 32 (multiple integrals), Ch. 36 (surface integrals).

# ∫₀¹ ∫₀¹ (x + y) dy dx   — note SciPy's argument order is f(y, x)
val, err = integrate.dblquad(lambda y, x: x + y, 0, 1, 0, 1)
# val → 1.0

# ∫₀¹ ∫₀¹ ∫₀¹ xyz dz dy dx — f(z, y, x), innermost variable first
val3, _ = integrate.tplquad(lambda z, y, x: x*y*z, 0, 1, 0, 1, 0, 1)
# val3 → 0.125   (= 1/8)

Common Pitfall — In dblquad/tplquad the integrand's arguments are ordered innermost variable first (f(y, x), f(z, y, x)), which is the reverse of how you read the integral. The inner bounds may be functions of the outer variables.

4.3 Solving ODEs — integrate.solve_ivp (modern API)

Used in: Ch. 19 (differential equations, SIR model), Ch. 39 (capstone modeling).

def decay(t: float, y: np.ndarray) -> list[float]:
    """dy/dt = -0.5 y."""
    return [-0.5 * y[0]]

sol = integrate.solve_ivp(
    decay, t_span=(0, 10), y0=[1.0],
    t_eval=np.linspace(0, 10, 100), method="RK45",
)
sol.y[0, -1]      # → 0.006737...  (≈ e^(-5), the exact solution at t=10)

Computational Note — Prefer solve_ivp (returns sol.t, sol.y) over the older odeint. The signature is f(t, y) for solve_ivp but f(y, t) for odeint — a frequent source of bugs.

4.4 The SIR system (anchor example)

Used in: Ch. 19 (introduction), Ch. 39 (capstone). Models susceptible $S$, infected $I$, recovered $R$.

def sir(t: float, y: np.ndarray, beta: float, gamma: float) -> list[float]:
    S, I, R = y
    N = S + I + R
    return [-beta*S*I/N, beta*S*I/N - gamma*I, gamma*I]

sol = integrate.solve_ivp(
    sir, (0, 160), [999, 1, 0], args=(0.3, 0.1),
    t_eval=np.linspace(0, 160, 320),
)
# sol.y -> shape (3, 320): rows are S(t), I(t), R(t)
# Infected peaks near t ≈ 45-55 for these parameters.

The system is $\dot S = -\beta SI/N,\ \dot I = \beta SI/N - \gamma I,\ \dot R = \gamma I$.

4.5 Root finding — optimize.fsolve and brentq

Used in: Ch. 11 (Newton's method), Ch. 10 (optimization critical points).

optimize.brentq(lambda x: x**3 - 2, 0, 2)          # → 1.2599210498948732  (∛2)
optimize.fsolve(lambda x: x**3 - 2, x0=1.0)        # → array([1.25992105])

# System of equations:
def system(v):
    x, y = v
    return [x**2 + y**2 - 1, x - y]
optimize.fsolve(system, [1, 1])                     # → array([0.70710678, 0.70710678])

4.6 Minimization — optimize.minimize

Used in: Ch. 10 (optimization), Ch. 31 (several variables), Ch. 30 (gradient descent anchor).

res = optimize.minimize(lambda v: (v[0]-3)**2 + (v[1]-4)**2, x0=[0, 0])
res.x        # → array([3., 4.])   (the minimizer)
res.fun      # → ~0.0              (the minimum value)

4.7 Constrained minimization

Used in: Ch. 31 (Lagrange multipliers, optimization under constraints), Economics track.

# Minimize x² + y² subject to x + y = 1  (expect x = y = 0.5)
cons = ({"type": "eq", "fun": lambda v: v[0] + v[1] - 1},)
res = optimize.minimize(lambda v: v[0]**2 + v[1]**2, [0, 0],
                        constraints=cons, method="SLSQP")
res.x        # → array([0.5, 0.5])
res.fun      # → 0.5

This solves $\min(x^2+y^2)$ subject to $x+y=1$ — the same problem Lagrange multipliers solve by hand in Ch. 31.

4.8 Linear algebra — scipy.linalg

Used in: Ch. 30–33 (gradients, Hessians, change of variables/Jacobian determinants).

A = np.array([[2., 1.], [1., 3.]])
b = np.array([1., 2.])
linalg.solve(A, b)        # → array([0.2, 0.6])   (solves A x = b)
linalg.det(A)             # → 5.0
vals, vecs = linalg.eig(A)
# vals → array([1.382+0.j, 3.618+0.j])  (eigenvalues; used for definiteness tests)

The determinant is exactly the Jacobian factor in change-of-variables integrals (Ch. 33).

5. Matplotlib — Visualization

5.1 Line plots

Used in: every chapter.

x = np.linspace(-5, 5, 400)
plt.plot(x, np.sin(x), label="sin x")
plt.plot(x, np.cos(x), "--", label="cos x")
plt.xlabel("x"); plt.ylabel("y")
plt.legend(); plt.grid(True)
plt.show()

5.2 Shading areas — fill_between

Used in: Ch. 13 (area under a curve), Ch. 18 (area between curves), normal-curve anchor.

x = np.linspace(0, 2, 400)
y = x**2
plt.plot(x, y)
plt.fill_between(x, y, where=(x >= 0.5) & (x <= 1.5),
                 alpha=0.3, color="tab:blue")     # shades ∫₀.₅^1.5 x² dx
plt.show()

5.3 Contour plots — contour, contourf

Used in: Ch. 29 (level curves), Ch. 30 (gradient ⟂ level curves), Ch. 31 (optimization landscapes).

xs = np.linspace(-2, 2, 200)
ys = np.linspace(-2, 2, 200)
X, Y = np.meshgrid(xs, ys)
Z = X**2 + Y**2
cs = plt.contour(X, Y, Z, levels=[0.5, 1, 2, 3])   # level curves
plt.clabel(cs)
plt.contourf(X, Y, Z, levels=20, cmap="viridis")    # filled version
plt.colorbar()
plt.show()

np.meshgrid turns two 1-D axes into the 2-D coordinate arrays every surface/contour/field plot needs.

5.4 3-D surfaces — plot_surface

Used in: Ch. 29 (functions of several variables), Ch. 32 (volume under a surface).

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")   # modern API; no Axes3D import needed
X, Y = np.meshgrid(np.linspace(-2, 2, 60), np.linspace(-2, 2, 60))
Z = np.sin(np.sqrt(X**2 + Y**2))
ax.plot_surface(X, Y, Z, cmap="viridis")
ax.set_xlabel("x"); ax.set_ylabel("y"); ax.set_zlabel("z")
plt.show()

5.5 Vector fields — quiver

Used in: Ch. 34 (vector fields), Ch. 35 (line integrals), Ch. 37 (flux).

xs = np.linspace(-2, 2, 15)
X, Y = np.meshgrid(xs, xs)
U, V = -Y, X                       # rotational field F = (-y, x)
plt.quiver(X, Y, U, V)
plt.gca().set_aspect("equal")
plt.show()

5.6 Streamlines — streamplot

Used in: Ch. 34 (flow visualization), Ch. 37 (divergence/curl intuition).

xs = np.linspace(-2, 2, 40)
X, Y = np.meshgrid(xs, xs)
U, V = -Y, X
plt.streamplot(X, Y, U, V, density=1.2)
plt.gca().set_aspect("equal")
plt.show()

5.7 Parametric and polar plots

Used in: Ch. 25 (parametric curves), Ch. 26 (polar coordinates), Ch. 28 (vector-valued functions).

# Parametric: a circle
t = np.linspace(0, 2*np.pi, 400)
plt.plot(np.cos(t), np.sin(t))
plt.gca().set_aspect("equal"); plt.show()

# Polar: a four-petal rose r = cos(2θ)
theta = np.linspace(0, 2*np.pi, 800)
r = np.cos(2*theta)
ax = plt.subplot(111, projection="polar")
ax.plot(theta, r)
plt.show()

6. Worked Mini-Recipes

These combine the pieces above into the complete workflows the chapters use. They embody the book's three-tier pattern (continuity §4): pose, solve by hand, verify by machine.

6.1 Verify a hand-computed derivative

Ch. 7. Confirm that $\frac{d}{dx}\big[x^2 e^x\big] = 2x e^x + x^2 e^x$.

x = sp.symbols('x')
hand = 2*x*sp.exp(x) + x**2*sp.exp(x)
machine = sp.diff(x**2 * sp.exp(x), x)
sp.simplify(hand - machine)        # → 0   ✓  (they agree)

6.2 Numerical vs. symbolic integral

Ch. 13–14. Cross-check $\int_0^1 e^{-x^2}\,dx$ three ways.

x = sp.symbols('x')
exact = sp.integrate(sp.exp(-x**2), (x, 0, 1))   # symbolic: sqrt(pi)*erf(1)/2
float(exact)                                     # → 0.7468241328124271

val, _ = integrate.quad(lambda x: np.exp(-x**2), 0, 1)   # SciPy
# val → 0.7468241328124271   (matches to ~1e-15)

riemann = np.sum(np.exp(-(np.linspace(0, 1, 1000) + 0.0005)**2)) * 0.001
# riemann ≈ 0.74682   (midpoint sum, agrees to ~4 decimals)

6.3 Solve and plot an ODE

Ch. 19. Logistic growth $\dot P = 0.8\,P\,(1 - P/100),\ P(0)=5$.

def logistic(t: float, P: np.ndarray) -> list[float]:
    return [0.8 * P[0] * (1 - P[0]/100)]

sol = integrate.solve_ivp(logistic, (0, 20), [5.0],
                          t_eval=np.linspace(0, 20, 200))
plt.plot(sol.t, sol.y[0])
plt.xlabel("t"); plt.ylabel("P(t)")
plt.show()
# Curve rises sigmoidally from 5 toward the carrying capacity 100.
# sol.y[0, -1] ≈ 99.99 (essentially saturated by t = 20)

6.4 Visualize Taylor-series convergence

Ch. 23. Watch partial sums of the Taylor series of $\sin x$ close in on the function.

x = sp.symbols('x')
xs = np.linspace(-2*np.pi, 2*np.pi, 400)
plt.plot(xs, np.sin(xs), 'k', lw=2, label="sin x")
for order in (3, 5, 7, 9):
    poly = sp.series(sp.sin(x), x, 0, order).removeO()
    f = sp.lambdify(x, poly, 'numpy')
    plt.plot(xs, f(xs), '--', label=f"order {order-1}")
plt.ylim(-2, 2); plt.legend(); plt.show()
# Higher-order curves hug sin x over a wider interval around 0.

6.5 Vector field with divergence and curl

Ch. 34, 37. Compute $\nabla\cdot\mathbf{F}$ and $\nabla\times\mathbf{F}$ symbolically, then plot.

x, y = sp.symbols('x y')
Fx, Fy = x*y, y**2 - x         # F = (xy, y² - x)

div = sp.diff(Fx, x) + sp.diff(Fy, y)            # → y + 2*y = 3*y
curl_z = sp.diff(Fy, x) - sp.diff(Fx, y)         # → -1 - x   (k-component)

# Plot the field:
xs = np.linspace(-2, 2, 15)
X, Y = np.meshgrid(xs, xs)
plt.quiver(X, Y, X*Y, Y**2 - X)
plt.gca().set_aspect("equal"); plt.show()

Here $\nabla\cdot\mathbf{F} = 3y$ measures local expansion and the scalar curl $\nabla\times\mathbf{F} = -1 - x$ measures local rotation.

7. Quick Function Index

Add to Your Modeling Portfolio — Whichever track you chose, your final notebook will use this appendix end-to-end: SymPy to derive the model's equations, SciPy solve_ivp or minimize to run it, NumPy to process the output, and Matplotlib to present it. Bookmark Recipe 6.3 (Biology/SIR), 6.5 (Physics fields), 4.6–4.7 (Economics/Data Science optimization).