# HW4 Companion — Least Squares, Interpolation, Newton for Systems & 2-D Integration > [!info]- How to use this companion > This is a scaffold, not a solution key. For each problem you get: > - **What to hand in** — what is a *written* calculation, what is a *code* calculation, and what to *report*. > - **Hints** — a ladder of collapsed callouts. Open **Hint 1** first; the next only if still stuck. > - **Code guidance** — a callout per language. Much of the machinery you built in lecture or class code > (Lagrange interpolation, `cubicInterp`/`cubicInterpPeriodic`, the finite-difference gradient/Hessian > Newton step, `mySimp`, a built-in integrator); the code here shows how to *drive* it, with `# TODO` > marks where the real work goes. > > Submit **one written PDF** (calculations, tables, plots, reflections) plus your **code files** and **figures**. Label every axis. --- ## Deliverables at a glance | # | Topic | By hand | In code | Report | |---|-------|---------|---------|--------| | 1 | Least squares (parallel lines) | set up $A,\vec b$; normal equations; the singularity | `lstsq`/pinv, residual | why singular, the $y=x$ midline, min-norm point | | 2 | Steam-table Lagrange | the basis polynomials, value at 210 | build & plot $L(x)$ | $h(210)$, the plot, big-table reflection | | 3 | Cubic spline end conditions | (read the BC rows) | natural + periodic spline of a closed curve | seam plot, corner + error, refinement table | | 4 | Newton for systems | (classification rule) | FD gradient/Hessian Newton | the surface plot, two critical points + types | | 5 | Double integral | separability check | built-in + nested Simpson | error, timing, $n^2$ discussion | --- ## Problem 1 — Least squares of two parallel lines **What to hand in.** - **By hand:** write the system $A\vec z=\vec b$, form the normal equations, and show $A^{\mathsf T}A$ is singular. - **In code:** a least-squares / pseudoinverse solve, reporting the minimum-norm point and the residual. - **Report:** which set of points minimizes the residual, where it sits relative to the two lines, and what the singular normal matrix means. > [!hint]- Hint 1 — turn the lines into a system > Each line is one linear equation in the unknown point $(x,y)$: $y=x+1\Rightarrow x-y=-1$ and $y=x-1\Rightarrow x-y=1$. Stack them: $A=\left[\begin{smallmatrix}1&-1\\1&-1\end{smallmatrix}\right]$, $\vec b=\left[\begin{smallmatrix}-1\\1\end{smallmatrix}\right]$. The two rows are identical on the left but disagree on the right — that is why no $(x,y)$ solves both. > [!hint]- Hint 2 — the normal equations, and the surprise > $A^{\mathsf T}A=\left[\begin{smallmatrix}2&-2\\-2&2\end{smallmatrix}\right]$ and $A^{\mathsf T}\vec b=\left[\begin{smallmatrix}0\\0\end{smallmatrix}\right]$. Compute $\det(A^{\mathsf T}A)$: it is $0$. A singular normal matrix means least squares has **no unique** answer here. > [!hint]- Hint 3 — read off the solution set > The one surviving equation is $2x-2y=0$, i.e. $x=y$. So every point on the line $y=x$ is a least-squares solution — and that line is exactly halfway between $y=x+1$ and $y=x-1$. The best the fit can do is split the difference between the offsets $+1$ and $-1$, landing on $0$; the residual is $\sqrt2$. > [!hint]- Hint 4 — what the code returns > A least-squares routine on a rank-deficient system returns the **minimum-norm** solution: the point on $y=x$ closest to the origin, namely $(0,0)$. The residual $\lVert A\vec z-\vec b\rVert=\sqrt2$ no matter which point on the line you pick. > [!example]- Python > > ```python > import numpy as np > A = np.array([[1., -1.], [1., -1.]]); b = np.array([-1., 1.]) > print(np.linalg.det(A.T @ A)) # 0 -> singular > z, *_ = np.linalg.lstsq(A, b, rcond=None) > print(z, np.linalg.norm(A @ z - b)) # ~[0,0], residual sqrt(2) > ``` > *Look for:* `det(A.T@A)` is $0$; `z` is essentially $(0,0)$; residual $\approx1.4142$. > [!example]- R > > ```r > A <- matrix(c(1,1,-1,-1), 2, 2); b <- c(-1, 1) > det(t(A) %*% A) # 0 > z <- MASS::ginv(A) %*% b # pseudoinverse -> min-norm > c(z, sqrt(sum((A %*% z - b)^2))) > ``` > *Look for:* determinant $0$; `z` $\approx(0,0)$; residual $\approx\sqrt2$. > [!example]- MATLAB > > ```matlab > A = [1 -1; 1 -1]; b = [-1; 1]; > det(A'*A) % 0 > z = pinv(A)*b; % min-norm (backslash warns on singular) > [z.' , norm(A*z - b)] > ``` > *Look for:* determinant $0$; `z` $\approx[0\ 0]$; residual $\approx\sqrt2$. > [!example]- Mathematica *(numerical only)* > > ```wolfram > A = {{1., -1.}, {1., -1.}}; b = {-1., 1.}; > Det[Transpose[A].A] (* 0 *) > z = PseudoInverse[A].b; (* min-norm *) > {z, Norm[A.z - b]} > ``` > *Look for:* determinant $0$; `z` $\approx\{0,0\}$; residual $\approx\sqrt2$. --- ## Problem 2 — Steam-table Lagrange interpolation **What to hand in.** - **By hand (or symbolically):** the three Lagrange basis polynomials for the nodes $180,200,220$, and the value at $210$. - **In code:** build $L(x)$ and plot it with the data over $[175,225]$. - **Report:** $h(210)$, the plot, and the reflection on why Lagrange is unsafe for a large table. > [!hint]- Hint 1 — the Lagrange form > With nodes $x_0,x_1,x_2=180,200,220$ and values $h_0,h_1,h_2=763,853,945$, $L(x)=\sum_i h_i\,\ell_i(x)$ where $\ell_i(x)=\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}$. Each $\ell_i$ is $1$ at its own node and $0$ at the others. > [!hint]- Hint 2 — evaluate at 210 > Plug $x=210$ into the three $\ell_i$ and combine. You should get $h(210)=898.75$ kJ/kg. Since $210$ sits between tabulated points and the data is gently curved, this is trustworthy. > [!hint]- Hint 3 — plotting > Evaluate $L$ on a fine grid across $[175,225]$ and overlay the three data points as markers. The quadratic through three points is smooth and monotone here — reasonable. > [!hint]- Hint 4 — why not a whole steam table > A single polynomial through $n{+}1$ points has degree $n$; over many nodes it oscillates wildly between them ([Runge's phenomenon](https://en.wikipedia.org/wiki/Runge%27s_phenomenon)). Practice uses low-degree **piecewise** interpolation instead — e.g. cubic splines (Problem 3). > [!example]- Python > > ```python > import numpy as np > xs = np.array([180., 200., 220.]); hs = np.array([763., 853., 945.]) > def L(x): > s = 0.0 > for i in range(3): > li = 1.0 > for j in range(3): > if j != i: li *= (x - xs[j]) / (xs[i] - xs[j]) # TODO: basis product > s += hs[i] * li > return s > print(L(210.0)) # 898.75 > ``` > *Look for:* `L(210)` = 898.75; `L(xs)` returns the data back exactly. > [!example]- R > > ```r > xs <- c(180,200,220); hs <- c(763,853,945) > L <- function(x) { s <- 0 > for (i in 1:3) { li <- 1 > for (j in 1:3) if (j != i) li <- li * (x - xs[j])/(xs[i]-xs[j]) > s <- s + hs[i]*li } > s } > L(210) # 898.75 > ``` > [!example]- MATLAB > > ```matlab > xs = [180 200 220]; hs = [763 853 945]; > L = @(x) 0; % build the sum of basis terms: > for i = 1:3 > li = @(x) 1; > for j = 1:3 > if j ~= i, li = @(x) li(x).*(x-xs(j))/(xs(i)-xs(j)); end > end > L = @(x) L(x) + hs(i)*li(x); > end > L(210) % 898.75 > ``` > [!example]- Mathematica *(numerical only)* > > ```wolfram > xs = {180., 200., 220.}; hs = {763., 853., 945.}; > L[x_] := Sum[hs[[i]] Product[If[j != i, (x - xs[[j]])/(xs[[i]] - xs[[j]]), 1], {j, 3}], {i, 3}]; > L[210] (* 898.75 *) > ``` --- ## Problem 3 — Cubic spline end conditions on a closed curve **What to hand in.** - **By hand:** identify which rows of the spline's linear system carry the end conditions, and what "natural" vs. "periodic" put there. - **In code:** spline $x(t)$ and $y(t)$ separately against $t$ with **natural** and **periodic** end conditions; plot both over the true curve; measure the seam corner and the max error; refine $n$. - **Report:** the seam plot, the corner-angle and max-distance numbers, and the refinement table. > [!hint]- Hint 1 — parametrize, then spline each coordinate > Sample $t$ at $n{+}1$ equally spaced values on $[0,2\pi]$ (the last repeats the first) and spline $x(t)$ and $y(t)$ *against $t$*, not against each other. The class file `CubicClosedCurveTest.m` is the exact template — start from it. > [!hint]- Hint 2 — what the two end conditions do at the seam > The **natural** spline sets the second derivative to zero at $t=0$ and $t=2\pi$ and treats them as two unrelated free ends — it has no way to know the curve closes, so it leaves a slope/curvature mismatch (a visible corner) at the seam. The **periodic** spline forces value, slope, and curvature to match across $t=0\equiv2\pi$, closing it smoothly. > [!hint]- Hint 3 — measuring the corner > Read the tangent heading just after the seam, $\operatorname{atan2}(y'(0^+),x'(0^+))$, and just before it, $\operatorname{atan2}(y'(2\pi^-),x'(2\pi^-))$; their difference (reduced to $(-180^\circ,180^\circ]$) is the corner angle. Periodic $\approx0$; natural leaves a real kink. For the fit error, take the max of $\sqrt{(x-x_{\text{true}})^2+(y-y_{\text{true}})^2}$ over a dense $t$. > [!hint]- Hint 4 — refine and connect to the derivation > Tabulate, for $n=7,11,17,25,41$, the max distance from the true curve *and* the seam corner angle. Both splines' max-distance error falls as you add points, but the **seam corner** is the qualitative tell: periodic holds it at $\approx0^\circ$ throughout, while the natural spline carries a real corner (tens of degrees at $n=7$) that only slowly shrinks. That corner is the fingerprint of the two boundary rows — swap natural for periodic and the closure changes completely. > [!example]- Python *(built-in spline carries the end condition)* > > ```python > import numpy as np > from scipy.interpolate import CubicSpline > xc = lambda t: 5*np.cos(t) + 0.6*np.cos(2*t) - 0.25*np.cos(5*t) > yc = lambda t: 3*np.sin(t) + 0.35*np.sin(3*t) + 0.15*np.sin(6*t) > n = 7; tk = np.linspace(0, 2*np.pi, n+1); td = np.linspace(0, 2*np.pi, 2000) > xk, yk = xc(tk), yc(tk) > xN = CubicSpline(tk, xk, bc_type="natural"); yN = CubicSpline(tk, yk, bc_type="natural") > xkp, ykp = xk.copy(), yk.copy(); xkp[-1] = xkp[0]; ykp[-1] = ykp[0] # periodic wants EXACT closure > xP = CubicSpline(tk, xkp, bc_type="periodic"); yP = CubicSpline(tk, ykp, bc_type="periodic") > errP = np.max(np.hypot(xP(td) - xc(td), yP(td) - yc(td))) # TODO: same for natural > ``` > *Look for:* the periodic curve closes at the seam; the natural one leaves a corner. **Gotcha:** the periodic solver rejects the data unless the last sampled value *exactly* equals the first — but $x(2\pi)$ differs from $x(0)$ by rounding, so copy the first value into the last (as above). R's `method="periodic"` wants the same. > [!example]- R *(built-in spline carries the end condition)* > > ```r > xc <- function(t) 5*cos(t)+0.6*cos(2*t)-0.25*cos(5*t) > yc <- function(t) 3*sin(t)+0.35*sin(3*t)+0.15*sin(6*t) > n <- 7; tk <- seq(0,2*pi,length.out=n+1); td <- seq(0,2*pi,length.out=2000) > xN <- splinefun(tk, xc(tk), method="natural"); xP <- splinefun(tk, xc(tk), method="periodic") > # same for yN,yP; evaluate xN(td), xP(td), ... ; compare to xc(td),yc(td) > ``` > [!example]- MATLAB *(your class code)* > > ```matlab > xc = @(t) 5*cos(t)+0.6*cos(2*t)-0.25*cos(5*t); > yc = @(t) 3*sin(t)+0.35*sin(3*t)+0.15*sin(6*t); > n = 7; tk = linspace(0,2*pi,n+1); td = linspace(0,2*pi,2000); > [xP,xPp] = cubicInterpPeriodic(tk, xc(tk), td); % periodic > [xN,xNp] = cubicInterp(tk, xc(tk), td); % natural > % same for y; CubicClosedCurveTest.m has the full driver, seam zoom, and corner/error prints > ``` > [!example]- Mathematica *(numerical only)* > > ```wolfram > xc[t_] := 5 Cos[t] + 0.6 Cos[2 t] - 0.25 Cos[5 t]; > yc[t_] := 3 Sin[t] + 0.35 Sin[3 t] + 0.15 Sin[6 t]; > n = 7; tk = Subdivide[0., 2 Pi, n]; > xP = Interpolation[Transpose[{tk, xc /@ tk}], PeriodicInterpolation -> True]; (* periodic *) > xN = Interpolation[Transpose[{tk, xc /@ tk}], InterpolationOrder -> 3]; (* non-periodic cubic *) > ``` > *Look for:* the periodic interpolant closes the seam; the plain cubic leaves it open. --- ## Problem 4 — Newton's method for systems (critical points) **What to hand in.** - **By hand:** the classification rule (Hessian eigenvalue signs $\to$ min / max / saddle). - **In code:** a surface/contour plot; Newton for systems with the gradient and Hessian from finite differences; two critical points and their types. - **Report:** the plot, the two critical points you converged to, and each one's classification. > [!hint]- Hint 1 — the gradient and Hessian from samples of $f$ only > Centered gradient $f_x\approx\frac{f(x+h,y)-f(x-h,y)}{2h}$, $f_y$ likewise; central second differences for $f_{xx},f_{yy}$; the four-corner mixed difference for $f_{xy}$. Assemble $H=\left[\begin{smallmatrix}f_{xx}&f_{xy}\\f_{xy}&f_{yy}\end{smallmatrix}\right]$. The class file `NewtonGrid_CriticalPoints2D_FD.m` (and its `.py/.R/.nb` ports) has `gradFD`/`hessFD` ready. > [!hint]- Hint 2 — the Newton-for-systems step > The scalar $x\leftarrow x-f'/f''$ becomes $\vec p\leftarrow\vec p-H^{-1}\nabla f$: solve $H\,\vec s=\nabla f$, then $\vec p\leftarrow\vec p-\vec s$. Iterate until $\lVert\nabla f\rVert$ is tiny. Geometrically you are fitting a paraboloid to the samples and jumping to its vertex. > [!hint]- Hint 3 — where to start > The plot shows four critical points near $(\pm1,\pm1)$. Start each Newton run *near* the one you want — e.g. $(0.8,0.8)$ lands on $(1,1)$, $(-0.8,0.8)$ on $(-1,1)$. Two well-chosen starts give you the two the problem asks for. > [!hint]- Hint 4 — classify from the eigenvalues > Evaluate $H$ at the landing point and take its eigenvalues: both positive $\Rightarrow$ local min, both negative $\Rightarrow$ local max, opposite signs $\Rightarrow$ saddle. (For this $f$ the Hessian is diagonal at each critical point, so the signs are easy to read.) > [!example]- Python > > ```python > import numpy as np > f = lambda x, y: np.exp(-(x**3/3 - x) - (y**3/3 - y)) > def gradFD(p, h=1e-4): > x, y = p > return np.array([(f(x+h,y)-f(x-h,y))/(2*h), (f(x,y+h)-f(x,y-h))/(2*h)]) > def hessFD(p, h=1e-4): > x, y = p > fxx=(f(x+h,y)-2*f(x,y)+f(x-h,y))/h**2; fyy=(f(x,y+h)-2*f(x,y)+f(x,y-h))/h**2 > fxy=(f(x+h,y+h)-f(x+h,y-h)-f(x-h,y+h)+f(x-h,y-h))/(4*h**2) > return np.array([[fxx,fxy],[fxy,fyy]]) > def newton(p, tol=1e-10): > p = np.array(p, float) > for _ in range(100): > g = gradFD(p) > if np.linalg.norm(g) < tol: break > p = p - np.linalg.solve(hessFD(p), g) # TODO: the Newton-for-systems step > return p > p = newton([0.8, 0.8]); print(p, np.linalg.eigvalsh(hessFD(p))) # -> (1,1), eig<0 : max > ``` > *Look for:* $(0.8,0.8)\to(1,1)$ (both eigenvalues negative, a max); $(-0.8,-0.8)\to(-1,-1)$ (both positive, a min). > [!example]- R > > ```r > f <- function(x,y) exp(-(x^3/3 - x) - (y^3/3 - y)) > # gradFD / hessFD / newton as in NewtonGrid_CriticalPoints2D_FD.R (class code); then: > # p <- newton(c(0.8,0.8)); eigen(hessFD(p))$values > ``` > [!example]- MATLAB *(your class code)* > > ```matlab > f = @(x,y) exp(-(x.^3/3 - x) - (y.^3/3 - y)); > % gradFD, hessFD, newtonND live in NewtonGrid_CriticalPoints2D_FD.m > p = newtonND(f, [0.8; 0.8], 1e-4, 1e-10, 100); % -> (1,1) > eig(hessFD(f, p, 1e-4)) % both < 0 : local max > ``` > [!example]- Mathematica *(numerical only)* > > ```wolfram > f[x_, y_] := Exp[-(x^3/3 - x) - (y^3/3 - y)]; > (* gradFD, hessFD, newton in NewtonGrid_CriticalPoints2D_FD.nb (class code) *) > (* p = newton[{0.8, 0.8}]; Eigenvalues[hessFD[p]] -> both < 0, a max *) > ``` --- ## Problem 5 — Double integration: built-in vs. a nested loop **What to hand in.** - **By hand:** note that $e^{-x^2-y^2}=e^{-x^2}e^{-y^2}$ separates, so the exact value is the square of the 1-D Gaussian integral. - **In code:** a built-in adaptive 2-D integral, and a tensor-product composite Simpson you build yourself; compare accuracy and time. - **Report:** each method's error against $(\sqrt\pi\,\mathrm{erf}\,1)^2$, the timings, and the $n^2$-cost discussion. > [!hint]- Hint 1 — the built-in > Call your language's 2-D adaptive integrator on $f=e^{-x^2-y^2}$ over $[-1,1]^2$ and compare to the exact $(\sqrt\pi\,\mathrm{erf}\,1)^2=2.2309851$. This is your accuracy yardstick. > [!hint]- Hint 2 — nest a 1-D rule > Fix a grid of $y_j$. For each, integrate in $x$ with composite Simpson to get $g(y_j)=\int_{-1}^1 f(x,y_j)\,dx$. Then apply Simpson again to the $g(y_j)$ in $y$. That double pass is a tensor-product Simpson rule (this is exactly what `DoubleIntegral.m` does with `mySimp`). > [!hint]- Hint 3 — compare fairly > Time both (start with, say, $101\times101$). The built-in adapts its points; your loop uses a fixed $n\times n$ grid. Report both errors and both run times. > [!hint]- Hint 4 — the $n^2$ lesson > A grid of $n$ points per axis costs $n^2$ evaluations, and in $d$ dimensions $n^d$ (the [curse of dimensionality](https://en.wikipedia.org/wiki/Curse_of_dimensionality)). For low dimension and smooth integrands the loop is simple and fine; as dimension or accuracy demand grows, adaptive built-ins (and, higher up, Monte Carlo) win. > [!example]- Python > > ```python > import numpy as np, time > from scipy import integrate > from scipy.special import erf > f = lambda x, y: np.exp(-x**2 - y**2) > exact = (np.sqrt(np.pi)*erf(1))**2 > val, _ = integrate.dblquad(lambda y, x: f(x, y), -1, 1, -1, 1) # built-in > def simp(g, a, b, n): > x = np.linspace(a, b, n+1); y = g(x); h = (b-a)/n > return h/3*(y[0] + y[-1] + 4*y[1:-1:2].sum() + 2*y[2:-1:2].sum()) > n = 100; ys = np.linspace(-1, 1, n+1) > g = np.array([simp(lambda x: f(x, yj), -1, 1, n) for yj in ys]) # TODO: inner integral per y_j > loop = simp(lambda t: np.interp(t, ys, g), -1, 1, n) # or Simpson directly on g > print(abs(val-exact), abs(loop-exact)) > ``` > *Look for:* both errors small; the built-in reaches machine accuracy, the fixed grid a bit coarser. > [!example]- R > > ```r > f <- function(x,y) exp(-x^2 - y^2); exact <- (sqrt(pi)*erf(1))^2 # erf via 2*pnorm(x*sqrt2)-1 > # built-in: pracma::integral2(f, -1,1,-1,1)$Q > # nested: for each y_j, Simpson in x (your mySimp / Day-14 rule), then Simpson the results in y > ``` > [!example]- MATLAB *(your class code)* > > ```matlab > f = @(x,y) exp(-x.^2 - y.^2); exact = (sqrt(pi)*erf(1))^2; > tic; Q = integral2(f, -1,1,-1,1, 'AbsTol',1e-14,'RelTol',1e-14); toc % built-in > % nested: g(j) = mySimp(-1,1,@(xx) f(xx,y(j)), nx); then Simpson g in y (see DoubleIntegral.m) > ``` > *Look for:* `integral2` error near machine level; the $101\times101$ loop a few digits behind, and slower per digit. > [!example]- Mathematica *(numerical only)* > > ```wolfram > f[x_, y_] := Exp[-x^2 - y^2]; exact = (Sqrt[Pi] Erf[1])^2; > NIntegrate[f[x, y], {x, -1, 1}, {y, -1, 1}] (* built-in adaptive *) > (* nested: Table of NIntegrate/Simpson in x per y_j, then Simpson in y *) > ``` --- *Companion to `assignment/MATH307Su26-HW4.tex`. Code here drives the machinery you built in lecture and the class-code files (Lagrange, `cubicInterp`/`cubicInterpPeriodic`, FD-Newton, `mySimp`) plus a built-in integrator. The `# TODO` marks and the hand calculations are yours. Submit one written PDF plus your code files and figures.*