Math for Data Science · MDS-8
Turn the math you know into algorithms a computer can actually run — fast, stable, and accurate enough to trust.
Real data science runs on finite-precision arithmetic, so the clean formulas from calculus and linear algebra rarely transfer to code unchanged. This course teaches the algorithms that scientific computing actually uses — how they introduce and amplify error, when they are stable, and how to make them fast. Every idea is implemented from scratch and then compared against the production routines in NumPy/SciPy or Julia.
Course Outline
Every lesson is self-contained: a full explanation, worked step-by-step examples (in Python and R/Julia), common mistakes, a try-it-yourself, and an interactive quiz. Jump to a unit:
Computers store numbers in a fixed number of bits, so the continuous real line must be approximated by a finite, unevenly spaced set of values.
The real numbers are uncountably infinite and dense, but any machine word has finitely many bits, so it can represent only finitely many values. Floating point picks a clever finite subset: numbers of the form m * 2^e where the mantissa m has fixed precision and the exponent e ranges over a wide span. This gives huge dynamic range but uneven spacing — the gap between consecutive representable numbers grows with magnitude. Most decimals (even 0.1) are not exactly representable, so storing them introduces representation error. Understanding that nearly every stored value is already slightly wrong is the foundation of numerical analysis: we never trust the last few digits, and we track how errors propagate through computation.
Worked Example 1
Problem. Show that 0.1 is not exactly representable in float64.
Answer. 0.1 is stored as 0.10000000000000000555..., off by ~5.55e-18; it is not exact.
Worked Example 2
Problem. Demonstrate uneven spacing: find the gap between consecutive floats near 1.0 versus near 1e9.
Answer. Spacing is ~2.22e-16 near 1.0 but ~1.19e-7 near 1e9 — representable numbers thin out as magnitude grows.
Problem. Does 0.1 + 0.2 == 0.3 evaluate to True in float64? Explain and give a robust check.
Solution. `0.1 + 0.2` yields 0.30000000000000004 because each operand and the sum carry representation error, so `0.1 + 0.2 == 0.3` is False. Use `np.isclose(0.1 + 0.2, 0.3)`, which compares within a relative+absolute tolerance and returns True. Answer: False with exact ==, True with np.isclose.
IEEE 754 fixes exactly how the bits of a float encode sign, exponent, and mantissa, which determines precision and the value of machine epsilon.
A float64 uses 64 bits: 1 sign bit, 11 exponent bits (biased by 1023), and 52 mantissa (fraction) bits. The value is (-1)^s * 1.f * 2^(e-1023), where the leading 1 is implicit for normal numbers, giving 53 bits of precision. Machine epsilon is the gap between 1.0 and the next representable float: 2^-52 ≈ 2.22e-16 for float64. It bounds the relative round-off of a single rounded operation. Special bit patterns encode +/-inf (all-ones exponent, zero mantissa) and nan (all-ones exponent, nonzero mantissa). Knowing the layout lets you predict precision (float32 has 23 mantissa bits, eps ≈ 1.19e-7) and explain why doubling precision roughly halves the digit count of error.
Worked Example 1
Problem. Confirm machine epsilon for float64 and float32 from NumPy.
Answer. float64 eps = 2**-52 ≈ 2.22e-16; float32 eps = 2**-23 ≈ 1.19e-7.
Worked Example 2
Problem. Decode the float64 bits of 1.0 and of the next float after 1.0.
Answer. 1.0 = 1.0 * 2^0; the next float is 1.0 + 2**-52, so the step equals machine epsilon.
Problem. How many decimal significant digits can float64 reliably represent, and how does that follow from 52 mantissa bits?
Solution. Precision is 53 bits (52 stored + 1 implicit). Decimal digits ≈ 53 * log10(2) ≈ 53 * 0.30103 ≈ 15.95. So float64 carries about 15-16 reliable significant decimal digits. Answer: ~15-16 digits.
Every operation rounds to the nearest float, and subtracting two nearly equal values can wipe out the leading correct digits.
When a real result lands between two representable floats, IEEE 754 rounds to nearest (ties to even), introducing a relative error bounded by eps/2. This per-operation round-off is usually harmless, but subtraction is dangerous: if x and y agree to k leading digits, x - y cancels those digits and leaves only the noisy trailing bits, magnifying relative error by roughly the ratio of magnitudes — catastrophic cancellation. The cure is algebraic reformulation that avoids the subtraction, such as rationalizing, using log1p/expm1, or the stable quadratic formula. The data was never lost in the inputs; cancellation merely exposes that the inputs themselves carried round-off in their low bits.
Worked Example 1
Problem. Evaluate (1 - cos(x)) / x**2 at x = 1e-8 directly and stably; the true limit is 0.5.
Answer. Direct gives 0.0 (totally wrong); the stable form gives 0.5, the correct value.
Worked Example 2
Problem. Compute exp(x) - 1 at x = 1e-15 accurately.
Answer. Use np.expm1(x) → ~1.0e-15; the naive exp(x)-1 loses nearly all significant digits.
Problem. Compute log(1 + x) for x = 1e-16 accurately and explain the fix.
Solution. `np.log(1 + 1e-16)` gives 0.0 because 1 + 1e-16 rounds to exactly 1.0 (1e-16 < eps). Use `np.log1p(1e-16)`, which gives 1e-16, computing log(1+x) without forming 1+x. Answer: 0.0 naively, 1e-16 with log1p.
Absolute error is the raw gap from the truth; relative error scales that gap by the magnitude and ties directly to significant digits.
If x_hat approximates true value x, the absolute error is |x_hat - x| and the relative error is |x_hat - x| / |x| (for x != 0). Absolute error answers 'how far off in raw units,' while relative error answers 'what fraction of the value is wrong,' which is scale-independent and comparable across magnitudes. A relative error near 10^-k corresponds to roughly k correct significant decimal digits. Because float64 round-off is bounded in relative terms (eps ≈ 2.22e-16), relative error is the natural currency of numerical analysis. Reporting only absolute error hides whether 0.001 is excellent (for values near 1e6) or terrible (for values near 0.001).
Worked Example 1
Problem. Approximate pi by 22/7. Find absolute and relative error and the number of correct significant digits.
Answer. Absolute ≈ 1.264e-3, relative ≈ 4.02e-4, about 3 correct significant digits.
Worked Example 2
Problem. Two estimates: A = 1000.5 for true 1000.0, B = 0.0005 for true 0.0. Which is relatively better?
Answer. A has relative error 5e-4 (~3 good digits); B's relative error is undefined near zero, so use absolute error there.
Problem. An estimate gives 2.71850 for e = 2.718281828. Find relative error and approximate correct significant digits.
Solution. Absolute error = |2.71850 - 2.718281828| = 0.000218172. Relative error = 0.000218172 / 2.718281828 ≈ 8.03e-5. Correct digits ≈ -log10(8.03e-5) ≈ 4.1, so about 4 significant digits. Answer: rel err ≈ 8.0e-5, ~4 digits.
Forward error measures how wrong the answer is; backward error measures how wrong an input would have to be to make that answer exact; conditioning links the two.
Given a problem y = f(x) and a computed y_hat, the forward error is |y_hat - y|, how far the output is from the truth. The backward error is the smallest |delta x| such that f(x + delta x) = y_hat exactly — how much we'd have to perturb the input to justify our answer. An algorithm is backward stable if it always has tiny backward error. The condition number of the problem amplifies backward error into forward error: forward error ≲ condition number * backward error. This separation is powerful: a backward-stable algorithm on an ill-conditioned problem can still produce a large forward error, and that is the problem's fault, not the algorithm's.
Worked Example 1
Problem. For f(x) = sqrt(x) at x = 2, a routine returns y_hat = 1.41421356. Find forward and backward error.
Answer. Forward error ≈ 2.37e-9; backward error ≈ 6.7e-9 (the input would need to be ~2.0000000067).
Worked Example 2
Problem. Estimate the condition number of f(x) = sqrt(x) at x = 2 and relate the two errors above.
Answer. Condition number is 1/2; it correctly predicts the forward error is about half the (relative) backward error.
Problem. A solver returns y_hat = 3.0001 for the problem y = f(x) whose true answer is 3.0. The problem's condition number is about 100. Roughly what backward error does this imply?
Solution. Relative forward error ≈ |3.0001 - 3.0| / 3.0 ≈ 3.33e-5. Since forward ≈ condition * backward, relative backward error ≈ forward/condition ≈ 3.33e-5 / 100 ≈ 3.33e-7. A tiny backward error on a moderately conditioned problem — the algorithm is fine. Answer: backward error ≈ 3.3e-7 (relative).
Summing many floats accumulates round-off, and compensated (Kahan) summation recovers most of the lost low-order bits.
Floating-point addition is not associative: each partial sum is rounded, so adding n numbers can accumulate an error that grows like n * eps in the worst case, and small terms can vanish entirely against a large running total. Naive left-to-right summation is therefore unreliable for long or ill-scaled sequences. Kahan (compensated) summation keeps a running correction term that captures the low-order bits lost in each addition and feeds them back on the next step, reducing the error bound to roughly 2 * eps independent of n. Sorting by magnitude or pairwise (tree) summation also helps; NumPy's np.sum uses pairwise summation by default, which is why it is more accurate than a Python loop.
Worked Example 1
Problem. Sum 1.0 plus ten million copies of 1e-8. The exact answer is 1.1.
Answer. Naive accumulation drifts from 1.1; np.sum (pairwise) and math.fsum recover 1.1.
Worked Example 2
Problem. Implement Kahan summation and compare to a naive sum on a hard case.
Answer. Kahan summation tracks a correction term and recovers 1.1 where the naive loop loses accuracy.
Problem. Why does np.sum often beat a hand-written Python for-loop for adding a large array, even with identical inputs?
Solution. A left-to-right loop rounds after every one of n additions, so error can grow like n*eps and small terms vanish into a large running sum. np.sum uses pairwise (tree) summation, which adds in balanced groups so error grows like log(n)*eps instead of n*eps, and runs in compiled code. Answer: pairwise summation gives a much smaller, log(n)-scaled error bound than naive accumulation.
Implement the textbook quadratic formula in float64 and find coefficients (such as a very small a or b much larger than 4ac) where it returns one badly inaccurate root. Use numpy.float64 and compare against a high-precision reference. Then reimplement the formula using the numerically stable variant that avoids cancellation and confirm the error disappears. Quantify both the absolute and relative error before and after.
Deliverable · A NumPy script plus a short write-up showing a failing input, the error magnitude, the stable rewrite, and a sentence explaining the cancellation that caused the failure.
1. What does machine epsilon for float64 represent?
Answer B. Machine epsilon is the relative spacing near 1.0; for float64 it is 2**-52, roughly 2.22e-16. The smallest positive value is a different (much smaller) subnormal quantity.
2. Catastrophic cancellation occurs when you:
Answer C. Subtracting nearly equal values cancels the leading correct digits and exposes the rounding error in the low-order bits, inflating relative error.
3. An algorithm has small backward error when:
Answer B. Backward error measures how much the input must change so the computed result is exact; a backward-stable algorithm keeps that perturbation tiny.
4. Relative error is preferred over absolute error mainly because it:
Answer B. Relative error divides by the true magnitude, so an error of 1 means something very different for a value near 1 than near 1e9; this makes it comparable across scales.
5. Summing many positive floats of widely different magnitudes naively tends to:
Answer B. Floating-point addition is not associative; small terms can be smaller than the rounding step of the running sum and vanish. Compensated (Kahan) summation reduces this.
I can explain how IEEE 754 represents numbers and where machine epsilon comes from.
I can identify and reformulate computations that suffer from catastrophic cancellation.
I can distinguish forward error, backward error, and conditioning when reasoning about an algorithm's accuracy.
Most equations have no closed-form root, so we replace exact algebra with an iterative process that drives a residual toward zero within a tolerance.
Solving f(x) = 0 numerically means producing a sequence of approximations x_0, x_1, ... whose residuals f(x_n) shrink toward machine precision, rather than an exact symbolic root that usually does not exist. We almost never reach the true root: we stop when the residual or the step size falls below a tolerance. Two ideas frame everything: a bracket is an interval [a, b] with f(a) and f(b) of opposite sign, which (by continuity) traps a root and gives a reliability guarantee; an open method like Newton starts from a guess and may converge fast or diverge. Numerically, beware that |f(x)| small does not always mean x is close to a root — for flat (ill-conditioned) functions a tiny residual can hide a large error.
Worked Example 1
Problem. Reformulate 'find where cos(x) equals x' as a root problem and verify a candidate.
Answer. f(x)=cos(x)-x has a root near 0.739085; sign change on (0,1) confirms a bracketed root exists.
Worked Example 2
Problem. Show that a small residual can mislead for a flat function f(x) = (x - 1)**3.
Answer. Residual 1e-6 but error 1e-2: at a flat/multiple root, small |f| does not imply a close x.
Problem. Turn 'x e^x = 1' into f(x) = 0, find a bracket, and confirm the root lies inside.
Solution. Let f(x) = x*np.exp(x) - 1. Evaluate f(0) = -1 and f(1) = e - 1 ≈ 1.718, opposite signs, so a root is bracketed in (0, 1). Checking x = 0.567143: `f=lambda x: x*np.exp(x)-1; print(f(0.567143))` gives about 1.8e-7. Answer: root ≈ 0.567143 (the Omega constant), bracketed in (0, 1).
Bisection repeatedly halves a sign-change interval, guaranteed by the intermediate value theorem to enclose a root.
If f is continuous on [a, b] and f(a) and f(b) have opposite signs, the intermediate value theorem guarantees a root inside. Bisection exploits this: evaluate the midpoint m = (a+b)/2, keep whichever half still shows a sign change, and repeat. The interval width halves each step, so after n steps the uncertainty is (b-a)/2^n — linear convergence with rate 1/2, gaining about one bit (0.3 decimal digits) per iteration. It is slow but bulletproof: it cannot diverge and needs no derivative, only continuity and a valid bracket. Numerically, compute the midpoint as a + (b-a)/2 rather than (a+b)/2 to avoid overflow, and test the sign with the product's sign (or np.sign) rather than multiplying, which can underflow.
Worked Example 1
Problem. Do two bisection steps on f(x) = x**2 - 2 over [1, 2] to bracket sqrt(2).
Answer. After two steps the bracket is [1.25, 1.5]; width halved twice from 1 to 0.25.
Worked Example 2
Problem. How many bisection steps are needed to locate a root in [0, 1] to within 1e-6?
Answer. About 20 halvings; n = 19 iterations guarantee error <= 1e-6 on a unit interval.
Worked Example 3
Problem. Implement bisection and find the root of f(x)=x**3 - x - 2 on [1, 2].
Answer. Root ≈ 1.5213797068, residual ~2e-10; bisection converged reliably.
Problem. Use bisection logic to bracket the root of f(x) = cos(x) - x in [0, 1] for three steps.
Solution. f(0)=1>0, f(1)=cos(1)-1≈-0.4597<0: root in [0,1]. m=0.5, f(0.5)=cos(0.5)-0.5≈0.3776>0, keep [0.5,1]. m=0.75, f(0.75)≈cos(0.75)-0.75≈-0.0183<0, keep [0.5,0.75]. m=0.625, f(0.625)≈0.1860>0, keep [0.625,0.75]. Answer: after three steps the root (~0.739085) is bracketed in [0.625, 0.75].
Rewriting f(x)=0 as x = g(x) and iterating x_{n+1} = g(x_n) converges when g is a local contraction, |g'| < 1.
A fixed point of g satisfies x* = g(x*). To solve f(x) = 0 we algebraically rearrange it into x = g(x) and iterate x_{n+1} = g(x_n) from a starting guess. Convergence is governed by the derivative at the fixed point: if |g'(x*)| < 1 the map is locally contracting and the error shrinks by roughly that factor each step (linear convergence with rate |g'(x*)|); if |g'(x*)| > 1 the iteration diverges no matter how close you start. The same equation can be rearranged many ways, and only some give a convergent g — choosing a good rearrangement (small |g'|) is the whole art. When g'(x*) happens to be 0, convergence becomes faster than linear, which is exactly what Newton's method engineers.
Worked Example 1
Problem. Solve x = cos(x) by fixed-point iteration from x_0 = 1; check the contraction condition.
Answer. Converges to x* ≈ 0.739085 because |g'(x*)| ≈ 0.67 < 1; convergence is linear.
Worked Example 2
Problem. For x**2 = 2, compare the divergent rearrangement g(x)=2/x with the convergent g(x)=(x + 2/x)/2.
Answer. g=2/x stalls (|g'|=1); the averaged form has g'=0 at the root and converges quadratically to sqrt(2).
Problem. To solve x = e^{-x}, will iterating x_{n+1} = e^{-x_n} from x_0 = 0.5 converge? Find the fixed point.
Solution. g(x) = e^{-x}, g'(x) = -e^{-x}. The fixed point is near x* ≈ 0.567; there |g'| = e^{-0.567} ≈ 0.567 < 1, so it converges (linearly). Iterating: `x=0.5\nfor _ in range(30): x=np.exp(-x)` gives x ≈ 0.5671432904. Answer: converges to x* ≈ 0.567143 since |g'(x*)| ≈ 0.567 < 1.
Newton's method linearizes f at the current guess and jumps to the line's root, doubling the number of correct digits each step near a simple root.
Newton's method approximates f near x_n by its tangent line and takes the next iterate where that line crosses zero: x_{n+1} = x_n - f(x_n)/f'(x_n). For a simple root (f(x*)=0, f'(x*) != 0) the error obeys e_{n+1} ≈ (f''/2f') e_n^2, so it converges quadratically — the number of correct digits roughly doubles each step. The cost is one f and one f' evaluation per step, plus fragility: a near-zero derivative sends the step to infinity, a bad initial guess can diverge or cycle, and at a multiple root convergence drops to linear. There is no bracket, so Newton offers speed without a reliability guarantee. In practice it is often safeguarded by falling back to bisection when a step leaves the bracket.
Worked Example 1
Problem. Apply two Newton steps to f(x) = x**2 - 2 from x_0 = 1 to approximate sqrt(2).
Answer. x_1=1.5, x_2≈1.41667, x_3≈1.4142157 — digits roughly double each step (quadratic).
Worked Example 2
Problem. Implement Newton's method for f(x)=x**3 - x - 2 from x_0=1.5 and observe quadratic convergence.
Answer. Converges to 1.5213797068045676 in ~4 steps; error squares each iteration.
Worked Example 3
Problem. Show Newton failing: f(x) = x**3 - 2x + 2 starting at x_0 = 0.
Answer. Newton enters a 2-cycle (0 → 1 → 0 → ...) and never converges; a poor starting point defeats it.
Problem. Use Newton's method to find a root of f(x) = cos(x) - x from x_0 = 0.5; do two steps.
Solution. f'(x) = -sin(x) - 1. x_1 = 0.5 - (cos(0.5)-0.5)/(-sin(0.5)-1) = 0.5 - (0.37758)/(-1.47943) ≈ 0.5 + 0.25522 = 0.75522. x_2 = 0.75522 - (cos(0.75522)-0.75522)/(-sin(0.75522)-1) ≈ 0.75522 - (-0.01865)/(-1.68507) ≈ 0.73914. Answer: x_2 ≈ 0.739085, converging quadratically to the true root.
The secant method replaces Newton's analytic derivative with a finite difference of the two latest iterates, needing no derivative.
When the derivative is unavailable or expensive, the secant method approximates f'(x_n) by the slope through the last two points: x_{n+1} = x_n - f(x_n) * (x_n - x_{n-1}) / (f(x_n) - f(x_{n-1})). It is Newton with a numerical derivative, so it needs two starting values and only one new function evaluation per step. Its convergence order is the golden ratio, about 1.618 — superlinear, slower than Newton's 2 but faster than bisection's 1, and often more efficient per function evaluation since it never computes f'. Like Newton it can diverge without a bracket, and it breaks down if two consecutive iterates give nearly equal f values (the denominator cancels). It underlies SciPy's derivative-free newton when no fprime is supplied.
Worked Example 1
Problem. Take two secant steps on f(x) = x**2 - 2 with x_0 = 1, x_1 = 2.
Answer. x_2 = 1.33333, x_3 ≈ 1.4 — superlinear approach to sqrt(2) using no derivative.
Worked Example 2
Problem. Implement the secant method for f(x)=x**3 - x - 2 with x_0=1, x_1=2.
Answer. Converges to 1.5213797068 in ~6 steps; order ≈ 1.618 (golden ratio), no derivative required.
Problem. Do one secant step for f(x) = e^x - 3 with x_0 = 1, x_1 = 1.5.
Solution. f(1) = e - 3 ≈ -0.28172, f(1.5) = e^1.5 - 3 ≈ 1.48169. x_2 = 1.5 - 1.48169*(1.5-1)/(1.48169-(-0.28172)) = 1.5 - 1.48169*0.5/1.76341 = 1.5 - 0.42010 = 1.07990. The true root is ln(3) ≈ 1.09861. Answer: x_2 ≈ 1.0799, already close to ln(3) ≈ 1.0986 after one step.
Convergence order quantifies how fast error shrinks, and sensible stopping rules combine step size, residual, and an iteration cap.
An iteration converges with order p if |e_{n+1}| ≈ C |e_n|^p near the root: p = 1 is linear (constant digit gain), p ≈ 1.618 is the secant's superlinear rate, p = 2 is Newton's quadratic (digits double). To decide when to stop without knowing the true root, you combine criteria: a step tolerance |x_{n+1} - x_n| < xtol, a residual tolerance |f(x_n)| < ftol, and always a maximum iteration count to bound runaway cases. Use a mixed absolute-plus-relative tolerance (xtol + rtol*|x|) so the test scales with the root's magnitude. Demanding a tolerance below machine epsilon times the value is futile — the iteration stalls on round-off. Estimating p empirically from log-error ratios is a good sanity check on an implementation.
Worked Example 1
Problem. From the Newton errors 1e-2, 1e-4, 1e-8, estimate the convergence order p.
Answer. p = 2: the errors square each step, confirming quadratic (Newton) convergence.
Worked Example 2
Problem. Design a robust stopping rule for a Newton loop on a root near x ≈ 1000.
Answer. Combine xtol + rtol*|x| with a hard iteration cap (e.g. 50); the relative term scales the test to the root's magnitude.
Problem. A method gives errors 0.2, 0.04, 0.008, 0.0016. Is it linear, and what is the rate?
Solution. Each error is the previous times 0.2 (0.04/0.2 = 0.2, 0.008/0.04 = 0.2, 0.0016/0.008 = 0.2). A constant ratio means linear convergence (p = 1) with rate C = 0.2. Order estimate: log(0.008/0.0016... ) approach gives ratio constant, so p=1. Answer: linear convergence, rate 0.2 (each step gains ~log10(1/0.2) ≈ 0.7 digits).
SciPy's brentq blends bisection's safety with fast interpolation, and root_scalar offers a unified interface to many methods.
Brent's method (scipy.optimize.brentq) is the workhorse for a known bracket: it combines guaranteed bisection convergence with faster inverse-quadratic interpolation and secant steps, automatically falling back to bisection whenever the fast step misbehaves. Given f(a), f(b) of opposite sign it cannot fail and typically converges superlinearly — the best of both worlds, with no derivative needed. scipy.optimize.root_scalar is a dispatcher: pass method='brentq'/'bisect' with a bracket, or method='newton'/'secant' with a starting point (and fprime for Newton), and it returns a result object with the root, iteration count, and convergence flag. Rule of thumb: if you can bracket the root, use brentq for robustness; reach for Newton only when the derivative is cheap and a good initial guess is available.
Worked Example 1
Problem. Solve x**3 - x - 2 = 0 with scipy.optimize.brentq on the bracket [1, 2].
Answer. brentq returns 1.5213797068045751 reliably; it stays bracketed and converges superlinearly.
Worked Example 2
Problem. Use scipy.optimize.root_scalar with Newton and with brentq, and inspect the result object.
Answer. Newton converged in 4 iterations; brentq took more steps but is bracket-safe — both return 1.52138 with converged=True.
Problem. Use scipy.optimize.brentq to find where cos(x) = x on a suitable bracket.
Solution. Define f(x) = np.cos(x) - x. A bracket is [0, 1] since f(0)=1>0 and f(1)=cos(1)-1<0. `from scipy.optimize import brentq; import numpy as np; r = brentq(lambda x: np.cos(x)-x, 0, 1); print(r)` gives 0.7390851332151607. Answer: brentq returns x ≈ 0.7390851332, the unique solution of cos(x)=x.
Pick a nonlinear function with a known root (for example f(x) = x**3 - x - 2) and implement bisection, Newton's method, and the secant method yourself. Track the error at each iteration and verify the expected convergence orders by plotting log error against iteration. Then solve the same problem with scipy.optimize.brentq and scipy.optimize.newton and compare iteration counts and final accuracy.
Deliverable · A script and short report with a convergence plot, a table of iterations-to-tolerance per method, and one sentence on why Newton converged fastest yet was least robust.
1. Bisection is guaranteed to converge as long as:
Answer B. Bisection only needs a sign change (a bracket) and continuity; it then halves the interval reliably regardless of smoothness.
2. Near a simple root, Newton's method typically converges:
Answer B. For a simple root with a nonzero derivative, the error roughly squares each step, giving quadratic (order 2) convergence.
3. Fixed-point iteration x_{n+1} = g(x_n) converges to a fixed point when, near it:
Answer A. The map is a local contraction when |g'| < 1, which guarantees convergence; |g'| > 1 makes the iteration diverge.
4. The main advantage of the secant method over Newton's method is that it:
Answer B. The secant method approximates f' with a finite difference of past iterates, so no analytic derivative is needed; its order (~1.618) is between linear and quadratic.
5. Why is scipy.optimize.brentq often preferred for a known bracket?
Answer B. Brent's method blends bisection's guaranteed convergence with faster secant/inverse-quadratic steps, staying robust while accelerating.
I can implement bisection, Newton's method, and the secant method from scratch and explain their convergence rates.
I can choose sensible tolerances and stopping criteria and recognize when a method diverges.
I can select an appropriate SciPy root finder for a given problem and justify the choice.
We solve linear systems by factoring and substitution, never by forming A^-1, which is slower, less accurate, and rarely needed.
The clean math says x = A^-1 b, but numerically computing the inverse is a mistake: it costs about three times the work of a direct solve, amplifies round-off, and discards structure (sparsity, symmetry). Instead we factor A once (e.g. into triangular factors) and solve by substitution, which is both cheaper and more accurate. NumPy's np.linalg.solve and SciPy's lu_solve do exactly this. Forming the inverse is justified only when you genuinely need the matrix itself (rare). A useful diagnostic is the residual r = b - A x_hat: it is cheap to compute, but for ill-conditioned A a small residual does not guarantee a small error in x. Always prefer solve(A, b) over inv(A) @ b.
Worked Example 1
Problem. Solve the 2x2 system [[3,2],[1,2]] x = [5,5] correctly in NumPy.
Answer. x = [0, 2.5]; solve(A, b) factors A and substitutes, no inverse formed.
Worked Example 2
Problem. Show that np.linalg.solve beats inv(A) @ b in accuracy and cost on a random system.
Answer. solve gives a smaller residual and runs faster; inv(A)@b does strictly more work for a worse answer.
Problem. Solve [[2,1,1],[1,3,2],[1,0,0]] x = [4,5,6] with NumPy and verify.
Solution. `import numpy as np; A=np.array([[2.,1.,1.],[1.,3.,2.],[1.,0.,0.]]); b=np.array([4.,5.,6.]); x=np.linalg.solve(A,b); print(x)` gives `[ 6. 15. -23.]`. Verify `A @ x` returns `[4. 5. 6.]`. Answer: x = [6, 15, -23], obtained by direct solve without inverting A.
Gaussian elimination systematically zeros out subdiagonal entries; recording the multipliers yields the LU factorization A = LU.
Gaussian elimination reduces A to upper-triangular U by subtracting multiples of each pivot row from the rows below. If you store the multipliers used (the factor by which each row was scaled) in a unit lower-triangular matrix L, you obtain A = LU. This costs about (2/3)n^3 flops — the dominant expense. Once factored, solving Ax = b becomes two cheap O(n^2) triangular sweeps: forward substitution solves Ly = b, then back substitution solves Ux = y. The big win is reuse: factor A once, then solve for many right-hand sides b cheaply. Numerically, LU without pivoting can fail if a pivot is zero or tiny (dividing by it magnifies error), which is why partial pivoting (the next lesson) is essential in practice.
Worked Example 1
Problem. Hand-compute the LU factorization of A = [[2,1],[6,4]] (no pivoting).
Answer. L = [[1,0],[3,1]], U = [[2,1],[0,1]]; their product reconstructs A.
Worked Example 2
Problem. Use the LU factors to solve [[2,1],[6,4]] x = [3, 10].
Answer. x = [1, 1], found by forward then back substitution on the LU factors.
Problem. Factor A = [[4,3],[8,7]] into LU by hand and verify with NumPy.
Solution. Multiplier m = 8/4 = 2. Row2 - 2*Row1 = [8,7]-[8,6] = [0,1]. So U = [[4,3],[0,1]], L = [[1,0],[2,1]]. Check: L@U = [[4,3],[8,7]] = A. In SciPy, `from scipy.linalg import lu; P,L,U = lu(A)` gives the pivoted version (here it pivots since 8>4). Answer: L=[[1,0],[2,1]], U=[[4,3],[0,1]] for the unpivoted factorization.
Partial pivoting swaps rows to put the largest-magnitude pivot on the diagonal, preventing small pivots from amplifying round-off.
Plain Gaussian elimination divides by the pivot; a tiny pivot produces huge multipliers that blow up rounding error, and a zero pivot halts the algorithm entirely. Partial pivoting fixes this: before eliminating each column, swap in the row with the largest absolute entry in that column to serve as the pivot. This keeps all multipliers <= 1 in magnitude, bounding error growth and making the factorization backward stable in practice. The result is recorded as PA = LU, where P is a permutation matrix tracking the row swaps. SciPy's scipy.linalg.lu_factor returns the LU factors and a pivot index array; lu_solve then applies the same permutation when solving. Partial pivoting is standard in every production solver — you almost never want LU without it.
Worked Example 1
Problem. Show how a tiny pivot wrecks accuracy: solve [[1e-20, 1],[1, 1]] x = [1, 2] with and without pivoting.
Answer. Pivoting gives the correct [1, 1]; the naive no-pivot version loses the data and returns ~[0, 1].
Worked Example 2
Problem. Use scipy.linalg.lu_factor / lu_solve (which pivot) to solve a system and inspect the pivots.
Answer. lu_factor pivots automatically (piv records swaps); lu_solve applies them to return an accurate x.
Problem. For A = [[0, 2],[1, 1]], why is partial pivoting required, and what swap occurs?
Solution. The (1,1) pivot is 0, so elimination would divide by zero. Partial pivoting scans column 1 for the largest magnitude entry, finds 1 in row 2, and swaps rows: PA = [[1,1],[0,2]], which is already upper triangular. P swaps rows 1 and 2. Answer: pivoting swaps the two rows to avoid the zero pivot, giving PA = LU with U = [[1,1],[0,2]].
For symmetric positive-definite matrices, Cholesky factors A = L L^T at half the cost of LU and with guaranteed stability.
When A is symmetric positive-definite (SPD) — symmetric with all positive eigenvalues, common for covariance matrices and normal-equation matrices A^T A — you should not use general LU. The Cholesky factorization writes A = L L^T with L lower-triangular and positive diagonal, exploiting symmetry to use about (1/3)n^3 flops, half of LU, and needs no pivoting because SPD guarantees stability. Solving then uses two triangular sweeps as usual. A bonus: Cholesky is a cheap positive-definiteness test — scipy.linalg.cholesky raises LinAlgError exactly when A is not SPD. Other specialized factorizations match other structure: QR for least squares and orthogonalization, and banded/sparse solvers for matrices with mostly-zero entries. Matching the factorization to the matrix structure is how you get both speed and accuracy.
Worked Example 1
Problem. Compute the Cholesky factor of A = [[4, 2],[2, 3]] by hand.
Answer. L = [[2, 0],[1, sqrt(2)]]; A = L L^T reconstructs [[4,2],[2,3]].
Worked Example 2
Problem. Use Cholesky to detect that A = [[1, 2],[2, 1]] is NOT positive-definite.
Answer. Cholesky raises LinAlgError because A has a negative eigenvalue; it is symmetric but not positive-definite.
Problem. Find the Cholesky factor of A = [[9, 3],[3, 5]].
Solution. l11 = sqrt(9) = 3. l21 = 3/3 = 1. l22 = sqrt(5 - 1^2) = sqrt(4) = 2. So L = [[3,0],[1,2]]. Verify L L^T = [[9,3],[3,5]] = A. Check: `import numpy as np; print(np.linalg.cholesky([[9.,3],[3,5]]))` gives `[[3. 0.] [1. 2.]]`. Answer: L = [[3,0],[1,2]].
Norms assign a single nonnegative size to a vector or matrix, and they are the language for measuring errors and sensitivity.
A vector norm measures length: the 2-norm (Euclidean) is sqrt(sum of squares), the 1-norm sums absolute values, and the infinity-norm takes the max absolute entry. A matrix norm measures how much a matrix can stretch a vector; the induced p-norm is max over x of ||Ax||_p / ||x||_p. The induced 2-norm equals the largest singular value of A, the 1-norm is the max absolute column sum, and the infinity-norm is the max absolute row sum. The Frobenius norm treats the matrix as a long vector (sqrt of sum of squared entries). Norms are the foundation of error analysis: relative error ||x_hat - x|| / ||x|| and the condition number ||A|| ||A^-1|| are both defined through them. Compute them with np.linalg.norm.
Worked Example 1
Problem. Compute the 1-, 2-, and inf-norms of v = [3, -4, 12].
Answer. 1-norm = 19, 2-norm = 13, inf-norm = 12.
Worked Example 2
Problem. Compute the induced 1-, inf-, and 2-norms of A = [[1, 2],[3, 4]].
Answer. 1-norm = 6, inf-norm = 7, 2-norm ≈ 5.465 (the largest singular value).
Problem. For A = [[2, 0],[0, -3]], find the induced 2-norm and the Frobenius norm.
Solution. Singular values are the absolute eigenvalues here: 2 and 3, so the 2-norm (largest singular value) is 3. Frobenius norm = sqrt(2^2 + 0 + 0 + (-3)^2) = sqrt(13) ≈ 3.606. Check: `import numpy as np; A=np.array([[2.,0],[0,-3]]); print(np.linalg.norm(A,2), np.linalg.norm(A))` gives `3.0 3.6055...`. Answer: 2-norm = 3, Frobenius ≈ 3.606.
The condition number kappa(A) bounds how much relative input error is amplified into relative solution error, predicting digit loss.
For Ax = b, the condition number kappa(A) = ||A|| ||A^-1|| (in the 2-norm, the ratio of largest to smallest singular value) bounds the worst-case error amplification: the relative error in x can be up to kappa(A) times the relative error in b or A. A rule of thumb is that you lose about log10(kappa) significant decimal digits when solving in float64, which carries ~16. So kappa ≈ 1e8 risks losing ~8 digits, and kappa near 1/eps ≈ 1e16 means the answer may have no correct digits. Conditioning is a property of the problem, not the algorithm — even a perfect solver cannot beat it. Compute it with np.linalg.cond. The Hilbert matrix is the classic ill-conditioned example, with kappa exploding as its size grows.
Worked Example 1
Problem. Compute the condition number of the 4x4 Hilbert matrix and predict digit loss.
Answer. kappa ≈ 1.55e4, so log10 ≈ 4.2: expect to lose about 4 significant digits.
Worked Example 2
Problem. Contrast a well-conditioned identity-like matrix with an ill-conditioned near-singular one.
Answer. kappa(I) = 1 (ideal); the near-singular 2x2 has kappa ≈ 4e4, amplifying input error ~40000-fold.
Problem. If a system has kappa(A) ≈ 1e6 and the input data is accurate to ~1e-10 relative, roughly how accurate is the solution, and how many digits remain?
Solution. Relative solution error <= kappa * relative input error ≈ 1e6 * 1e-10 = 1e-4. Correct digits ≈ -log10(1e-4) = 4, but you started with ~16 from float64 and lose log10(1e6)=6 just from conditioning. Answer: solution accurate to ~1e-4 relative, about 4 correct significant digits remain.
Production libraries wrap LAPACK: use scipy.linalg for factorizations and Julia's backslash operator for an automatic best-method solve.
In practice you call optimized LAPACK routines rather than hand-rolled loops. SciPy exposes them: scipy.linalg.lu_factor / lu_solve for reusable LU with pivoting, scipy.linalg.cholesky / cho_solve for SPD systems, scipy.linalg.solve for a one-shot solve (with assume_a='pos' or 'sym' to pick the right algorithm), and scipy.linalg.qr for least squares. NumPy's np.linalg.solve covers the common case. Julia makes this idiomatic: the backslash operator A \ b inspects A's structure and dispatches to the appropriate factorization (Cholesky if SPD, LU otherwise, QR for non-square) automatically, and you can precompute factorize(A) to reuse it. The lesson is to let the library choose: pass structural hints so it can select Cholesky over LU, and reuse factorizations across right-hand sides.
Worked Example 1
Problem. Solve an SPD system with scipy.linalg.cho_factor / cho_solve and reuse the factor for two b vectors.
Answer. cho_factor once, cho_solve per right-hand side: efficient SPD solving with the factor reused across b's.
Worked Example 2
Problem. Use scipy.linalg.solve with a structural hint and show the Julia equivalent.
Answer. solve(A, b, assume_a='pos') -> [-0.125, 0.75] via Cholesky; Julia's A \ b auto-selects the same factorization.
Problem. Solve A x = b for A = [[3,1],[1,2]] (SPD) and b = [9, 8] using SciPy's Cholesky path; what is the Julia one-liner?
Solution. `from scipy.linalg import cho_factor, cho_solve; import numpy as np; A=np.array([[3.,1],[1,2]]); c=cho_factor(A); x=cho_solve(c, [9.,8]); print(x)` gives `[2. 3.]`. Verify A@x = [9, 8]. Julia: `[3.0 1.0; 1.0 2.0] \ [9.0, 8.0]` returns `[2.0, 3.0]`. Answer: x = [2, 3].
Build a Hilbert matrix of size n (use scipy.linalg.hilbert) for n from 2 to 12 and solve H x = b for a known true x. Use scipy.linalg.lu_factor / lu_solve and record both the condition number (numpy.linalg.cond) and the relative solution error as n grows. Show that the error tracks the condition number, and contrast with solving a well-conditioned random system.
Deliverable · A script and short report plotting condition number and relative error versus n, plus a sentence relating the observed digit loss to log10 of the condition number.
1. Why solve Ax = b with LU decomposition instead of computing A^-1 and multiplying?
Answer B. Explicitly inverting is more expensive and less accurate; an LU factorization can be reused to solve many right-hand sides via cheap triangular solves.
2. Partial pivoting is used during Gaussian elimination to:
Answer B. Choosing the largest-magnitude pivot in the column prevents division by tiny numbers, which would otherwise magnify rounding error and destabilize the factorization.
3. A condition number of about 1e8 for a float64 system suggests you may lose roughly:
Answer B. Relative solution error is bounded by about kappa times machine epsilon; log10(kappa) approximates the digits lost, so ~8 of the ~16 float64 digits.
4. Cholesky decomposition applies to matrices that are:
Answer B. Cholesky factors A = L L^T and requires symmetry and positive-definiteness; in return it is about twice as fast as general LU.
5. A small residual b - A x_hat guarantees a small solution error only when A is:
Answer C. For ill-conditioned A the residual can be tiny while the actual error is large; only good conditioning links a small residual to a small error.
I can solve Ax = b with LU decomposition and explain why this beats forming the inverse.
I can apply partial pivoting and explain how it improves numerical stability.
I can compute and interpret a condition number to predict how sensitive a solution is to input error.
Interpolation forces a curve through every data point exactly; approximation (regression) seeks a best overall fit that may miss every point.
Interpolation builds a function passing exactly through all n data points — appropriate when the data is essentially noise-free (e.g. table lookups, smooth physical samples). Approximation, typically least-squares regression, fits a lower-dimensional model that minimizes total squared error without touching every point — the right choice when data is noisy, because forcing the curve through noisy points fits the noise. The key trade-off: interpolation has zero residual at the nodes but can oscillate wildly between them (especially high-degree polynomials), while regression smooths through noise but leaves residuals. Numerically, choose by asking whether your data is exact (interpolate) or noisy (regress). Overfitting — using too flexible a model on noisy data — is the failure mode of treating an approximation problem as an interpolation problem.
Worked Example 1
Problem. Given exact points (0,1), (1,2), (2,5), find the interpolating quadratic.
Answer. y = x^2 + 1 passes exactly through all three points (zero residual).
Worked Example 2
Problem. For noisy data y ≈ 2x + noise, contrast interpolating vs. regressing 5 points.
Answer. Interpolation fits the noise with wiggles; the degree-1 regression recovers the true slope ~2 — better for noisy data.
Problem. You have 50 sensor readings with measurement noise and want the underlying linear trend. Interpolate or regress, and why?
Solution. Regress. The data is noisy, so forcing a curve through all 50 points (a degree-49 interpolant) would fit the noise and oscillate uselessly. A least-squares line `m, b = np.polyfit(x, y, 1)` minimizes total squared error and recovers the trend while ignoring noise. Answer: use regression (low-degree least squares) because the data is noisy; interpolation would overfit.
A unique degree n-1 polynomial passes through n points, but on equally spaced nodes high degree causes huge edge oscillations — the Runge phenomenon.
Through n points there is exactly one polynomial of degree at most n-1; it can be built via the Lagrange or Newton form, or fit with np.polyfit. The trap is that raising the degree does not always improve accuracy. On equally spaced nodes, interpolating a smooth function like Runge's f(x) = 1/(1+25x^2) with a high-degree polynomial produces enormous oscillations near the interval ends — the error grows, not shrinks, as the degree increases. The cause is that equally spaced nodes give a poorly conditioned interpolation; the cure is to cluster nodes toward the endpoints using Chebyshev points (denser near +/-1), which tame the oscillation and give near-optimal interpolation. Alternatively, abandon a single global polynomial for piecewise splines.
Worked Example 1
Problem. Fit a degree-2 polynomial through (-1, 1), (0, 0), (1, 1) and evaluate at x = 0.5.
Answer. The interpolating polynomial is p(x) = x^2; p(0.5) = 0.25.
Worked Example 2
Problem. Demonstrate the Runge phenomenon: interpolate 1/(1+25x^2) at 15 equally spaced nodes and check the max error.
Answer. Equally spaced degree-14 interpolation has max error ~7 near the edges; Chebyshev nodes fix it — classic Runge phenomenon.
Problem. Why does interpolating a smooth function at 20 equally spaced points often give a worse result than at 5 points near the edges?
Solution. Higher degree on equally spaced nodes triggers the Runge phenomenon: the interpolating polynomial of degree 19 develops large oscillations near +/-1 because the node distribution makes the problem ill-conditioned. With 5 points the degree is only 4, so the oscillation is mild. The fix is Chebyshev nodes (clustered at the ends) or piecewise splines. Answer: more equally spaced nodes raise the degree and amplify edge oscillations, so the error grows.
Splines interpolate with low-degree polynomial pieces stitched together smoothly, avoiding high-degree oscillation entirely.
Instead of one global high-degree polynomial, a spline uses a separate low-degree polynomial on each interval between knots, joined for smoothness. A linear spline connects points with straight segments (continuous but kinked). The cubic spline is the standard: each piece is a cubic, and the pieces match in value, first derivative, and second derivative at the interior knots, giving a curve that looks smooth to the eye and has no Runge oscillation. The continuity conditions plus boundary conditions (natural: second derivative zero at the ends; or not-a-knot, SciPy's default) determine all the cubic coefficients via a tridiagonal linear system that solves in O(n). Use scipy.interpolate.CubicSpline. Splines are the go-to interpolant for smooth-looking, well-behaved curves through many points.
Worked Example 1
Problem. Build a cubic spline through (0,0), (1,1), (2,0), (3,1) and evaluate at x = 1.5.
Answer. The cubic spline interpolates the points and gives cs(1.5) ≈ 0.344, smoothly between knots.
Worked Example 2
Problem. Compare a cubic spline to a high-degree polynomial on Runge's function (max error).
Answer. The cubic spline's max error is ~0.02 versus ~7 for the degree-14 polynomial — splines defeat the Runge problem.
Problem. Fit a natural cubic spline to (0,1), (1,3), (2,2) and evaluate at x = 0.5.
Solution. `from scipy.interpolate import CubicSpline; import numpy as np; x=np.array([0,1,2.]); y=np.array([1,3,2.]); cs=CubicSpline(x,y,bc_type='natural'); print(cs(0.5))`. The natural spline (zero second derivative at the ends) gives cs(0.5) ≈ 2.0625. It passes through all three points: cs([0,1,2]) = [1,3,2]. Answer: cs(0.5) ≈ 2.06 with the natural boundary condition.
When there are more equations than unknowns, no exact solution exists; least squares finds the x minimizing ||Ax - b||_2.
Fitting a model with p parameters to m > p data points gives an overdetermined system Ax = b that generally has no exact solution. Least squares instead finds the x that minimizes the residual norm ||Ax - b||_2^2 — the closest fit in the Euclidean sense. Geometrically, A x* is the orthogonal projection of b onto the column space of A, so the residual b - A x* is perpendicular to every column of A. This orthogonality gives the normal equations A^T A x = A^T b, whose solution is the least-squares estimate. The columns of A are the basis functions of the model (for a polynomial fit, powers of x — a Vandermonde matrix). NumPy's np.linalg.lstsq solves it stably; the residual sum of squares measures fit quality.
Worked Example 1
Problem. Fit a line y = a x + b to (1,1), (2,2), (3,2) by least squares.
Answer. Best-fit line y ≈ 0.5 x + 0.667; least squares gives slope 0.5, intercept 0.667.
Worked Example 2
Problem. Fit y = c (a single constant) to data [4, 6, 8, 10] by least squares and show it is the mean.
Answer. The least-squares constant fit is c = 7, exactly the mean of the data.
Problem. Set up and solve the least-squares line fit for (0,1), (1,1), (2,3), (3,4).
Solution. Design matrix A = [[0,1],[1,1],[2,1],[3,1]], b = [1,1,3,4]. `import numpy as np; A=np.array([[0,1],[1,1],[2,1],[3,1.]]); b=np.array([1,1,3,4.]); print(np.linalg.lstsq(A,b,rcond=None)[0])` gives `[1.1 0.5]`. So y ≈ 1.1 x + 0.5. Answer: best-fit line y = 1.1 x + 0.5 (slope 1.1, intercept 0.5).
The normal equations are simple but square the condition number; QR factorization solves least squares stably without forming A^T A.
The normal equations A^T A x = A^T b give the least-squares solution directly, but forming A^T A roughly squares the condition number (kappa(A^T A) = kappa(A)^2), so a moderately ill-conditioned A becomes severely so and accuracy collapses. The stable alternative factors A = QR with Q orthogonal and R upper-triangular; since orthogonal transforms preserve norms, ||Ax - b|| = ||Rx - Q^T b||, and you solve R x = Q^T b (the relevant top block) by back substitution — all at condition number kappa(A), not its square. np.linalg.lstsq and scipy.linalg.lstsq use a stable SVD/QR-based method internally, which is why you should call them rather than hand-coding the normal equations. Reserve the normal equations for small, well-conditioned problems.
Worked Example 1
Problem. Solve a least-squares fit via QR for A=[[1,1],[1,2],[1,3]], b=[1,2,2].
Answer. QR gives x = [0.667, 0.5]; it solves R x = Q^T b stably without squaring the condition number.
Worked Example 2
Problem. Show the normal equations lose accuracy where QR does not, on an ill-conditioned fit.
Answer. Normal equations blow up (kappa squared) on the near-collinear design; lstsq's QR/SVD path stays accurate.
Problem. Explain why QR is preferred over the normal equations when A has condition number 1e4, in float64.
Solution. The normal equations form A^T A, whose condition number is kappa(A)^2 = (1e4)^2 = 1e8. In float64 (~16 digits) that risks losing ~8 digits. QR works directly with A at condition number 1e4, risking only ~4 digits lost, because orthogonal transformations preserve the 2-norm and never square the conditioning. Answer: QR keeps kappa at 1e4 (lose ~4 digits) while the normal equations inflate it to 1e8 (lose ~8 digits).
Raising model flexibility lowers training error but risks fitting noise; regularization penalizes large coefficients to restore generalization.
Fitting a degree-d polynomial by least squares (np.polyfit) gets more flexible as d grows. Past a point, extra flexibility fits the noise rather than the signal — overfitting — yielding tiny training error but large error on new data and huge, oscillating coefficients. The standard remedy is regularization: ridge (L2) regression minimizes ||Ax - b||^2 + lambda ||x||^2, shrinking coefficients toward zero and trading a little bias for much lower variance. The penalty lambda is a hyperparameter chosen by validation (e.g. cross-validation). Ridge also improves conditioning: the regularized normal equations (A^T A + lambda I) x = A^T b are always well-posed even when A^T A is singular. Lasso (L1) instead drives some coefficients exactly to zero, doubling as feature selection.
Worked Example 1
Problem. Fit degree-1 and degree-9 polynomials to 10 noisy points and compare coefficient sizes.
Answer. Degree-1 keeps small coefficients (~2); degree-9 has huge oscillating coefficients despite zero training error — overfitting.
Worked Example 2
Problem. Apply ridge regularization to the regularized normal equations and watch coefficients shrink.
Answer. Adding lambda*I shrinks the coefficient magnitudes and makes the system well-conditioned, reducing overfitting.
Problem. You fit a degree-12 polynomial to 13 noisy points and get zero training error but terrible predictions. Diagnose and fix.
Solution. This is overfitting: a degree-12 polynomial interpolates all 13 points exactly (zero residual) but fits the noise, so coefficients are huge and predictions oscillate. Fix by (1) lowering the degree to match the true signal complexity, and/or (2) adding ridge regularization: solve (A^T A + lambda I) x = A^T b with lambda chosen by cross-validation to shrink coefficients. Answer: reduce model flexibility (lower degree) and/or regularize with L2 to restore generalization.
SciPy offers ready interpolators (CubicSpline, interp1d) and nonlinear curve fitting (curve_fit) so you rarely hand-code these.
For interpolation, scipy.interpolate provides CubicSpline (smooth C2 cubic interpolant), interp1d (linear/cubic/nearest with a callable result), and PchipInterpolator (shape-preserving, no overshoot). For approximation, numpy.polyfit handles polynomial least squares, np.linalg.lstsq handles general linear least squares, and scipy.optimize.curve_fit fits arbitrary nonlinear models y = f(x; params) by minimizing squared residuals via Levenberg-Marquardt, returning best-fit parameters and a covariance matrix for uncertainty estimates. The decision tree: noise-free smooth data -> CubicSpline; linear-in-parameters model on noisy data -> polyfit/lstsq; nonlinear model (exponentials, sigmoids) -> curve_fit with a good initial guess p0. Always supply p0 for curve_fit and sanity-check the returned covariance for parameter reliability.
Worked Example 1
Problem. Use scipy.optimize.curve_fit to fit y = a*exp(b*x) to sampled data.
Answer. curve_fit recovers a ≈ 3.0, b ≈ 0.5; pcov's diagonal square-roots give the parameter standard errors.
Worked Example 2
Problem. Build a callable linear interpolant with scipy.interpolate.interp1d and evaluate between samples.
Answer. interp1d returns a callable; linear interp gives f(1.5) = 2.5, cubic gives a smoother value.
Problem. Fit y = a / (1 + exp(-b*(x - c))) (a logistic) to data with scipy.optimize.curve_fit; what must you supply?
Solution. You must give curve_fit the model function and an initial guess p0 for [a, b, c]: `from scipy.optimize import curve_fit; model=lambda x,a,b,c: a/(1+np.exp(-b*(x-c))); popt,pcov=curve_fit(model, x, y, p0=[max(y), 1.0, np.median(x)])`. A reasonable p0 (a near the data max, c near the midpoint) is essential so Levenberg-Marquardt converges to the right minimum. Answer: supply the model and a sensible p0; popt holds the fitted parameters, pcov their covariance.
Sample Runge's function f(x) = 1/(1 + 25 x**2) at equally spaced points and interpolate with both a single high-degree polynomial (numpy.polyfit) and a cubic spline (scipy.interpolate.CubicSpline). Plot both against the true function to expose the Runge oscillations. Then add Gaussian noise to the samples and fit a low-degree least-squares polynomial, comparing the normal-equations solution to numpy.linalg.lstsq.
Deliverable · A script and report with two figures (interpolation comparison and noisy regression fit) and a paragraph on when you would prefer a spline, a regression, or an exact interpolant.
1. The Runge phenomenon is best avoided by:
Answer B. High-degree polynomials on equally spaced nodes oscillate near the edges; piecewise splines or Chebyshev-clustered nodes suppress this.
2. A cubic spline interpolant is constructed to be:
Answer B. Cubic splines join cubic pieces so that the function and its first two derivatives are continuous, giving a smooth interpolant without high-degree oscillation.
3. Solving least squares via the normal equations is numerically risky because it:
Answer B. Forming A^T A roughly squares kappa(A), so a moderately conditioned A can become severely ill-conditioned; QR avoids this by never forming A^T A.
4. Interpolation differs from regression in that interpolation:
Answer B. Interpolation forces the curve through each point exactly, while regression trades exactness for a best overall fit, which is preferable with noisy data.
5. A high-degree polynomial that fits training points perfectly but predicts new points poorly is exhibiting:
Answer B. Fitting the noise as if it were signal is overfitting; it shows low training error but high error on unseen data, often mitigated by regularization.
I can interpolate data with polynomials and splines and explain when high-degree polynomials fail.
I can set up and solve a least-squares problem and explain why QR is more stable than the normal equations.
I can recognize overfitting and choose between interpolation and regression for a given dataset.
Taylor expansion turns a derivative into a difference quotient and reveals exactly how big the truncation error is.
A finite-difference formula approximates a derivative using nearby function values, and the Taylor series both derives the formula and quantifies its error. Expanding f(x+h) = f(x) + h f'(x) + (h^2/2) f''(x) + ... and solving for f'(x) gives the forward difference f'(x) ≈ (f(x+h) - f(x))/h with leading error (h/2) f''(x), so it is first-order accurate, O(h). Combining expansions cleverly cancels low-order terms to raise the order: subtracting f(x-h) from f(x+h) cancels the f'' term, yielding the second-order central difference. The Taylor remainder makes the error term explicit, which is how we predict that halving h shrinks a p-th order error by 2^p. This is the workhorse for numerical derivatives and the building block of ODE and PDE solvers.
Worked Example 1
Problem. Derive the forward-difference formula for f'(x) and its leading error from Taylor's series.
Answer. Forward difference is first-order, O(h); the e^x test gives 1.00501, error ~h/2 = 0.005 as predicted.
Worked Example 2
Problem. Derive the second-order central difference by combining f(x+h) and f(x-h).
Answer. Central difference is second-order, O(h^2); the f' and odd terms reinforce while f'' cancels — error ~h^2/6 f'''.
Problem. Use a forward difference with h = 0.1 to approximate f'(2) for f(x) = x^2, and compare to the exact value.
Solution. Forward difference: (f(2.1) - f(2))/0.1 = (4.41 - 4)/0.1 = 0.41/0.1 = 4.1. Exact f'(2) = 2*2 = 4. The error is 0.1, equal to (h/2)*f'' = (0.1/2)*2 = 0.1, confirming O(h) accuracy. Answer: approximation 4.1, exact 4, error 0.1 = (h/2)*f''.
The three basic stencils trade accuracy for which points they use; central is second-order, forward and backward are first-order.
Three standard first-derivative stencils exist. Forward difference (f(x+h)-f(x))/h uses the point ahead; backward difference (f(x)-f(x-h))/h uses the point behind; both are first-order O(h). The central difference (f(x+h)-f(x-h))/(2h) straddles x symmetrically, cancelling the second-order Taylor term to give second-order O(h^2) accuracy at the same cost. Use forward/backward at domain boundaries where you cannot reach both sides (e.g. the start or end of a time series), and central everywhere in the interior for the extra accuracy. The same idea extends to the second derivative: (f(x+h) - 2f(x) + f(x-h))/h^2 is a second-order central approximation to f''(x). These stencils underpin finite-difference solvers for differential equations.
Worked Example 1
Problem. Approximate f'(1) for f(x)=ln(x) with forward, backward, and central differences, h=0.1; exact is 1.
Answer. Forward 0.953, backward 1.054, central 1.003; central's O(h^2) error (0.003) beats the O(h) one-sided errors (~0.05).
Worked Example 2
Problem. Approximate f''(0) for f(x)=cos(x) with the central second-difference, h=0.01; exact is -1.
Answer. Second central difference gives ≈ -1.0000, matching f''(0) = -1 to ~1e-5 (O(h^2) error).
Problem. For f(x) = x^3 at x = 1 with h = 0.1, compute the central difference and compare to f'(1) = 3.
Solution. Central difference = (f(1.1) - f(0.9))/(2*0.1) = (1.331 - 0.729)/0.2 = 0.602/0.2 = 3.01. Exact f'(1) = 3*1^2 = 3. The error is 0.01 = (h^2/6)*f'''(1) = (0.01/6)*6 = 0.01, confirming O(h^2). Answer: central difference = 3.01, exact 3, error 0.01.
Truncation error falls as h shrinks while round-off rises, so total error is V-shaped with an optimal nonzero h.
Two competing errors govern finite differences. Truncation error comes from cutting off the Taylor series and decreases as h shrinks (like h for forward, h^2 for central). Round-off error comes from cancellation: subtracting two nearly equal function values f(x+h) and f(x), each carrying relative error ~eps, leaves an absolute error ~eps that gets divided by h, so it grows like eps/h as h shrinks. Total error is therefore U/V-shaped: minimized at an optimal h. For the forward difference, balancing (h/2)|f''| against eps|f|/h gives h_opt ≈ sqrt(eps) ≈ 1.5e-8 in float64; for the central difference the balance gives h_opt ≈ eps^(1/3) ≈ 6e-6. Picking h far below these (e.g. 1e-15) makes the answer worse, not better — a classic beginner trap.
Worked Example 1
Problem. Estimate the optimal h for a forward-difference derivative in float64 and explain the scaling.
Answer. Forward-difference h_opt ≈ sqrt(eps) ≈ 1.5e-8, giving a best achievable error of about sqrt(eps).
Worked Example 2
Problem. Show the V curve numerically for the central difference of sin at x=1 across h.
Answer. Error is V-shaped: best near h ≈ eps^(1/3) ≈ 6e-6 for the central difference; smaller h worsens it via round-off.
Problem. For the central difference, why is h = 1e-12 a worse choice than h = 1e-6 in float64?
Solution. Central-difference total error ≈ (h^2/6)|f'''| + eps|f|/h. At h=1e-6 truncation is ~1.7e-13 and round-off ~2.2e-10, near the V's bottom. At h=1e-12 truncation is negligible but round-off is ~eps/h ≈ 2.2e-4, far larger. So h=1e-12 sits on the round-off-dominated side of the V. Answer: h=1e-12 amplifies cancellation round-off (~2e-4) while h=1e-6 is near optimal (~2e-10).
Newton-Cotes rules integrate by replacing the integrand with an interpolating polynomial: lines (trapezoid) or parabolas (Simpson).
Newton-Cotes quadrature approximates an integral by integrating a polynomial that interpolates the integrand at equally spaced points. The trapezoidal rule fits a straight line on each subinterval, giving (h/2)(f0 + f1) per panel and O(h^2) accuracy. Simpson's rule fits a parabola over each pair of subintervals, giving (h/3)(f0 + 4f1 + f2) and O(h^4) accuracy — a big jump because the symmetric parabola also integrates cubics exactly (it has one extra order of precision than its degree suggests). Higher-order Newton-Cotes rules exist (Simpson's 3/8, Boole's) but the weights eventually go negative and become unstable, so in practice trapezoid and Simpson dominate, applied in composite form (next lesson). Simpson requires an even number of subintervals.
Worked Example 1
Problem. Estimate the integral of x^2 from 0 to 2 with the single trapezoidal rule; exact is 8/3.
Answer. Single trapezoid gives 4 vs. exact 8/3 ≈ 2.667; the convex curve makes the line overestimate.
Worked Example 2
Problem. Estimate the same integral of x^2 from 0 to 2 with Simpson's rule (one parabola); exact is 8/3.
Answer. Simpson's rule gives exactly 8/3 ≈ 2.667; it integrates the quadratic with zero error.
Problem. Use Simpson's rule with nodes at 0, 1, 2 to integrate f(x) = x^3 from 0 to 2; exact is 4.
Solution. Simpson: (h/3)(f(0) + 4 f(1) + f(2)) with h=1 = (1/3)(0 + 4*1 + 8) = (1/3)(12) = 4. Exact integral of x^3 from 0 to 2 is 2^4/4 = 4. They match exactly because Simpson is exact for cubics (degree of precision 3). Answer: Simpson gives 4, exactly the true integral.
Composite rules apply the basic stencil on many small panels; the error inherits the panel order, so halving h cuts error by 2^p.
A single trapezoid or Simpson panel over a wide interval is crude; the composite rules split [a,b] into n subintervals of width h = (b-a)/n and sum the basic rule over each, dramatically improving accuracy. The composite trapezoid rule sums (h/2)(f_i + f_{i+1}) and is globally O(h^2); composite Simpson sums parabolic panels and is globally O(h^4). Order of accuracy is the exponent p in error ~ C h^p: it tells you the payoff of refinement. For a p-th order rule, halving h (doubling n) multiplies the error by (1/2)^p — a factor of 4 for trapezoid (p=2) and 16 for Simpson (p=4). You can verify the order empirically by computing the error at h and h/2 and checking their ratio against 2^p; this is the standard convergence test.
Worked Example 1
Problem. Integrate sin(x) from 0 to pi with composite trapezoid at n=4 and n=8; exact is 2.
Answer. Composite trapezoid: 1.896 (n=4), 1.974 (n=8); the 4x error reduction confirms second-order accuracy.
Worked Example 2
Problem. Integrate sin(x) from 0 to pi with composite Simpson at n=4 and n=8; verify O(h^4).
Answer. Composite Simpson: 2.00456 (n=4), 2.000269 (n=8); the ~16x error reduction confirms fourth-order accuracy.
Problem. Composite trapezoid gives error 0.08 at n=10 and 0.02 at n=20 for some integral. What is the order of accuracy?
Solution. Doubling n from 10 to 20 halves h. The error ratio is 0.08/0.02 = 4 = 2^p, so 2^p = 4 gives p = 2. This is second-order accuracy, exactly what the composite trapezoidal rule predicts. Answer: order p = 2 (second-order), since halving h reduced the error by a factor of 4.
Adaptive quadrature concentrates sample points where the integrand is hard, refining only until a local error estimate meets the tolerance.
Uniform composite rules waste effort: they use the same fine spacing everywhere, even where the integrand is nearly straight. Adaptive quadrature instead estimates the local error on each subinterval (typically by comparing a coarse and a fine rule, e.g. Simpson on the whole interval vs. on its two halves) and recursively subdivides only the panels whose estimated error exceeds their share of the tolerance. This places many points near singularities, peaks, or rapid oscillation and few where the function is smooth, reaching a target accuracy with far fewer evaluations. scipy.integrate.quad implements an adaptive Gauss-Kronrod scheme and returns both the integral and an estimated absolute error. It handles infinite limits and integrable singularities far better than fixed-grid rules, which is why it is the default general-purpose integrator.
Worked Example 1
Problem. Integrate the peaked function 1/(1+25x^2) from -1 to 1 with scipy.integrate.quad and read the error estimate.
Answer. quad returns 0.549360 with estimated error ~7e-12, matching (2/5)arctan(5); it auto-refined near the peak.
Worked Example 2
Problem. Integrate a function with an integrable singularity, 1/sqrt(x) from 0 to 1; exact is 2.
Answer. quad handles the 1/sqrt(x) singularity and returns ≈ 2.0; adaptive refinement near x=0 makes it accurate where a fixed grid struggles.
Problem. Use scipy.integrate.quad to integrate exp(-x^2) from 0 to infinity and compare to the known value sqrt(pi)/2.
Solution. `from scipy.integrate import quad; import numpy as np; val, err = quad(lambda x: np.exp(-x**2), 0, np.inf); print(val, err)` gives about (0.8862269254527579, 7.1e-09). The exact value is sqrt(pi)/2 = `np.sqrt(np.pi)/2` ≈ 0.8862269254527580. Answer: quad returns ≈ 0.886227 with tiny error, matching sqrt(pi)/2; it handles the infinite limit via adaptive transformation.
SciPy and NumPy provide ready tools: np.trapz, scipy.integrate.simpson/quad for integrals, and finite-difference helpers for derivatives.
For integration, NumPy's np.trapz integrates sampled data with the composite trapezoidal rule, scipy.integrate.simpson does composite Simpson on samples, and scipy.integrate.quad adaptively integrates a callable function (the default choice when you have f as a function and want accuracy with an error estimate). For higher dimensions there are dblquad/tplquad and nquad. For numerical differentiation, scipy.differentiate.derivative (modern SciPy) or scipy.optimize.approx_fprime computes finite-difference gradients with sensible step control, and numpy.gradient differentiates sampled arrays with second-order central differences in the interior. The practical rule: integrate sampled data with trapz/simpson, integrate a function with quad; and for derivatives prefer a well-implemented helper that manages the optimal step size rather than hand-picking h.
Worked Example 1
Problem. Integrate a function with quad and the same data sampled with np.trapz; compare.
Answer. quad gives 2.0 (machine accuracy); np.trapz on 1001 samples gives 1.9999983 — use quad for callables, trapz for fixed samples.
Worked Example 2
Problem. Differentiate sampled data with numpy.gradient and a callable with scipy.optimize.approx_fprime.
Answer. np.gradient returns the central-difference derivative array; approx_fprime gives the finite-difference gradient of a callable (cos(1) ≈ 0.5403).
Problem. You have measured arrays t and v (velocity) and want total distance traveled. Which SciPy/NumPy tool, and how?
Solution. Distance is the integral of velocity over time, and you have sampled data (not a callable), so use the trapezoidal rule: `import numpy as np; distance = np.trapz(v, t)`, or `from scipy.integrate import simpson; distance = simpson(v, x=t)` for higher accuracy if the samples are smooth and evenly spaced. Use quad only if you had v as a function. Answer: np.trapz(v, t) (or scipy.integrate.simpson(v, x=t)) integrates the sampled velocity to give distance.
For a known function such as f(x) = sin(x), approximate f'(x) at a point with the central-difference formula across step sizes h from 1e-1 down to 1e-12. Plot the absolute error versus h on log-log axes and show the characteristic V shape where truncation error dominates for large h and round-off dominates for small h. Then estimate the integral of f over an interval with numpy.trapz, scipy.integrate.simpson, and scipy.integrate.quad and compare their errors.
Deliverable · A script and report with the error-vs-h log-log plot, the empirically best h, and a small table comparing the three integration methods against the exact integral.
1. The central-difference approximation to a first derivative has truncation error of order:
Answer B. The symmetric central difference cancels the first-order Taylor term, leaving a leading error proportional to h^2 (second-order accurate).
2. As the step size h is made extremely small, the error in a finite-difference derivative eventually:
Answer B. Subtracting nearly equal function values for tiny h causes cancellation; the round-off term grows like 1/h, so total error rises again below an optimal h.
3. Compared with the trapezoidal rule, Simpson's rule on smooth functions is:
Answer B. Simpson fits parabolas and integrates them exactly, giving O(h^4) accuracy versus the trapezoid's O(h^2) for smooth integrands.
4. Adaptive quadrature like scipy.integrate.quad improves efficiency by:
Answer B. Adaptive schemes estimate local error and refine only the hard regions, reaching a tolerance with far fewer evaluations than uniform sampling.
5. A method described as second-order accurate means that halving the step size should reduce the error by roughly a factor of:
Answer B. Second order means error ~ O(h^2); halving h multiplies the error by (1/2)^2 = 1/4, a fourfold reduction.
I can derive finite-difference formulas from Taylor expansions and state their order of accuracy.
I can choose a step size that balances truncation error against round-off error.
I can integrate functions numerically with trapezoid, Simpson, and adaptive quadrature and estimate the error.
An initial value problem pairs an ODE y' = f(t, y) with a starting point; its solution is the trajectory the system follows from that point.
An initial value problem (IVP) asks: given the rate of change y' = f(t, y) and a starting state y(t0) = y0, what is y(t) for later times? The function f defines a slope field; the solution is the unique curve threading that field from the initial point, often called the flow of the system. Existence and uniqueness are guaranteed when f is Lipschitz continuous in y (a smoothness condition). Higher-order ODEs are handled by rewriting them as a first-order system: y'' = g(t, y, y') becomes a vector ODE for the state [y, y']. Numerical solvers march forward in steps, evaluating f to advance the state; they never produce the closed-form solution, only sampled trajectory points. Casting any model as y' = f(t, y) with a state vector is the first step before calling a solver.
Worked Example 1
Problem. Write the second-order ODE y'' + y = 0, y(0)=1, y'(0)=0 as a first-order system.
Answer. State u=[y,y'], f(t,u)=[u[1],-u[0]], u(0)=[1,0]; the true solution is y(t)=cos(t).
Worked Example 2
Problem. Verify that y(t) = e^{-2t} solves the IVP y' = -2y, y(0) = 1.
Answer. y(t)=e^{-2t} satisfies both y'=-2y and y(0)=1, so it is the unique solution of the IVP.
Problem. Convert the pendulum equation theta'' = -sin(theta) into a first-order system suitable for solve_ivp.
Solution. Introduce the state u = [theta, omega] where omega = theta'. Then u[0]' = omega and u[1]' = theta'' = -sin(theta) = -sin(u[0]). So f(t, u) = [u[1], -np.sin(u[0])]. In code: `import numpy as np; f = lambda t, u: [u[1], -np.sin(u[0])]`, with an initial state like u0 = [theta0, omega0]. Answer: state [theta, omega], f(t,u) = [u[1], -sin(u[0])].
Euler's method takes a single tangent-line step per interval; its per-step (local) error is O(h^2) but the accumulated (global) error is only O(h).
The forward Euler method is the simplest IVP solver: from (t_n, y_n) it follows the slope f for one step, y_{n+1} = y_n + h f(t_n, y_n). Geometrically it walks along the tangent line, so it drifts off curved solutions. The local truncation error — the error made in one step from exact data — is O(h^2), from the dropped Taylor term. But over the interval [t0, T] you take N = (T-t0)/h steps, and the per-step errors accumulate, so the global error is O(h): Euler is first-order. This one-power drop from local to global order is a general pattern. First-order accuracy is poor — halving h only halves the error — which motivates the higher-order Runge-Kutta methods. Euler is also only conditionally stable, limiting h for stiff problems.
Worked Example 1
Problem. Take two Euler steps for y' = y, y(0) = 1 with h = 0.5; compare to e^1.
Answer. Euler gives y(1) ≈ 2.25 vs. exact 2.718; the big first-order error reflects the crude tangent steps.
Worked Example 2
Problem. Empirically confirm Euler is first-order on y' = -y, y(0)=1 by halving h.
Answer. Errors ~0.017 (h=0.1) and ~0.009 (h=0.05), ratio ~2 — first-order global accuracy confirmed.
Problem. Take one Euler step for y' = t + y, y(0) = 1 with h = 0.1, and give y(0.1).
Solution. Euler: y_1 = y_0 + h*f(t_0, y_0) = 1 + 0.1*(0 + 1) = 1 + 0.1 = 1.1. So the estimate at t=0.1 is y(0.1) ≈ 1.1. (The exact solution is y(t) = 2e^t - t - 1, giving y(0.1) ≈ 1.1103, so the one-step error is ~0.01 = O(h^2).) Answer: y(0.1) ≈ 1.1.
Second-order Runge-Kutta methods (midpoint, Heun) sample the slope twice per step to cancel Euler's leading error term.
Euler uses only the slope at the start of the interval, which is biased for curved solutions. Second-order Runge-Kutta methods evaluate f twice to get a better average slope, achieving global O(h^2) accuracy. The midpoint method takes a half Euler step to t_n + h/2, evaluates the slope there, and uses that midpoint slope for the full step. Heun's method (improved Euler) is a predictor-corrector: it predicts y* with an Euler step, evaluates the slope at the endpoint, and corrects using the average of the start and end slopes — k1 = f(t_n, y_n), k2 = f(t_n+h, y_n+h k1), y_{n+1} = y_n + (h/2)(k1 + k2). Both cost two function evaluations per step but halving h now reduces error fourfold, a clear improvement over Euler for the modest extra cost.
Worked Example 1
Problem. Take one Heun step for y' = y, y(0)=1, h=0.5; compare to Euler and to e^0.5.
Answer. Heun gives 1.625 vs. exact 1.6487 (error 0.024), much closer than Euler's 1.5 — second-order pays off.
Worked Example 2
Problem. Take one midpoint-method step for y' = -2y, y(0)=1, h=0.1; exact y(0.1)=e^{-0.2}.
Answer. Midpoint gives 0.82 vs. exact 0.81873 (error ~0.001), ten times better than Euler at the same h.
Problem. Take one Heun step for y' = t - y, y(0) = 1 with h = 0.2.
Solution. k1 = f(0, 1) = 0 - 1 = -1. Predictor: y* = 1 + 0.2*(-1) = 0.8. k2 = f(0.2, 0.8) = 0.2 - 0.8 = -0.6. Corrector: y_1 = 1 + (0.2/2)(-1 + -0.6) = 1 + 0.1*(-1.6) = 1 - 0.16 = 0.84. (Exact y(0.2) ≈ 0.8375.) Answer: Heun gives y(0.2) ≈ 0.84.
The classic RK4 evaluates the slope four times per step in a weighted average, achieving fourth-order accuracy — the standard non-stiff workhorse.
The classic fourth-order Runge-Kutta method (RK4) samples the slope at four points and combines them with carefully chosen weights so that error terms through h^4 cancel: k1 = f(t_n, y_n); k2 = f(t_n+h/2, y_n+h k1/2); k3 = f(t_n+h/2, y_n+h k2/2); k4 = f(t_n+h, y_n+h k3); then y_{n+1} = y_n + (h/6)(k1 + 2k2 + 2k3 + k4). The midpoint slopes k2, k3 get double weight, mimicking Simpson's rule. RK4 is globally O(h^4), so halving h reduces error sixteenfold, at the cost of four function evaluations per step. This excellent accuracy-to-cost ratio makes it the default for smooth, non-stiff problems. Higher-order and embedded (adaptive) Runge-Kutta variants build on the same multi-stage idea.
Worked Example 1
Problem. Take one RK4 step for y' = y, y(0)=1, h=0.5; compare to e^0.5.
Answer. RK4 gives 1.64844 vs. exact 1.64872 (error ~3e-4) in one step; fourth-order accuracy in action.
Worked Example 2
Problem. Implement RK4 for y' = -y, y(0)=1 over [0,1] and check fourth-order convergence by halving h.
Answer. Errors ~8e-7 (h=0.1) and ~5e-8 (h=0.05), ratio ~16 — fourth-order global accuracy confirmed.
Problem. Why does RK4 typically beat Euler even when both are given the same total number of function evaluations?
Solution. RK4 uses 4 evaluations per step but is O(h^4); Euler uses 1 per step but is O(h). For the same total evaluations, Euler takes 4x more steps (h/4), giving error ~ (h/4) = h/4. RK4 keeps step h with error ~ h^4. For small h, h^4 is vastly smaller than h/4, so RK4 wins decisively. Answer: RK4's fourth-order accuracy overwhelms its 4x per-step cost; its error (h^4) shrinks far faster than Euler's (h).
Adaptive solvers estimate each step's local error from two embedded methods and grow or shrink h to keep it near a tolerance.
A fixed step size is wasteful: it uses tiny steps everywhere even where the solution is smooth. Adaptive Runge-Kutta methods estimate the local error cheaply by computing two solutions of different order from shared stage evaluations — an embedded pair like RK45 (Dormand-Prince) produces a 4th- and 5th-order estimate from the same k's. The difference between them approximates the local error. If that error exceeds the tolerance (atol + rtol*|y|), the step is rejected and retried with a smaller h; if it is comfortably below, h is enlarged for the next step. This automatically takes small steps through rapid transients and large steps through smooth stretches, hitting the requested accuracy with minimal evaluations. scipy.integrate.solve_ivp uses this with method='RK45' by default; you set rtol and atol.
Worked Example 1
Problem. Solve y' = -y, y(0)=1 on [0,5] with scipy.integrate.solve_ivp and inspect the adaptive steps.
Answer. solve_ivp (RK45) hits y(5) ≈ e^{-5} = 0.006738 with few adaptive steps, enlarging h as the solution smooths out.
Worked Example 2
Problem. Tighten tolerances and observe more steps and higher accuracy.
Answer. Tighter rtol means more adaptive steps and smaller error; the solver self-adjusts h to meet the requested tolerance.
Problem. In scipy.integrate.solve_ivp, what do rtol and atol control, and how is the per-step error test formed?
Solution. rtol (relative tolerance) and atol (absolute tolerance) set the accuracy target. The solver accepts a step when the estimated local error for each component is below atol + rtol*|y|. rtol scales the allowed error with the solution's magnitude (good for large values), while atol provides a floor so components near zero are not held to an impossible relative standard. If the estimate exceeds this, the step is rejected and h is reduced. Answer: rtol/atol define the threshold atol + rtol*|y| that the embedded error estimate must satisfy each step.
Stiff systems have fast and slow components together; explicit methods need impractically tiny steps, so implicit solvers are required.
A stiff ODE has dynamics on widely separated time scales — a fast-decaying transient alongside slow evolution. Explicit methods (Euler, RK45) have a bounded stability region, so even after the fast transient has died, they must keep h tiny just to stay stable, making them painfully slow. Implicit methods solve for y_{n+1} on both sides of the update equation, e.g. backward Euler y_{n+1} = y_n + h f(t_{n+1}, y_{n+1}), which requires solving a (possibly nonlinear) system each step but is unconditionally stable (A-stable) for decaying problems, allowing large steps. The cost per step is higher (a Newton solve, needing the Jacobian), but far fewer steps are needed. SciPy's solve_ivp offers implicit methods 'BDF' and 'Radau' for stiff problems; the explicit 'RK45' will crawl or fail.
Worked Example 1
Problem. Show backward Euler is stable for the test equation y' = -50y, y(0)=1, h=0.1, where forward Euler is not.
Answer. Forward Euler factor 4 (unstable, blows up); backward Euler factor 1/6 (stable) — implicit wins for stiff decay.
Worked Example 2
Problem. Solve a stiff system with solve_ivp using 'BDF' and note why 'RK45' struggles.
Answer. BDF solves the stiff system efficiently; RK45 needs far more steps (stability-limited by the fast mode) — use implicit methods when stiff.
Problem. For y' = -1000(y - cos(t)) - sin(t), why will solve_ivp with method='RK45' be slow, and what should you use?
Solution. The factor -1000 makes the system stiff: deviations from cos(t) decay on a time scale ~1/1000, while the solution itself varies slowly with t. Explicit RK45's stability region forces h ~ 1/1000 throughout, so it takes enormous numbers of tiny steps even where the solution is smooth. Use an implicit stiff solver: method='BDF' or method='Radau', which are stable at large h. Answer: RK45 is stability-limited by the -1000 mode; use method='BDF' or 'Radau'.
solve_ivp is the modern unified interface: choose a method, set tolerances, request dense output and event detection.
scipy.integrate.solve_ivp(fun, t_span, y0, method=..., t_eval=..., rtol=..., atol=..., events=..., dense_output=...) is the recommended IVP solver. fun(t, y) returns the derivative vector, t_span = (t0, tf), and y0 is the initial state. Pick the method by character: 'RK45' (default) or 'DOP853' for smooth non-stiff problems, 'BDF' or 'Radau' for stiff ones, 'LSODA' to auto-switch. Use t_eval to get the solution at specific points, dense_output=True to obtain a continuous interpolant sol.sol(t), and events for stopping or recording when a condition (like a zero crossing) occurs. The returned object exposes sol.t, sol.y (shape n_states x n_times), sol.success, and sol.t_events. Always check sol.success and read the message if it failed.
Worked Example 1
Problem. Solve the harmonic oscillator y'' + y = 0, y(0)=1, y'(0)=0 and evaluate at chosen times.
Answer. sol.y[0] ≈ [1, 0, -1, 0, 1], matching cos(t); t_eval pins the output to the requested times.
Worked Example 2
Problem. Use an event to detect when a falling object reaches the ground (height = 0).
Answer. solve_ivp stops at the ground via the terminal event; sol.t_events[0] ≈ 4.515 s, matching sqrt(200/9.81).
Problem. Solve y' = -2y + 1, y(0) = 0 on [0, 3] with solve_ivp and compare the endpoint to the exact solution.
Solution. The exact solution is y(t) = 0.5(1 - e^{-2t}). `from scipy.integrate import solve_ivp; import numpy as np; sol = solve_ivp(lambda t, y: -2*y + 1, [0, 3], [0.0], rtol=1e-9); print(sol.y[0, -1], 0.5*(1-np.exp(-6)))` gives about (0.4987..., 0.4987...). The numerical endpoint matches the exact value 0.5(1 - e^{-6}) ≈ 0.49876. Answer: y(3) ≈ 0.4988, matching the closed-form solution.
Take an IVP with a known closed-form solution (for example y' = -y, y(0) = 1) and implement Euler's method and classic RK4 yourself. Plot the global error versus step size on log-log axes for both and confirm the expected slopes (order 1 vs. order 4). Then re-solve the problem with scipy.integrate.solve_ivp using an explicit method, and apply it to a stiff system to show why a stiff solver (such as 'Radau' or 'BDF') is needed.
Deliverable · A script and report with the order-of-accuracy plot, the measured convergence slopes, and a short explanation of when you switched solvers for the stiff system.
1. Euler's method is which order of accuracy in its global error?
Answer B. Euler has local error O(h^2) per step but accumulates over O(1/h) steps to give first-order global error O(h).
2. The classic RK4 method is popular mainly because it:
Answer B. RK4 uses four stages per step to reach O(h^4) global accuracy, an excellent balance of cost and accuracy for non-stiff problems.
3. A stiff ODE system is characterized by:
Answer B. Stiffness comes from fast and slow components coexisting; explicit methods must use tiny steps for stability, so implicit solvers are preferred.
4. Adaptive step-size solvers decide how big each step should be by:
Answer B. Embedded methods compare two estimates to gauge local error, then adjust h up or down to stay near the requested tolerance efficiently.
5. For a stiff problem, which scipy.integrate.solve_ivp method is appropriate?
Answer C. BDF and Radau are implicit methods designed for stiff systems; the explicit RK options (RK45, RK23, DOP853) struggle with stiffness.
I can implement Euler and Runge-Kutta methods and explain the difference between local and global error.
I can recognize a stiff system and explain why an implicit solver is needed.
I can solve an initial value problem with scipy.integrate.solve_ivp and choose an appropriate method.
Stability is a property of an algorithm — whether it amplifies rounding errors — distinct from conditioning, which is a property of the problem.
Revisiting Unit 1's vocabulary: conditioning describes how sensitive a problem's true answer is to input perturbations (it is intrinsic to the problem), while stability describes whether an algorithm amplifies the rounding errors it makes along the way. A stable (specifically backward-stable) algorithm produces the exact answer to a slightly perturbed input, so its forward error is bounded by roughly condition number times machine epsilon — the best one can hope for. An unstable algorithm injects extra error growth beyond what the conditioning forces. The practical test: when an answer is bad, ask whether the problem is ill-conditioned (nothing to fix in the code) or the algorithm is unstable (reformulate it). The same mathematical result can have several algorithms of differing stability — e.g. classical vs. modified Gram-Schmidt, or the naive vs. stable quadratic formula.
Worked Example 1
Problem. Show classical Gram-Schmidt loses orthogonality where modified Gram-Schmidt (or QR) does not.
Answer. Classical GS yields a non-orthogonal Q on ill-conditioned input; modified GS / np.linalg.qr keep Q^T Q ≈ I — same math, different stability.
Worked Example 2
Problem. Distinguish conditioning from stability for a sample sum and a sample subtraction.
Answer. The cancellation in 1-cos is a conditioning issue; the naive variance formula is an algorithm-stability issue fixed by a two-pass method (np.var).
Problem. An algorithm gives 12 correct digits on a problem with condition number ~1e10 in float64. Is the algorithm stable?
Solution. float64 carries ~16 digits. With condition number ~1e10, even a perfectly backward-stable algorithm loses about log10(1e10) = 10 digits, leaving ~6 correct. Getting 12 correct digits would be better than the conditioning allows, which is impossible — so either the condition number is overestimated or the problem is locally better conditioned. A stable algorithm should give ~6 correct digits here; 12 is suspiciously good, not a sign of instability. Answer: a stable algorithm should yield ~6 correct digits (16 - 10); 12 correct digits exceeds what kappa=1e10 permits, so re-examine the conditioning estimate.
Each Python loop iteration pays interpreter dispatch, dynamic type checks, and object boxing overhead that compiled array operations avoid.
Python is interpreted and dynamically typed, so a simple a + b inside a loop is expensive: the interpreter dispatches bytecode, checks the runtime types of a and b, finds the right __add__, allocates a new boxed object for the result, and manages reference counts — all per element. A C loop over a typed array does none of this; it runs a single machine instruction per element on contiguous memory. That overhead is why a pure-Python numeric loop can be 10-100x slower than the equivalent compiled operation. NumPy sidesteps it by storing data in homogeneous, contiguous C arrays and performing whole-array operations in compiled loops (ufuncs), paying the per-element overhead once at the C level instead of once per Python iteration. The fix for slow numeric Python is almost always to push the loop down into compiled code via vectorization.
Worked Example 1
Problem. Time summing a million numbers with a Python loop versus np.sum.
Answer. np.sum is ~100x faster than the Python loop; the speedup is the eliminated per-element interpreter and boxing overhead.
Worked Example 2
Problem. Explain why list-of-floats arithmetic is slower and more memory-heavy than a NumPy array.
Answer. Lists hold scattered boxed objects (slow, memory-heavy); NumPy arrays hold contiguous typed values (fast, compact), enabling compiled loops.
Problem. A colleague's element-by-element Python loop computing c[i] = a[i]*b[i] for a million elements is slow. What is the root cause and the fix?
Solution. The root cause is per-iteration interpreter overhead: each step dispatches bytecode, type-checks a[i] and b[i], finds __mul__, and allocates a boxed result — millions of times. The fix is vectorization: with NumPy arrays, `c = a * b` performs the whole multiplication in one compiled C loop over contiguous memory, typically 50-100x faster. Answer: interpreter/dynamic-typing overhead per element; replace the loop with the vectorized c = a * b.
Vectorization expresses computation as whole-array operations (ufuncs) that run in compiled loops, replacing explicit Python iteration.
Vectorization means writing array expressions like a*b + np.sin(c) instead of element loops; NumPy executes each as a universal function (ufunc) — a compiled, type-specialized loop over contiguous memory. ufuncs apply elementwise (np.add, np.exp, np.maximum), support reductions (sum, mean, max along an axis), and accumulate (np.cumsum). The wins are speed (compiled C, often SIMD-vectorized) and clarity. Reductions take an axis argument to collapse specific dimensions. Replace conditional loops with np.where(cond, a, b) or boolean-mask indexing. Many problems that look inherently sequential can be vectorized by reformulating them as array algebra — pairwise distances via broadcasting, running statistics via cumulative ufuncs. The mental shift is from 'loop over elements' to 'operate on whole arrays at once.'
Worked Example 1
Problem. Vectorize the computation of sum of squares of 1..n and compare to a loop.
Answer. Both give 333833500; the vectorized (a*a).sum() does it in compiled ufuncs with no Python loop.
Worked Example 2
Problem. Replace an if/else loop that clips negatives to zero with a vectorized ufunc.
Answer. np.maximum(x, 0) (or np.where(x<0,0,x)) gives [0,3,0,5] in one vectorized call, replacing the if/else loop.
Worked Example 3
Problem. Vectorize column-wise normalization of a matrix (subtract column mean, divide by column std).
Answer. (X - X.mean(axis=0)) / X.std(axis=0) standardizes every column at once via reductions plus broadcasting.
Problem. Vectorize computing the Euclidean norm of each row of a matrix X (shape m x n).
Solution. Square elementwise, sum along the columns (axis=1), and take the square root — all vectorized: `import numpy as np; norms = np.sqrt((X**2).sum(axis=1))`, or more directly `norms = np.linalg.norm(X, axis=1)`. For X = [[3,4],[6,8]] this gives [5, 10]. No Python loop is needed; the reduction over axis=1 collapses each row to one number. Answer: np.linalg.norm(X, axis=1) (or np.sqrt((X**2).sum(axis=1))).
Broadcasting lets NumPy combine arrays of different shapes by virtually stretching size-1 dimensions, eliminating loops without copying data.
Broadcasting is NumPy's rule for operating on arrays of different shapes. Aligning shapes from the right, two dimensions are compatible if they are equal or one of them is 1; a size-1 (or missing) dimension is virtually stretched to match the other. No data is copied — the stretched axis is read with a zero stride — so it is both fast and memory-light. This turns nested loops into single expressions: add a (m,1) column vector to a (1,n) row vector to get an (m,n) outer sum; subtract a (n,) mean from an (m,n) matrix to center it. Insert new axes with arr[:, None] or np.newaxis to control how shapes line up. The classic use is pairwise computations: X[:, None, :] - Y[None, :, :] forms all difference vectors at once. Mismatched shapes that satisfy no rule raise a ValueError.
Worked Example 1
Problem. Build a multiplication table (outer product) for 1..4 without loops.
Answer. a[:,None]*a[None,:] broadcasts (4,1) and (1,4) into the (4,4) multiplication table.
Worked Example 2
Problem. Compute all pairwise squared distances between rows of X (m x d) and Y (n x d) via broadcasting.
Answer. X[:,None,:]-Y[None,:,:] then sum over axis 2 yields the full (m,n) pairwise squared-distance matrix without any loop.
Problem. You have a (1000, 3) array of points and want to subtract a single (3,) centroid from every point. How does broadcasting handle it?
Solution. Just write `centered = points - centroid`. Aligning from the right, points is (1000, 3) and centroid is (3,), treated as (1, 3). The size-1 leading dimension is stretched to 1000, so the (3,) centroid is virtually subtracted from each of the 1000 rows — no loop, no copy of the centroid. Answer: points - centroid broadcasts the (3,) centroid across all 1000 rows automatically.
Array element order (C vs. Fortran) and the view/copy distinction govern cache efficiency and whether writes alias the original data.
NumPy arrays store elements in a flat buffer with a memory layout: C-order (row-major, the default) keeps each row contiguous; Fortran-order (column-major) keeps each column contiguous. Accessing memory contiguously is far faster because the CPU loads whole cache lines, so iterate along the contiguous axis (rows for C-order) and prefer operations that stride sequentially. A slice or reshape often returns a view — a new array object sharing the original buffer — so writing through it modifies the parent (aliasing); fancy/boolean indexing returns a copy. Check arr.flags['C_CONTIGUOUS'] and use arr.base or np.shares_memory to tell views from copies. Operations on non-contiguous arrays may silently copy. The performance lesson: keep data contiguous in the access direction, and the correctness lesson: know when you hold a view that can mutate shared data.
Worked Example 1
Problem. Show that a slice is a view that aliases the parent, while fancy indexing copies.
Answer. Slicing returns an aliasing view (writes hit the parent); fancy indexing returns an independent copy.
Worked Example 2
Problem. Demonstrate that summing along the contiguous axis is faster than along the strided axis.
Answer. A.sum(axis=1) (along contiguous rows) is faster than A.sum(axis=0) for a C-ordered array, due to cache-line reuse.
Problem. You take b = a.reshape(3, 4) and set b[0, 0] = 100; a is a 12-element array. Does a change, and why?
Solution. Yes, a changes: reshape returns a view (when the data is contiguous, as a fresh arange is), so b shares a's buffer. Writing b[0,0]=100 sets the first element of the shared buffer, which is a[0]. You can confirm with `np.shares_memory(a, b)` returning True. To avoid aliasing, use a.reshape(3,4).copy(). Answer: yes — reshape returns a view sharing memory, so a[0] becomes 100.
Measure before optimizing: time whole runs (timeit), find hot functions (cProfile), and inspect line-level costs to target the real bottleneck.
Optimization without measurement wastes effort on code that is not the bottleneck. Benchmarking times a snippet repeatedly to get a stable estimate: use timeit (or %timeit in IPython) for microbenchmarks, taking the minimum to reduce noise. Profiling attributes time across a whole program: cProfile reports per-function call counts and cumulative time, revealing which functions dominate; line_profiler (the @profile / kernprof tool) shows time per source line for the suspected hotspot. The standard workflow is profile -> find the function consuming most time -> drill into its lines -> optimize only that -> re-measure to confirm a real speedup. Beware microbenchmark pitfalls: warm up caches, avoid timing one-off allocations, and use realistic input sizes. Premature optimization guided by intuition rather than data is the classic mistake.
Worked Example 1
Problem. Benchmark two ways to square an array with timeit and pick the faster.
Answer. timeit shows a*a and a**2 are nearly identical; the data, not intuition, settles which to use.
Worked Example 2
Problem. Use cProfile to find the dominant function in a small pipeline.
Answer. cProfile's cumulative-time column flags slow() as the bottleneck, directing optimization to the function that actually costs the most.
Problem. Your program is slow and you suspect a particular function. What is the correct sequence of tools and steps?
Solution. First run cProfile on the whole program (cProfile.run('main()')) and read the cumulative-time column to confirm which function truly dominates — do not assume. Then apply line_profiler (@profile with kernprof) to that one function to find the costliest lines. Optimize only those lines (often by vectorizing), then re-run timeit/cProfile to verify the speedup is real. Answer: profile with cProfile to locate the hot function, drill in with line_profiler, optimize that, then re-measure to confirm.
When an algorithm is inherently sequential or hard to vectorize, JIT compilation (Numba) or a compiled language (Julia) recovers C-level speed.
Vectorization solves most performance problems, but some algorithms resist it — tight sequential dependencies (recurrences, iterative solvers with branching), element-wise control flow, or cases where the vectorized form would allocate huge temporaries. Then you reach for compilation. Numba JIT-compiles a Python function to machine code with the @njit decorator; written as plain loops over NumPy arrays, the compiled loop runs at near-C speed without leaving Python, and supports parallelism and SIMD. Julia takes this further: it compiles every function by default via LLVM, so idiomatic loops are already fast and there is no two-language split. The decision order is: vectorize with NumPy first; if that is impossible or memory-prohibitive, JIT the hotspot with Numba; if the whole project is performance-critical numerics, consider writing it in Julia. Always profile to confirm the hotspot before compiling.
Worked Example 1
Problem. Speed up a sequential recurrence with Numba's @njit.
Answer. @njit compiles the sequential recurrence to machine code, giving a large speedup where vectorization is impossible.
Worked Example 2
Problem. Decide between vectorization, Numba, and Julia for three scenarios.
Answer. Vectorize the array math; Numba-JIT the branchy sequential hotspot; choose Julia when the whole numerics project demands compiled speed.
Problem. An iterative solver with data-dependent branching is your profiled hotspot and cannot be vectorized without enormous temporaries. What do you do?
Solution. Since vectorization is ruled out (it would allocate huge temporaries) and profiling already confirmed this is the hotspot, JIT-compile it. Wrap the function with Numba's @njit and write it as plain loops over NumPy arrays; the data-dependent branches stay, but the loop now runs at near-C speed. If the entire project is performance-critical, porting to Julia (compiled by default) is the longer-term option. Answer: apply Numba @njit to the hotspot (loops over NumPy arrays); consider Julia if the whole codebase needs compiled speed.
Take a numeric kernel written as an explicit Python double loop (for example a pairwise distance matrix or a finite-difference stencil) and rewrite it using NumPy vectorization and broadcasting. Benchmark both versions with timeit across growing input sizes and report the speedup. Verify the two implementations agree with numpy.allclose, and discuss whether memory layout or broadcasting created any hidden array copies.
Deliverable · A script and report with a timing table or plot of loop vs. vectorized runtime, the measured speedup factor, and one paragraph on a case where vectorizing would cost too much memory.
1. A numerically stable algorithm is one that:
Answer B. Stability is about error growth during the computation; conditioning is a property of the problem, while stability is a property of the algorithm.
2. NumPy vectorization is faster than an equivalent Python loop mainly because it:
Answer B. Each Python loop iteration carries interpreter and dynamic-typing overhead; ufuncs push the loop into optimized compiled code over contiguous memory.
3. Broadcasting in NumPy lets you:
Answer A. Broadcasting virtually stretches dimensions of size 1 to match, enabling elementwise operations across different shapes without materializing large temporaries.
4. Before optimizing performance, the recommended first step is to:
Answer B. Optimization should be guided by measurement; profiling reveals where time is actually spent so effort targets the real hotspot.
5. A NumPy slice that returns a view rather than a copy means that:
Answer A. A view shares the parent's underlying buffer, so writes through the view change the original; this saves memory but can cause aliasing surprises.
I can distinguish a stable algorithm from an unstable one and explain the difference from conditioning.
I can rewrite scalar Python loops as vectorized NumPy operations and measure the speedup.
I can profile code to find bottlenecks and decide when compiled acceleration is worth it.
Where this leads
Part of the open-source Crunch Academy tech-education path. Pair it with the data-science and ML tracks when you're ready.
All Math for Data Science courses Crunch Academy home