Python Lab Reference

Turn a Python list (or list of lists) into a vector or matrix.

v = np.array([1.0, 2.0, 3.0])
A = np.array([[2., 1.], [1., 3.]])   # a 2×2 matrix

num evenly spaced points including both endpoints — your standard sample grid.

x = np.linspace(0, 1, 101)   # 0.00, 0.01, … 1.00

Values from start up to (but excluding) stop, spaced by step.

k = np.arange(0, 10)      # 0..9
h = np.arange(0.1, 1.0, 0.1)

num points spaced evenly in log10 — perfect for step-size / error sweeps.

hs = np.logspace(-1, -6, 20)   # 1e-1 down to 1e-6

Pre-allocate an array of zeros or ones. shape can be n or (rows, cols).

x = np.zeros(n)          # solution vector
A = np.zeros((n, n))     # coefficient matrix to fill in

Same shape & dtype as an existing array — great inside loops that build results.

errs = np.zeros_like(hs)   # one slot per step size

The n×n identity matrix.

I = np.eye(3)

Make a diagonal matrix from a vector — or extract the diagonal of a matrix.

D = np.diag([1., 2., 3.])   # vector → diagonal matrix
d = np.diag(A)              # matrix → its diagonal

Vandermonde matrix — the design matrix for polynomial fitting/interpolation.

V = np.vander(x, len(x))   # solve V c = y for polynomial coeffs

Lower / upper triangular part of a matrix (rest zeroed) — for LU-style work.

L = np.tril(A, -1) + np.eye(n)   # unit lower-triangular

Stack arrays as rows into a taller matrix.

M = np.vstack([row1, row2, row3])

Absolute value — the backbone of every error measurement.

err = np.abs(approx - exact)

e^x and natural log (ln). np.log10 for base 10.

f = np.exp(-x**2)
rate = np.log(err[1:] / err[:-1])

Square root, elementwise.

rms = np.sqrt(np.sum(r**2) / n)

Trig functions (radians) — common test integrands and ODE right-hand sides.

y = np.sin(np.pi * x)

The constants π and e. No parentheses — they are values, not functions.

area = np.pi * r**2

Round to a number of decimal places (display / comparison).

print(np.round(A, 4))

Clamp values into [lo, hi] — guards against out-of-range arguments.

theta = np.arccos(np.clip(c, -1, 1))

Largest / smallest element — max is how you report worst-case error.

max_err = np.max(np.abs(u - u_exact))

INDEX of the smallest / largest element (not the value).

i = np.argmin(errors)   # which step size was best

Sum of elements — quadrature rules, norms, residuals.

integral = h * np.sum(f)   # basic Riemann/trapezoid core

Consecutive differences — finite-difference derivatives and spacing checks.

slopes = np.diff(y) / np.diff(x)

Are arrays equal within tolerance? Use this, never == on floats.

assert np.allclose(A @ x, b)   # did the solver actually solve it?

Sort values / find the insertion index in a sorted array (piecewise lookup).

i = np.searchsorted(nodes, x)   # which interval x falls in

Fast piecewise-linear interpolation at points x from samples (xp, fp).

y = np.interp(0.37, xp, fp)

Solve the linear system A x = b directly. The right way to “apply A⁻¹”.

x = np.linalg.solve(A, b)

The matrix inverse A⁻¹.

Ainv = np.linalg.inv(A)

Vector or matrix norm — measure residuals and error magnitudes.

residual = np.linalg.norm(A @ x - b)
rel_err  = np.linalg.norm(x - x_true) / np.linalg.norm(x_true)

Condition number κ(A) — how much A amplifies input error. The course’s key diagnostic.

print(f"cond = {np.linalg.cond(A):.2e}")

Eigenvalues of A (values only). Use eig() if you also need eigenvectors.

lam = np.linalg.eigvals(A)
spectral_radius = np.max(np.abs(lam))

Least-squares solution to an overdetermined system (best fit).

coef, residuals, rank, sv = np.linalg.lstsq(A, b, rcond=None)

Cholesky factor L (A = L Lᵀ) for symmetric positive-definite matrices.

L = np.linalg.cholesky(A)

QR factorization A = Q R (orthogonal × upper-triangular).

Q, R = np.linalg.qr(A)

Machine epsilon — the gap between 1.0 and the next float (≈ 2.22e-16 for float64).

eps = np.finfo(float).eps
print(f"machine epsilon = {eps:.2e}")

Positive infinity — sentinel for “no bound yet” or a diverged value.

best = np.inf   # then keep the minimum found

Fix the random seed so results are reproducible (and match the grader).

np.random.seed(42)

Uniform random values in [0, 1) with the given shape.

A = np.random.rand(n, n)

Standard-normal (mean 0, std 1) random values — noise for perturbation tests.

noise = 1e-6 * np.random.randn(n)

Create a figure and one (or a grid of) axes. The standard opening line.

fig, ax = plt.subplots(figsize=(7, 4))
fig, (ax1, ax2) = plt.subplots(1, 2)   # two side-by-side

Auto-adjust spacing so labels/titles don’t overlap or clip.

plt.tight_layout()

Write the figure to a file — this is how the lab captures your plot.

plt.savefig('convergence.png', dpi=150, bbox_inches='tight')

Render the figure. In the lab runner, pair with savefig so output is captured.

plt.show()

Line/marker plot. Third arg is a style string; label feeds the legend.

ax.plot(x, y, 'b.-', label='approx')

Linear x, log y — shows exponential decay of error as a straight line.

ax.semilogy(iterations, residuals)

Log-log plot — the slope IS the order of accuracy p (err ∝ hᵖ).

ax.loglog(hs, errors, "o-")

Label the axes and title the plot.

ax.set_xlabel('step size h'); ax.set_ylabel('error'); ax.set_title('Convergence')

Show the legend built from each plot’s label=.

ax.legend()

Add grid lines — helpful for reading slopes off log plots.

ax.grid(True, which='both', alpha=0.3)

Draw a horizontal reference line (e.g. a tolerance or exact value).

ax.axhline(eps, color='r', ls='--', label='machine eps')

Fix the y-axis range so tiny errors don’t collapse to nothing.

ax.set_ylim(1e-16, 1)

Build a smooth natural cubic spline through the data points, then call it like a function.

cs = CubicSpline(x, y)
y_fine = cs(x_fine)          # evaluate anywhere
d1 = cs(x_fine, 1)           # first derivative

Robustly find a root of f in the bracket [a, b] where f changes sign.

root = brentq(lambda x: x**2 - 2, 0, 2)   # → √2

Adaptive high-accuracy definite integral — the reference “exact” value for quadrature labs.

val, err = quad(lambda x: np.exp(-x**2), 0, 1)

Integrate an ODE system dy/dt = fun(t, y) — the reference for your Euler/RK code.

def f(t, y):
    return -y                      # y' = -y
sol = solve_ivp(f, [0, 5], [1.0], t_eval=np.linspace(0, 5, 100))
y = sol.y[0]                       # solution over time

Factor A once, then solve for many right-hand sides cheaply.

lu, piv = lu_factor(A)
x1 = lu_solve((lu, piv), b1)
x2 = lu_solve((lu, piv), b2)   # reuse the factorization

The n×n Hilbert matrix — the classic ill-conditioned test case.

H = hilbert(6)
print(np.linalg.cond(H))   # enormous

Show a set number of decimal places — for tables of approximations.

print(f"x = {x:.6f}")

Scientific notation — the natural format for errors and condition numbers.

print(f"error = {err:.3e},  cond = {kappa:.2e}")

Pad to a width for aligned columns in convergence tables.

print(f"{n:>4}{err:>14.3e}{ratio:>10.2f}")

Define a function. 4-space indented body; remember to return.

def newton(f, df, x0, tol=1e-10):
    x = x0
    while abs(f(x)) > tol:
        x = x - f(x) / df(x)
    return x

A one-line anonymous function — handy for passing f to a solver.

root = brentq(lambda x: x**3 - x - 2, 1, 2)

Iterate a fixed number of times, or until a tolerance is met.

for k in range(max_iter):
    x_new = g(x)
    if abs(x_new - x) < tol:
        break
    x = x_new

Scalar math helpers — e.g. estimating how many iterations bisection needs.

n_steps = math.ceil(math.log((b - a) / tol, 2))

Gauss–Legendre quadrature nodes & weights on [-1, 1].

x, w = np.polynomial.legendre.leggauss(5)
# map [-1,1] → [a,b] then integral ≈ np.sum(w * f(xm)) * (b-a)/2

Absolute error is |approx − exact|; relative divides by |exact| (scale-free).

abs_err = np.abs(approx - exact)
rel_err = abs_err / np.abs(exact)

If error ∝ hᵖ, then log(err) vs log(h) is a line of slope p.

p = np.polyfit(np.log(hs), np.log(errs), 1)[0]
print(f"observed order ≈ {p:.2f}")

Halving h and watching the error ratio reveals the order (ratio ≈ 2ᵖ).

ratio = errs[:-1] / errs[1:]   # ~2 ⇒ 1st order, ~4 ⇒ 2nd order

Always verify the residual after solving Ax=b.

assert np.linalg.norm(A @ x - b) < 1e-8

Rounding makes exact equality almost always False.

# WRONG: if x == 0.3
if np.isclose(x, 0.3): ...

Subtracting two nearly equal numbers destroys significant digits.

# unstable:  (1 - np.cos(x)) / x**2   for small x
# stable:    0.5 * (np.sinc(x/np.pi))**2  or a Taylor form

Forming the inverse is slower and less accurate than solving directly.

# avoid: x = np.linalg.inv(A) @ b
x = np.linalg.solve(A, b)

A big condition number means the answer is untrustworthy no matter the algorithm.

if np.linalg.cond(A) > 1e12:
    print("warning: nearly singular")

np.arange excludes the stop and can drop the last point due to rounding.

# safer for a known count:
x = np.linspace(a, b, num)

Python loops over arrays are slow and error-prone — operate on whole arrays.

# slow: for i in range(n): y[i] = f(x[i])
y = f(x)   # if f uses np.* it already vectorizes