# HW1 Helper — Normal Distributions, Taylor Series & Finite Differences
> [!info]- How to use
> Read **what the problem asks** and **the mathematics**, then open the **Code Bases** callout for
> the language you are working in: each one walks the code one small piece at a time, with the
> explanation sitting right next to the lines it describes.
## Problem 1 — The normal density and the Gaussian integral
**What the problem actually asks.** Three things, in order:
1. **Prove the normalization** $\int_{\mathbb{R}}\rho(x;\mu,\sigma)\,dx=1$ for the
[normal distribution](https://en.wikipedia.org/wiki/Normal_distribution), by first establishing the
**[Gaussian integral](https://en.wikipedia.org/wiki/Gaussian_integral)**
$I=\int_{-\infty}^{\infty}e^{-x^2}\,dx=\sqrt{\pi}$ with the polar-coordinate
trick, then recovering the general case by shifting and rescaling.
2. **Define `normal_pdf`** in your code base with adjustable $\mu$ and $\sigma$.
3. **Make two list plots** of the density family and comment.
The proof is the heart of the problem; the code defines the function and visualizes it.
**The mathematics — the polar trick.** You cannot antidifferentiate $e^{-x^2}$ in elementary terms,
so instead of $I$ you work with $I^2$ and let a change to [polar coordinates](https://en.wikipedia.org/wiki/Polar_coordinate_system) do the work:
$I^2=\Big(\!\int_{\mathbb{R}}e^{-x^2}\,dx\Big)\Big(\!\int_{\mathbb{R}}e^{-y^2}\,dy\Big)
=\iint_{\mathbb{R}^2}e^{-(x^2+y^2)}\,dA
=\int_{0}^{2\pi}\!\!\int_{0}^{\infty}e^{-r^2}\,r\,dr\,d\theta .$
The one new ingredient is the [Jacobian](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant) factor $r$ from $dA=r\,dr\,d\theta$ — and it is exactly what
makes the radial integral elementary: with $u=r^2$, $\int_0^{\infty}r\,e^{-r^2}\,dr=\tfrac12$. The
angle contributes $2\pi$, so $I^2=2\pi\cdot\tfrac12=\pi$, and since $I>0$, $I=\sqrt{\pi}$. (Rigor
needs [Fubini's theorem](https://en.wikipedia.org/wiki/Fubini's_theorem) to turn the product into the double integral, a limit/[squeeze](https://en.wikipedia.org/wiki/Squeeze_theorem) argument for the improper
bounds, and the positivity of $I$.) Then $\rho(x;0,1/\sqrt2)=e^{-x^2}/\sqrt{\pi}$ gives
$\int\rho\,dx=I/\sqrt{\pi}=1$, and for general $\mu,\sigma$ a $u$-sub recenters and a $v$-sub rescales
back to the standard $e^{-x^2}$.
### Code Bases
> [!example]- Python
>
> **Step 1 — Define the density.** This is deliverable (2): a `normal_pdf` you can call with any mean and spread. It is written to evaluate on a scalar *or* a whole array of `x` at once (vectorized).
> ```python
> import numpy as np
> import matplotlib.pyplot as plt
>
> def normal_pdf(x, mu=0.0, sigma=1.0):
> return np.exp(-(x - mu)**2 / (2*sigma**2)) / (sigma*np.sqrt(2*np.pi))
> ```
> *Look for:* `normal_pdf(0,0,1)` returns $1/\sqrt{2\pi}\approx 0.399$.
>
> **Step 2 — Choose a sample grid and the two families.** A "list plot" is built from discrete samples, so we fix a grid of `x` and list the $(\mu,\sigma)$ pairs the problem gives: one family slides the mean, the other changes the spread.
> ```python
> xList = np.linspace(-6, 6, 49) # discrete samples for the list plots
> muFamily = [(-2, 1), (0, 1), (2, 1)] # fixed spread, shifting the mean
> sigmaFamily = [(0, 0.5), (0, 1), (0, 2)] # fixed mean, changing the spread
> ```
> *Look for:* `xList` is 49 evenly spaced samples on $[-6,6]$ — fine enough to trace the bell smoothly.
>
> **Step 3 — Plot the mean family (deliverable 3a).** Loop over the three $(\mu,\sigma)$ pairs with $\sigma=1$ and overlay them. Drawing the samples as points joined by line segments makes it a *list-line* plot (points only would be a plain list plot).
> ```python
> plt.figure()
> for mu, sigma in muFamily:
> plt.plot(xList, normal_pdf(xList, mu, sigma), marker="o", ms=3,
> label=f"$\\mu={mu},\\ \\sigma={sigma}
quot;)
> plt.xlabel("x"); plt.ylabel(r"$\rho(x;\mu,\sigma)quot;)
> plt.title("Normal density: shifting the mean")
> plt.legend(); plt.grid(True)
> plt.savefig("hw1_p1_pdf_mu.png", dpi=150, bbox_inches="tight") # saved to the cwd
> ```
> *Look for:* the same bell shape translated to be centered at $-2$, $0$, and $2$.
>
> **Step 4 — Plot the spread family (deliverable 3b).** Same idea with $\mu=0$ fixed and $\sigma\in\{\tfrac12,1,2\}$. This is the comment the problem wants: the area is always $1$, so a smaller $\sigma$ forces a taller, narrower peak.
> ```python
> plt.figure()
> for mu, sigma in sigmaFamily:
> plt.plot(xList, normal_pdf(xList, mu, sigma), marker="o", ms=3,
> label=f"$\\mu={mu},\\ \\sigma={sigma}quot;)
> plt.xlabel("x"); plt.ylabel(r"$\rho(x;\mu,\sigma)quot;)
> plt.title("Normal density: changing the spread")
> plt.legend(); plt.grid(True)
> plt.savefig("hw1_p1_pdf_sigma.png", dpi=150, bbox_inches="tight") # saved to the cwd
> plt.show() # display both figures
> ```
> *Look for:* peak height $=1/(\sigma\sqrt{2\pi})$ — about $0.80,\,0.40,\,0.20$ for $\sigma=\tfrac12,1,2$.
>
> **Step 5 (optional, peeks ahead) — confirm the normalization numerically.** The proof already settles $\int\rho=1$; here we just *watch* it happen with a trapezoid sum. Numerical integration is a later topic, so treat this as a convenience check, not part of the main thread.
> ```python
> def trapz(y, x):
> return np.sum((y[:-1] + y[1:]) / 2 * np.diff(x))
>
> xx = np.linspace(-20, 20, 200001)
> print(trapz(np.exp(-xx**2), xx)) # ~ sqrt(pi) = 1.7724539
> print(trapz(normal_pdf(xx, 0.0, 1/np.sqrt(2)), xx)) # ~ 1
> ```
> *Look for:* $I\approx 1.7724539$ and area $\approx 1.000$.
>
> [!example]- MATLAB
>
> **Step 1 — Define the density.** This is deliverable (2): a `normal_pdf` you can call with any mean and spread. It is written to evaluate on a scalar *or* a whole array of `x` at once (vectorized).
> ```matlab
> normal_pdf = @(x, mu, sigma) ...
> exp(-(x - mu).^2 ./ (2*sigma.^2)) ./ (sigma .* sqrt(2*pi));
> ```
> *Look for:* `normal_pdf(0,0,1)` returns $1/\sqrt{2\pi}\approx 0.399$.
>
> **Step 2 — Choose a sample grid and the two families.** A "list plot" is built from discrete samples, so we fix a grid of `x` and list the $(\mu,\sigma)$ pairs the problem gives: one family slides the mean, the other changes the spread.
> ```matlab
> xList = linspace(-6, 6, 49);
> muFamily = [-2 1; 0 1; 2 1]; % fixed spread, shifting the mean
> sigmaFamily = [0 0.5; 0 1; 0 2]; % fixed mean, changing the spread
> ```
> *Look for:* `xList` is 49 evenly spaced samples on $[-6,6]$ — fine enough to trace the bell smoothly.
>
> **Step 3 — Plot the mean family (deliverable 3a).** Loop over the three $(\mu,\sigma)$ pairs with $\sigma=1$ and overlay them. Drawing the samples as points joined by line segments makes it a *list-line* plot (points only would be a plain list plot).
> ```matlab
> figure; hold on
> for r = 1:size(muFamily,1)
> plot(xList, normal_pdf(xList, muFamily(r,1), muFamily(r,2)), '-o', 'MarkerSize', 3, ...
> 'DisplayName', sprintf('\\mu=%g, \\sigma=%g', muFamily(r,1), muFamily(r,2)));
> end
> xlabel('x'); ylabel('\rho(x;\mu,\sigma)');
> title('Normal density: shifting the mean'); legend show; grid on
> exportgraphics(gcf, 'hw1_p1_pdf_mu.png', 'Resolution', 150); % figure shows on screen; PNG saved to pwd
> hold off
> ```
> *Look for:* the same bell shape translated to be centered at $-2$, $0$, and $2$.
>
> **Step 4 — Plot the spread family (deliverable 3b).** Same idea with $\mu=0$ fixed and $\sigma\in\{\tfrac12,1,2\}$. This is the comment the problem wants: the area is always $1$, so a smaller $\sigma$ forces a taller, narrower peak.
> ```matlab
> figure; hold on
> for r = 1:size(sigmaFamily,1)
> plot(xList, normal_pdf(xList, sigmaFamily(r,1), sigmaFamily(r,2)), '-o', 'MarkerSize', 3, ...
> 'DisplayName', sprintf('\\mu=%g, \\sigma=%g', sigmaFamily(r,1), sigmaFamily(r,2)));
> end
> xlabel('x'); ylabel('\rho(x;\mu,\sigma)');
> title('Normal density: changing the spread'); legend show; grid on
> exportgraphics(gcf, 'hw1_p1_pdf_sigma.png', 'Resolution', 150); % figure shows on screen; PNG saved to pwd
> hold off
> ```
> *Look for:* peak height $=1/(\sigma\sqrt{2\pi})$ — about $0.80,\,0.40,\,0.20$ for $\sigma=\tfrac12,1,2$.
>
> **Step 5 (optional, peeks ahead) — confirm the normalization numerically.** The proof already settles $\int\rho=1$; here we just *watch* it happen with a trapezoid sum. Numerical integration is a later topic, so treat this as a convenience check, not part of the main thread.
> ```matlab
> trapz_ = @(y, x) sum((y(1:end-1) + y(2:end))/2 .* diff(x));
> xx = linspace(-20, 20, 200001);
> fprintf('I = %.12f (sqrt(pi) = %.12f)\n', trapz_(exp(-xx.^2), xx), sqrt(pi));
> fprintf('int rho = %.12f (should be 1)\n', trapz_(normal_pdf(xx, 0, 1/sqrt(2)), xx));
> ```
> *Look for:* $I\approx 1.7724539$ and area $\approx 1.000$.
>
> [!example]- R
>
> **Step 1 — Define the density.** This is deliverable (2): a `normal_pdf` you can call with any mean and spread. It is written to evaluate on a scalar *or* a whole array of `x` at once (vectorized).
> ```r
> normal_pdf <- function(x, mu = 0, sigma = 1)
> exp(-(x - mu)^2 / (2 * sigma^2)) / (sigma * sqrt(2 * pi))
> ```
> *Look for:* `normal_pdf(0,0,1)` returns $1/\sqrt{2\pi}\approx 0.399$.
>
> **Step 2 — Choose a sample grid and the two families.** A "list plot" is built from discrete samples, so we fix a grid of `x` and list the $(\mu,\sigma)$ pairs the problem gives: one family slides the mean, the other changes the spread.
> ```r
> xList <- seq(-6, 6, length.out = 49)
> muFamily <- list(c(-2, 1), c(0, 1), c(2, 1)) # shifting the mean
> sigmaFamily <- list(c(0, 0.5), c(0, 1), c(0, 2)) # changing the spread
> ```
> *Look for:* `xList` is 49 evenly spaced samples on $[-6,6]$ — fine enough to trace the bell smoothly.
>
> **Step 3 — Plot the mean family (deliverable 3a).** Loop over the three $(\mu,\sigma)$ pairs with $\sigma=1$ and overlay them. Drawing the samples as points joined by line segments makes it a *list-line* plot (points only would be a plain list plot).
> ```r
> cols <- c("black", "red", "blue")
> draw_mu <- function() {
> plot(NA, xlim = range(xList), ylim = c(0, 0.45),
> xlab = "x", ylab = expression(rho(x)), main = "Normal density: shifting the mean")
> for (k in seq_along(muFamily)) {
> p <- muFamily[[k]]
> lines(xList, normal_pdf(xList, p[1], p[2]), type = "o", pch = 16, cex = 0.5, col = cols[k])
> }
> legend("topright", col = cols, lty = 1, pch = 16,
> legend = sapply(muFamily, function(p) sprintf("mu=%g, sigma=%g", p[1], p[2])))
> }
> if (interactive()) draw_mu() # show on screen (RStudio)
> png("hw1_p1_pdf_mu.png", width = 800, height = 600); draw_mu(); dev.off() # save to getwd()
> ```
> *Look for:* the same bell shape translated to be centered at $-2$, $0$, and $2$.
>
> **Step 4 — Plot the spread family (deliverable 3b).** Same idea with $\mu=0$ fixed and $\sigma\in\{\tfrac12,1,2\}$. This is the comment the problem wants: the area is always $1$, so a smaller $\sigma$ forces a taller, narrower peak.
> ```r
> draw_sigma <- function() {
> plot(NA, xlim = range(xList), ylim = c(0, 0.85),
> xlab = "x", ylab = expression(rho(x)), main = "Normal density: changing the spread")
> for (k in seq_along(sigmaFamily)) {
> p <- sigmaFamily[[k]]
> lines(xList, normal_pdf(xList, p[1], p[2]), type = "o", pch = 16, cex = 0.5, col = cols[k])
> }
> legend("topright", col = cols, lty = 1, pch = 16,
> legend = sapply(sigmaFamily, function(p) sprintf("mu=%g, sigma=%g", p[1], p[2])))
> }
> if (interactive()) draw_sigma() # show on screen (RStudio)
> png("hw1_p1_pdf_sigma.png", width = 800, height = 600); draw_sigma(); dev.off() # save to getwd()
> ```
> *Look for:* peak height $=1/(\sigma\sqrt{2\pi})$ — about $0.80,\,0.40,\,0.20$ for $\sigma=\tfrac12,1,2$.
>
> **Step 5 (optional, peeks ahead) — confirm the normalization numerically.** The proof already settles $\int\rho=1$; here we just *watch* it happen with a trapezoid sum. Numerical integration is a later topic, so treat this as a convenience check, not part of the main thread.
> ```r
> trapz <- function(y, x) sum((head(y, -1) + tail(y, -1)) / 2 * diff(x))
> xx <- seq(-20, 20, length.out = 200001)
> cat(sprintf("I = %.12f (sqrt(pi) = %.12f)\n", trapz(exp(-xx^2), xx), sqrt(pi)))
> cat(sprintf("int rho = %.12f (should be 1)\n", trapz(normal_pdf(xx, 0, 1/sqrt(2)), xx)))
> ```
> *Look for:* $I\approx 1.7724539$ and area $\approx 1.000$.
>
> [!example]- Mathematica
>
> **Step 1 — Define the density.** This is deliverable (2): a `normal_pdf` you can call with any mean and spread. It is written to evaluate on a scalar *or* a whole array of `x` at once (vectorized).
> ```wolfram
> normalPDF[x_, mu_: 0, sigma_: 1] := Exp[-(x - mu)^2/(2 sigma^2)]/(sigma Sqrt[2 Pi]);
> ```
> *Look for:* `normal_pdf(0,0,1)` returns $1/\sqrt{2\pi}\approx 0.399$.
>
> **Step 2 — Choose a sample grid and the two families.** A "list plot" is built from discrete samples, so we fix a grid of `x` and list the $(\mu,\sigma)$ pairs the problem gives: one family slides the mean, the other changes the spread.
> ```wolfram
> xList = Subdivide[-6., 6., 48]; (* 49 samples *)
> muFamily = {{-2, 1}, {0, 1}, {2, 1}};
> sigmaFamily = {{0, 1/2}, {0, 1}, {0, 2}};
> ```
> *Look for:* `xList` is 49 evenly spaced samples on $[-6,6]$ — fine enough to trace the bell smoothly.
>
> **Step 3 — Plot the mean family (deliverable 3a).** Loop over the three $(\mu,\sigma)$ pairs with $\sigma=1$ and overlay them. Drawing the samples as points joined by line segments makes it a *list-line* plot (points only would be a plain list plot).
> ```wolfram
> muData = Table[Table[{x, normalPDF[x, p[[1]], p[[2]]]}, {x, xList}], {p, muFamily}];
> pMu = ListLinePlot[muData, Mesh -> All, AxesLabel -> {"x", "\[Rho]"},
> PlotLabel -> "Normal density: shifting the mean",
> PlotLegends -> (Row[{"\[Mu]=", #[[1]], ", \[Sigma]=", #[[2]]}] & /@ muFamily)];
> Export["hw1_p1_pdf_mu.png", pMu]; (* save to Directory[]; run SetDirectory[NotebookDirectory[]] first to use this notebook's folder *)
> pMu (* and display it in the notebook *)
> ```
> *Look for:* the same bell shape translated to be centered at $-2$, $0$, and $2$.
>
> **Step 4 — Plot the spread family (deliverable 3b).** Same idea with $\mu=0$ fixed and $\sigma\in\{\tfrac12,1,2\}$. This is the comment the problem wants: the area is always $1$, so a smaller $\sigma$ forces a taller, narrower peak.
> ```wolfram
> sigmaData = Table[Table[{x, normalPDF[x, p[[1]], p[[2]]]}, {x, xList}], {p, sigmaFamily}];
> pSigma = ListLinePlot[sigmaData, Mesh -> All, AxesLabel -> {"x", "\[Rho]"},
> PlotLabel -> "Normal density: changing the spread",
> PlotLegends -> (Row[{"\[Mu]=", #[[1]], ", \[Sigma]=", #[[2]]}] & /@ sigmaFamily)];
> Export["hw1_p1_pdf_sigma.png", pSigma]; (* save to Directory[]; SetDirectory[NotebookDirectory[]] to use this notebook's folder *)
> pSigma (* and display it in the notebook *)
> ```
> *Look for:* peak height $=1/(\sigma\sqrt{2\pi})$ — about $0.80,\,0.40,\,0.20$ for $\sigma=\tfrac12,1,2$.
>
> **Step 5 (optional, peeks ahead) — confirm the normalization numerically.** The proof already settles $\int\rho=1$; here we just *watch* it happen with a trapezoid sum. Numerical integration is a later topic, so treat this as a convenience check, not part of the main thread.
> ```wolfram
> trapz[y_, x_] := Total[(Most[y] + Rest[y])/2 Differences[x]];
> xx = Subdivide[-20., 20., 200000];
> Print["I = ", trapz[Exp[-xx^2], xx], " (Sqrt[Pi] = ", N[Sqrt[Pi]], ")"];
> Print["int rho = ", trapz[normalPDF[xx, 0, 1/Sqrt[2]], xx], " (should be 1)"];
> ```
> *Look for:* $I\approx 1.7724539$ and area $\approx 1.000$.
>
> [!check]- Checks (analytic truth)
> - $I=\int_{\mathbb{R}}e^{-x^2}\,dx=\sqrt{\pi}\approx 1.7724539$
> - $\int_{\mathbb{R}}\rho(x;0,1/\sqrt2)\,dx = 1$
> - peak height $=1/(\sigma\sqrt{2\pi})$; the area stays $1$ for every $\mu,\sigma$
*The two list plots the problem asks for:*
![[hw1_p1_pdf_mu.png]]
![[hw1_p1_pdf_sigma.png]]
## Problem 2 — First and second moments
**What the problem actually asks.** A pen-and-paper derivation, no code required: for the
standardized density ($\mu=0,\ \sigma=1/\sqrt2$, so $\rho=e^{-x^2}/\sqrt\pi$), show the **first
moment** $\mathbb{E}[X]$ by a
[$u$-substitution](https://en.wikipedia.org/wiki/Integration_by_substitution), and the
**second moment** $V[X]$ by writing $x^2=x\cdot x$ and
[integrating by parts](https://en.wikipedia.org/wiki/Integration_by_parts).
**The mathematics.**
*First moment.* By definition $\mathbb{E}[X]=\int_{\mathbb{R}} x\,\rho\,dx
=\tfrac1{\sqrt\pi}\int_{\mathbb{R}} x\,e^{-x^2}\,dx$. The integrand is odd, so it must vanish; to see
it mechanically, substitute $u=x^2$ ($du=2x\,dx$) to get the antiderivative
$\int x\,e^{-x^2}\,dx=-\tfrac12 e^{-x^2}$, which takes the same value ($0$) at both $\pm\infty$.
Hence
$\mathbb{E}[X]=0,$
i.e. the distribution sits on its mean $\mu=0$, as it should.
*Second moment / variance.* With $\mu=0$, $V[X]=\mathbb{E}[X^2]=\tfrac1{\sqrt\pi}\int_{\mathbb{R}}
x^2 e^{-x^2}\,dx$. Split $x^2=x\cdot x$ and integrate by parts with
$u=x,\qquad dv=x\,e^{-x^2}\,dx\ \Longrightarrow\ du=dx,\qquad v=-\tfrac12 e^{-x^2}.$
The point of the split is that $dv$ is *exactly the first-moment integrand* from above, whose
antiderivative we already have. Then
$\int_{\mathbb{R}} x^2 e^{-x^2}\,dx
=\underbrace{\Big[-\tfrac12 x\,e^{-x^2}\Big]_{-\infty}^{\infty}}_{=\,0}
+\tfrac12\int_{\mathbb{R}} e^{-x^2}\,dx
=\tfrac12\,I=\tfrac{\sqrt\pi}{2},$
where the boundary term dies because $x e^{-x^2}\to0$ at $\pm\infty$, and the leftover integral is
the **same $I=\sqrt\pi$ from Problem 1**. Therefore
$V[X]=\tfrac1{\sqrt\pi}\cdot\tfrac{\sqrt\pi}{2}=\tfrac12=\sigma^2 .$
No new integral was needed — only $I=\sqrt\pi$ reused. That is the real lesson here: this same
integration by parts turns the 4th moment into the 2nd, the 6th into the 4th, and so on, a
**moment recursion** for the Gaussian.
**Physical reading.** These are the same quantities from physics, applied to a probability density
instead of a mass density. The first moment is a [center of mass](https://en.wikipedia.org/wiki/Center_of_mass):
$\mathbb{E}[X]$ is the balance point of the density, exactly as $\bar x=\tfrac1M\int x\,dm$ balances
a rod. The second central moment is the [moment of inertia](https://en.wikipedia.org/wiki/Moment_of_inertia)
about that balance point — the variance measures how far the "mass" spreads from the center, just as
inertia measures resistance to spinning. In statistical language these are the
[expected value](https://en.wikipedia.org/wiki/Expected_value) and the
[variance](https://en.wikipedia.org/wiki/Variance); the umbrella idea (zeroth moment $=$ total mass,
first $=$ center of mass / mean, second $=$ inertia / variance) is the theory of
[moments](https://en.wikipedia.org/wiki/Moment_(mathematics)).
### Code Bases
Nothing here is required — the work above is the answer. But Problem 1's **optional** step already
handed us a quadrature routine (the `trapz` integrator), so we *might as well* reuse it to confirm
both moments numerically. Like that step, this is reassurance rather than a task, and it leans on
material that is itself a peek ahead.
> [!example]- Python
>
> **Step 1 — Set up the standardized density.** Reuse `normal_pdf` from Problem 1 and fix the standardized parameters $\mu=0,\ \sigma=1/\sqrt2$, so $\rho=e^{-x^2}/\sqrt\pi$.
> ```python
> import numpy as np
>
> def normal_pdf(x, mu=0.0, sigma=1.0):
> return np.exp(-(x - mu)**2 / (2*sigma**2)) / (sigma*np.sqrt(2*np.pi))
>
> mu, sigma = 0.0, 1/np.sqrt(2) # the standardized case
> ```
> *Look for:* this is the same density as Problem 1 with $\sigma=1/\sqrt2$.
>
> **Step 2 — Reuse the trapezoid integrator on a fine grid.** Moments are integrals, so we lay down a dense grid and the same `trapz` from Problem 1.
> ```python
> def trapz(y, x):
> return np.sum((y[:-1] + y[1:]) / 2 * np.diff(x))
>
> x = np.linspace(-20, 20, 400001)
> rho = normal_pdf(x, mu, sigma)
> ```
> *Look for:* a grid dense enough that the moment integrals are accurate.
>
> **Step 3 — First moment.** $\mathbb{E}[X]=\int x\,\rho\,dx$. The integrand $x\rho$ is odd, so this should come out to zero.
> ```python
> EX = trapz(x*rho, x)
> print(EX) # ~ 0
> ```
> *Look for:* a value near $0$ (about $10^{-17}$, i.e. roundoff) — the numerical echo of the odd-integrand argument.
>
> **Step 4 — Second moment (variance).** $V[X]=\int (x-\mu)^2\rho\,dx$ with $\mu=0$. The integration-by-parts derivation predicts $\tfrac12$.
> ```python
> VX = trapz((x - mu)**2 * rho, x)
> print(VX) # ~ 0.5
> ```
> *Look for:* $0.5=\sigma^2$, matching the hand calculation.
>
> [!example]- MATLAB
>
> **Step 1 — Set up the standardized density.** Reuse `normal_pdf` from Problem 1 and fix the standardized parameters $\mu=0,\ \sigma=1/\sqrt2$, so $\rho=e^{-x^2}/\sqrt\pi$.
> ```matlab
> normal_pdf = @(x, mu, sigma) ...
> exp(-(x - mu).^2 ./ (2*sigma.^2)) ./ (sigma .* sqrt(2*pi));
> mu = 0; sigma = 1/sqrt(2); % the standardized case
> ```
> *Look for:* this is the same density as Problem 1 with $\sigma=1/\sqrt2$.
>
> **Step 2 — Reuse the trapezoid integrator on a fine grid.** Moments are integrals, so we lay down a dense grid and the same `trapz` from Problem 1.
> ```matlab
> trapz_ = @(y, x) sum((y(1:end-1) + y(2:end))/2 .* diff(x));
> x = linspace(-20, 20, 400001);
> rho = normal_pdf(x, mu, sigma);
> ```
> *Look for:* a grid dense enough that the moment integrals are accurate.
>
> **Step 3 — First moment.** $\mathbb{E}[X]=\int x\,\rho\,dx$. The integrand $x\rho$ is odd, so this should come out to zero.
> ```matlab
> EX = trapz_(x .* rho, x); % first moment ~ 0 (odd integrand)
> fprintf('E[X] = %.3e (exact 0)\n', EX);
> ```
> *Look for:* a value near $0$ (about $10^{-17}$, i.e. roundoff) — the numerical echo of the odd-integrand argument.
>
> **Step 4 — Second moment (variance).** $V[X]=\int (x-\mu)^2\rho\,dx$ with $\mu=0$. The integration-by-parts derivation predicts $\tfrac12$.
> ```matlab
> VX = trapz_((x - mu).^2 .* rho, x); % second moment ~ 1/2
> fprintf('V[X] = %.12f (exact sigma^2 = %.12f)\n', VX, sigma^2);
> ```
> *Look for:* $0.5=\sigma^2$, matching the hand calculation.
>
> [!example]- R
>
> **Step 1 — Set up the standardized density.** Reuse `normal_pdf` from Problem 1 and fix the standardized parameters $\mu=0,\ \sigma=1/\sqrt2$, so $\rho=e^{-x^2}/\sqrt\pi$.
> ```r
> normal_pdf <- function(x, mu = 0, sigma = 1)
> exp(-(x - mu)^2 / (2 * sigma^2)) / (sigma * sqrt(2 * pi))
> mu <- 0; sigma <- 1/sqrt(2) # the standardized case
> ```
> *Look for:* this is the same density as Problem 1 with $\sigma=1/\sqrt2$.
>
> **Step 2 — Reuse the trapezoid integrator on a fine grid.** Moments are integrals, so we lay down a dense grid and the same `trapz` from Problem 1.
> ```r
> trapz <- function(y, x) sum((head(y, -1) + tail(y, -1)) / 2 * diff(x))
> x <- seq(-20, 20, length.out = 400001)
> rho <- normal_pdf(x, mu, sigma)
> ```
> *Look for:* a grid dense enough that the moment integrals are accurate.
>
> **Step 3 — First moment.** $\mathbb{E}[X]=\int x\,\rho\,dx$. The integrand $x\rho$ is odd, so this should come out to zero.
> ```r
> EX <- trapz(x * rho, x) # first moment ~ 0 (odd integrand)
> cat(sprintf("E[X] = %.3e (exact 0)\n", EX))
> ```
> *Look for:* a value near $0$ (about $10^{-17}$, i.e. roundoff) — the numerical echo of the odd-integrand argument.
>
> **Step 4 — Second moment (variance).** $V[X]=\int (x-\mu)^2\rho\,dx$ with $\mu=0$. The integration-by-parts derivation predicts $\tfrac12$.
> ```r
> VX <- trapz((x - mu)^2 * rho, x) # second moment ~ 1/2
> cat(sprintf("V[X] = %.12f (exact sigma^2 = %.12f)\n", VX, sigma^2))
> ```
> *Look for:* $0.5=\sigma^2$, matching the hand calculation.
>
> [!example]- Mathematica
>
> **Step 1 — Set up the standardized density.** Reuse `normal_pdf` from Problem 1 and fix the standardized parameters $\mu=0,\ \sigma=1/\sqrt2$, so $\rho=e^{-x^2}/\sqrt\pi$.
> ```wolfram
> normalPDF[x_, mu_: 0, sigma_: 1] := Exp[-(x - mu)^2/(2 sigma^2)]/(sigma Sqrt[2 Pi]);
> mu = 0; sigma = 1/Sqrt[2]; (* the standardized case *)
> ```
> *Look for:* this is the same density as Problem 1 with $\sigma=1/\sqrt2$.
>
> **Step 2 — Reuse the trapezoid integrator on a fine grid.** Moments are integrals, so we lay down a dense grid and the same `trapz` from Problem 1.
> ```wolfram
> trapz[y_, x_] := Total[(Most[y] + Rest[y])/2 Differences[x]];
> xx = Subdivide[-20., 20., 400000];
> rho = normalPDF[xx, mu, sigma];
> ```
> *Look for:* a grid dense enough that the moment integrals are accurate.
>
> **Step 3 — First moment.** $\mathbb{E}[X]=\int x\,\rho\,dx$. The integrand $x\rho$ is odd, so this should come out to zero.
> ```wolfram
> EX = trapz[xx rho, xx]; (* first moment ~ 0 *)
> Print["E[X] = ", EX, " (exact 0)"];
> ```
> *Look for:* a value near $0$ (about $10^{-17}$, i.e. roundoff) — the numerical echo of the odd-integrand argument.
>
> **Step 4 — Second moment (variance).** $V[X]=\int (x-\mu)^2\rho\,dx$ with $\mu=0$. The integration-by-parts derivation predicts $\tfrac12$.
> ```wolfram
> VX = trapz[(xx - mu)^2 rho, xx]; (* second moment ~ 1/2 *)
> Print["V[X] = ", VX, " (exact sigma^2 = ", N[sigma^2], ")"];
> ```
> *Look for:* $0.5=\sigma^2$, matching the hand calculation.
>
> [!check]- Checks (optional — we are only confirming the hand calculation)
> - $\mathbb{E}[X]\approx 0$ (to roundoff)
> - $V[X]\approx 0.5=\sigma^2$
## Problem 3 — The CDF by a Maclaurin series
**What the problem actually asks.** Approximate the
[CDF](https://en.wikipedia.org/wiki/Cumulative_distribution_function) value
$F(1)=\int_{-\infty}^{1}\rho(x;0,1/\sqrt2)\,dx$ by expanding the integrand as a
[Maclaurin series](https://en.wikipedia.org/wiki/Taylor_series) and integrating term by term, then
say **how many terms** are needed to reach the printed value $0.92135039647485743467$.
**The mathematics.**
*The split (recap).* There is no elementary antiderivative, and you cannot integrate a power series
from $-\infty$ — the pieces $\int_{-\infty}^{1}t^{2k}\,dt$ diverge. So peel off the half you already
know and expand only the finite part:
$F(1)=\underbrace{\int_{-\infty}^{0}\rho\,dx}_{=\,1/2}+\int_0^1\rho\,dx
=\tfrac12+\tfrac1{\sqrt\pi}\int_0^1 e^{-t^2}\,dt,$
where the left half is $\tfrac12$ by symmetry of the normalized density.
*The Maclaurin series (the fact you may not remember).* The exponential series is
$e^{u}=\sum_{k\ge0}u^k/k!$, valid for **every** $u$. Substituting $u=-t^2$,
$e^{-t^2}=\sum_{k=0}^{\infty}\frac{(-t^2)^k}{k!}=\sum_{k=0}^{\infty}\frac{(-1)^k t^{2k}}{k!}
=1-t^2+\tfrac{t^4}{2!}-\tfrac{t^6}{3!}+\cdots$
The factorial in the denominator outgrows any power of $t$, so this converges for all $t$ (the
function is *[entire](https://en.wikipedia.org/wiki/Entire_function)*) — which is exactly what makes
integrating term by term on $[0,1]$ legal. Using
$\int_0^1 t^{2k}\,dt=\tfrac1{2k+1}$,
$\int_0^1 e^{-t^2}\,dt=\sum_{k=0}^{\infty}\frac{(-1)^k}{k!\,(2k+1)}
=1-\tfrac13+\tfrac1{2!\cdot5}-\tfrac1{3!\cdot7}+\cdots,$
and therefore $F(1)=\tfrac12+\tfrac1{\sqrt\pi}\sum_{k\ge0}\frac{(-1)^k}{k!(2k+1)}$. (That sum is
precisely $\tfrac12(1+\operatorname{erf}1)$ — you have just rebuilt the
[error function](https://en.wikipedia.org/wiki/Error_function).)
*How many terms? — the alternating-series fact.* The terms $a_k=\frac1{k!(2k+1)}$ are positive,
strictly **decreasing** to $0$ (the $k!$ alone forces this), and the signs alternate. For any such
[alternating series](https://en.wikipedia.org/wiki/Alternating_series_test), the **Leibniz bound**
says the error from stopping after $k=K$ is no bigger than
the *first dropped term*:
$\Big|\,F(1)-\Big[\tfrac12+\tfrac1{\sqrt\pi}\sum_{k=0}^{K}\frac{(-1)^k}{k!(2k+1)}\Big]\Big|
\;\le\;\frac{a_{K+1}}{\sqrt\pi}=\frac{1}{\sqrt\pi\,(K+1)!\,(2K+3)}.$
To pin $20$ digits we just need that bound below $5\times10^{-21}$. Because of the factorial, $a_k$
collapses fast — each term is roughly $k$ times smaller than the one before — so only about
**20 terms** are needed. That factorial decay is why this "slow" method is in fact quick here.
*One catch (a preview of Problems 4–5).* All of that assumes exact arithmetic. In ordinary floating
point you cannot actually display 20 correct digits: the partial sums stop improving near
$10^{-16}$ because [round-off error](https://en.wikipedia.org/wiki/Round-off_error) sets a floor.
Truncation error keeps shrinking; roundoff does not.
### Code Bases
Unlike Problem 2, the code here *is* the deliverable — counting the terms is the question. It runs
the sum twice: once in extended precision (to answer "how many terms for 20 digits?") and once in
plain double precision (to expose the roundoff floor).
> [!example]- Python
>
> **Step 1 — Work in extended precision.** We are chasing 20 digits, past what ordinary double precision can hold, so we use the language's arbitrary-precision arithmetic. Set the working precision *before* reading in the 20-digit target, or it is silently rounded to the default.
> ```python
> import numpy as np
> from math import factorial, erf
> import mpmath as mp
>
> mp.mp.dps = 30 # set precision FIRST
> target = mp.mpf("0.92135039647485743467") # then parse the 20-digit target
> ```
> *Look for:* the working precision is raised before the target value is read in.
>
> **Step 2 — Sanity-check against the closed form.** The series equals $\tfrac12(1+\operatorname{erf}1)$, so compare to that first.
> ```python
> print("(1+erf 1)/2 =", mp.nstr((1 + mp.erf(1)) / 2, 22))
> print("target =", mp.nstr(target, 22))
> ```
> *Look for:* both lines read $0.9213503964748574346\ldots$
>
> **Step 3 — Sum the Maclaurin series and count terms.** Add the terms $k=0,\dots,K$ of $\tfrac12+\tfrac1{\sqrt\pi}\sum(-1)^k/(k!(2k+1))$ and stop at the first $K$ whose partial sum matches the target to 20 digits.
> ```python
> def F1_hp(K):
> s = mp.mpf(0)
> for k in range(K + 1):
> s += mp.mpf((-1)**k) / (mp.factorial(k) * (2*k + 1))
> return mp.mpf("0.5") + s / mp.sqrt(mp.pi)
>
> for K in range(40):
> if abs(F1_hp(K) - target) < mp.mpf("5e-21"):
> print(f"terms needed for the full 20 digits: {K + 1} (k = 0..{K})")
> break
> ```
> *Look for:* **20 terms** ($k=0..19$) — the factorial makes it converge that fast.
>
> **Step 4 — Repeat in double precision to see the roundoff floor.** The same series in ordinary floats cannot reach 20 digits: rounding error sets a floor near $10^{-16}$.
> ```python
> def F1_dp(K):
> s = sum((-1.0)**k / (factorial(k) * (2*k + 1)) for k in range(K + 1))
> return 0.5 + s / np.sqrt(np.pi)
>
> print("double-precision partial sums (note the ~1e-16 floor):")
> for K in [3, 5, 8, 10, 12, 15, 20]:
> print(f" terms={K+1:2d} F(1)~{F1_dp(K):.16f} |err|={abs(F1_dp(K)-float(target)):.1e}")
> ```
> *Look for:* `|err|` falls to about $10^{-16}$ by ~16 terms and then **stops improving** — more terms cannot beat roundoff. This previews Problems 4–5.
>
> [!example]- MATLAB
>
> **Step 1 — Work in extended precision.** We are chasing 20 digits, past what ordinary double precision can hold, so we use the language's arbitrary-precision arithmetic. Set the working precision *before* reading in the 20-digit target, or it is silently rounded to the default.
> ```matlab
> digits(40);
> target = vpa('0.92135039647485743467'); % needs Symbolic Math Toolbox
> ```
> *Look for:* the working precision is raised before the target value is read in.
>
> **Step 2 — Sanity-check against the closed form.** The series equals $\tfrac12(1+\operatorname{erf}1)$, so compare to that first.
> ```matlab
> fprintf('(1+erf 1)/2 = %s\n', char(vpa((1 + erf(sym(1)))/2, 22)));
> fprintf('target = %s\n', char(vpa(target, 22)));
> ```
> *Look for:* both lines read $0.9213503964748574346\ldots$
>
> **Step 3 — Sum the Maclaurin series and count terms.** Add the terms $k=0,\dots,K$ of $\tfrac12+\tfrac1{\sqrt\pi}\sum(-1)^k/(k!(2k+1))$ and stop at the first $K$ whose partial sum matches the target to 20 digits.
> ```matlab
> s = vpa(0);
> for k = 0:40
> s = s + vpa((-1)^k) / (factorial(sym(k)) * (2*k + 1));
> F1 = vpa(1)/2 + s / sqrt(vpa(pi));
> if abs(F1 - target) < vpa('5e-21')
> fprintf('terms needed for 20 digits: %d (k = 0..%d)\n', k+1, k);
> break
> end
> end
> ```
> *Look for:* **20 terms** ($k=0..19$) — the factorial makes it converge that fast.
>
> **Step 4 — Repeat in double precision to see the roundoff floor.** The same series in ordinary floats cannot reach 20 digits: rounding error sets a floor near $10^{-16}$.
> ```matlab
> fprintf('double-precision partial sums:\n');
> for K = [3 5 8 10 12 15 20]
> kk = 0:K;
> F = 0.5 + sum(((-1).^kk) ./ (factorial(kk) .* (2*kk + 1))) / sqrt(pi);
> fprintf(' terms=%2d F=%.16f err=%.1e\n', K+1, F, abs(F - 0.92135039647485743467));
> end
> ```
> *Look for:* `|err|` falls to about $10^{-16}$ by ~16 terms and then **stops improving** — more terms cannot beat roundoff. This previews Problems 4–5.
>
> [!example]- R
>
> **Step 1 — Work in extended precision.** We are chasing 20 digits, past what ordinary double precision can hold, so we use the language's arbitrary-precision arithmetic. Set the working precision *before* reading in the 20-digit target, or it is silently rounded to the default.
> ```r
> has_mpfr <- requireNamespace("Rmpfr", quietly = TRUE) # install.packages("Rmpfr") for 20 digits
> if (has_mpfr) {
> library(Rmpfr)
> prec <- 120 # bits (~36 digits)
> target <- mpfr("0.92135039647485743467", prec)
> }
> ```
> *Look for:* the working precision is raised before the target value is read in.
>
> **Step 2 — Sanity-check against the closed form.** The series equals $\tfrac12(1+\operatorname{erf}1)$, so compare to that first.
> ```r
> if (has_mpfr)
> cat("(1+erf 1)/2 =", format((1 + erf(mpfr(1, prec)))/2, digits = 22),
> "\ntarget =", format(target, digits = 22), "\n")
> ```
> *Look for:* both lines read $0.9213503964748574346\ldots$
>
> **Step 3 — Sum the Maclaurin series and count terms.** Add the terms $k=0,\dots,K$ of $\tfrac12+\tfrac1{\sqrt\pi}\sum(-1)^k/(k!(2k+1))$ and stop at the first $K$ whose partial sum matches the target to 20 digits.
> ```r
> if (has_mpfr) {
> s <- mpfr(0, prec)
> for (k in 0:40) {
> s <- s + mpfr((-1)^k, prec) / (factorialMpfr(k) * (2*k + 1))
> F1 <- mpfr(0.5, prec) + s / sqrt(Const("pi", prec))
> if (abs(F1 - target) < mpfr("5e-21", prec)) {
> cat(sprintf("terms needed for 20 digits: %d (k = 0..%d)\n", k + 1, k)); break
> }
> }
> } else cat("(20-digit count needs Rmpfr; showing the double-precision floor below)\n")
> ```
> *Look for:* **20 terms** ($k=0..19$) — the factorial makes it converge that fast.
>
> **Step 4 — Repeat in double precision to see the roundoff floor.** The same series in ordinary floats cannot reach 20 digits: rounding error sets a floor near $10^{-16}$.
> ```r
> F1_dp <- function(K) 0.5 + sum((-1)^(0:K) / (factorial(0:K) * (2*(0:K) + 1))) / sqrt(pi)
> cat("double-precision partial sums:\n")
> for (K in c(3, 5, 8, 10, 12, 15, 20))
> cat(sprintf(" terms=%2d F=%.16f err=%.1e\n",
> K + 1, F1_dp(K), abs(F1_dp(K) - 0.92135039647485743467)))
> ```
> *Look for:* `|err|` falls to about $10^{-16}$ by ~16 terms and then **stops improving** — more terms cannot beat roundoff. This previews Problems 4–5.
>
> [!example]- Mathematica
>
> **Step 1 — Work in extended precision.** We are chasing 20 digits, past what ordinary double precision can hold, so we use the language's arbitrary-precision arithmetic. Set the working precision *before* reading in the 20-digit target, or it is silently rounded to the default.
> ```wolfram
> target = 0.92135039647485743467`30; (* a 30-digit number *)
> ```
> *Look for:* the working precision is raised before the target value is read in.
>
> **Step 2 — Sanity-check against the closed form.** The series equals $\tfrac12(1+\operatorname{erf}1)$, so compare to that first.
> ```wolfram
> Print["(1+erf 1)/2 = ", N[(1 + Erf[1])/2, 22]];
> Print["target = ", N[target, 22]];
> ```
> *Look for:* both lines read $0.9213503964748574346\ldots$
>
> **Step 3 — Sum the Maclaurin series and count terms.** Add the terms $k=0,\dots,K$ of $\tfrac12+\tfrac1{\sqrt\pi}\sum(-1)^k/(k!(2k+1))$ and stop at the first $K$ whose partial sum matches the target to 20 digits.
> ```wolfram
> F1[K_] := N[1/2 + (1/Sqrt[Pi]) Sum[(-1)^k/(k! (2 k + 1)), {k, 0, K}], 30];
> Do[If[Abs[F1[K] - target] < 5*^-21,
> Print["terms needed for 20 digits: ", K + 1, " (k = 0..", K, ")"]; Break[]], {K, 0, 40}];
> ```
> *Look for:* **20 terms** ($k=0..19$) — the factorial makes it converge that fast.
>
> **Step 4 — Repeat in double precision to see the roundoff floor.** The same series in ordinary floats cannot reach 20 digits: rounding error sets a floor near $10^{-16}$.
> ```wolfram
> F1dp[K_] := 0.5` + Total[Table[(-1.)^k/(k! (2. k + 1)), {k, 0, K}]]/Sqrt[N[Pi]];
> Print["double-precision partial sums:"];
> Do[Print[" terms=", K + 1, " err=",
> Abs[F1dp[K] - 0.92135039647485743467]], {K, {3, 5, 8, 10, 12, 15, 20}}];
> ```
> *Look for:* `|err|` falls to about $10^{-16}$ by ~16 terms and then **stops improving** — more terms cannot beat roundoff. This previews Problems 4–5.
>
> [!check]- Checks (analytic truth)
> - $\tfrac12(1+\operatorname{erf}1)=0.92135039647485743467$
> - $\approx 20$ terms ($k=0..19$) pin all 20 digits in extended precision
> - in double precision the error floors near $10^{-16}$ by ~16 terms — more terms do not help
## Problem 4 — Finite differences from a Vandermonde system
**What the problem actually asks.** Use the
[Vandermonde](https://en.wikipedia.org/wiki/Vandermonde_matrix) formulation to **recreate the
second-order accurate _forward_ [finite-difference](https://en.wikipedia.org/wiki/Finite_difference)
approximation to $f''(x_0)$**, and write portable code that builds it.
**The mathematics.**
*Where the equation comes from.* Sample $f$ at $x_i=x_0+i\Delta x$ and look for weights $c_i$ so that
$\sum_i c_i f(x_i)$ isolates $f^{(m)}(x_0)$.
[Taylor-expand](https://en.wikipedia.org/wiki/Taylor_series) each sample about $x_0$,
$f(x_i)=f(x_0+i\Delta x)=\sum_{k=0}^{\infty}\frac{f^{(k)}(x_0)}{k!}\,(i\Delta x)^k,$
multiply by $c_i$, and sum. The only $i$-dependence is in $i^{\,k}$, so group by derivative order $k$:
$\sum_i c_i f(x_i)=\sum_{k=0}^{\infty}\frac{(\Delta x)^k}{k!}\Big(\underbrace{\textstyle\sum_i c_i\,i^{\,k}}_{=\,\mu_k}\Big)f^{(k)}(x_0).$
We want the right-hand side to be exactly $f^{(m)}(x_0)\,\dfrac{(\Delta x)^m}{m!}$ — i.e. the $k=m$
term should survive with coefficient $1$ and **every other order should cancel**. That is the set of
*moment conditions* (the right-hand side $\delta_{km}$ is the
[Kronecker delta](https://en.wikipedia.org/wiki/Kronecker_delta))
$\mu_k=\sum_i c_i\,i^{\,k}=\delta_{km}=\begin{cases}1,&k=m\\[2pt]0,&k\neq m.\end{cases}$
With $n+1$ nodes there are $n+1$ unknowns $c_i$, so we can impose this for $k=0,1,\dots,n$. Stacking
those rows is the **Vandermonde system** $V\mathbf c=\mathbf e_m$,
$\begin{bmatrix}1&1&\cdots&1\\ 0&1&\cdots&n\\ 0&1^2&\cdots&n^2\\ \vdots&&&\vdots\\ 0&1^n&\cdots&n^n\end{bmatrix}
\begin{bmatrix}c_0\\c_1\\\vdots\\c_n\end{bmatrix}=\mathbf e_m,$
where **row $k$ enforces the order-$k$ condition**, **column $i$ belongs to node $i$**, and the
right side is the unit vector with its $1$ in row $m$.
*A concrete example — three forward nodes.* Take $m=2$ and $i\in\{0,1,2\}$. The three rows read
$k=0:\ c_0+c_1+c_2=0,\qquad k=1:\ c_1+2c_2=0,\qquad k=2:\ c_1+4c_2=1.$
Subtract the $k=1$ row from the $k=2$ row: $2c_2=1$, so $c_2=\tfrac12$; then $c_1=-2c_2=-1$ and
$c_0=-(c_1+c_2)=\tfrac12$. Scaling by $m!=2$ and dividing by $\Delta x^2$,
$f''(x_0)\approx\frac{2}{\Delta x^2}\Big(\tfrac12 f_0-f_1+\tfrac12 f_2\Big)=\frac{f_0-2f_1+f_2}{\Delta x^2}.$
That is the familiar $(1,-2,1)/\Delta x^2$ — obtained not by memory but by **solving three linear
equations**. (Because these three nodes all sit on one side of $x_0$, this particular formula is only
first-order accurate; see the order note next.)
*Getting second order (what the problem wants).* Run the identical procedure with a fourth node,
$i\in\{0,1,2,3\}$ — a $4\times4$ Vandermonde with right side $\mathbf e_2$ — and you get
$(2,-5,4,-1)/\Delta x^2$. The general rule (the assignment's footnote) is that the
[order of accuracy](https://en.wikipedia.org/wiki/Order_of_accuracy) is $p=n+1-m$ with $n+1$ the
number of samples: three forward nodes give order $1$, four give order
$2$; centering the nodes about $x_0$ earns one extra order for free.
### Code Bases
The code base is the general solver plus the one call that reproduces the forward stencil.
> [!example]- Python
>
> **Step 1 — Build the general solver.** This is the matrix from the derivation, assembled and solved directly: row $k$ of `V` is $i^k$ across the nodes, the right-hand side is the unit vector $\mathbf e_m$, and the built-in linear solver returns the weights. The $m!$ factor is applied here; the caller still divides by $\Delta x^m$.
> ```python
> import numpy as np
> from math import factorial
>
> def fd_weights(offsets, m):
> s = np.asarray(offsets, dtype=float)
> n = len(s)
> V = np.vstack([s**k for k in range(n)]) # V[k, i] = s_i^k (Vandermonde)
> rhs = np.zeros(n); rhs[m] = 1.0 # unit vector e_m
> c = np.linalg.solve(V, rhs)
> return factorial(m) * c # scale by m! (still need /dx^m)
> ```
> *Look for:* the matrix `V` is exactly the Vandermonde from the derivation.
>
> **Step 2 — Recreate the forward second derivative.** One call with the four forward offsets and $m=2$ reproduces the stencil the problem asks for.
> ```python
> print(fd_weights([0, 1, 2, 3], m=2)) # -> [ 2. -5. 4. -1.]
> ```
> *Look for:* $(2,-5,4,-1)$, applied as $f''(x_0)\approx(2f_0-5f_1+4f_2-f_3)/\Delta x^2$.
>
> [!example]- MATLAB
>
> **Step 1 — Build the general solver.** This is the matrix from the derivation, assembled and solved directly: row $k$ of `V` is $i^k$ across the nodes, the right-hand side is the unit vector $\mathbf e_m$, and the built-in linear solver returns the weights. The $m!$ factor is applied here; the caller still divides by $\Delta x^m$.
> ```matlab
> % (in a script, define functions at the END of the file)
> function w = fd_weights(offsets, m)
> s = offsets(:).';
> n = numel(s);
> V = zeros(n, n);
> for k = 0:n-1
> V(k+1, :) = s.^k; % V(k,i) = s_i^k (Vandermonde)
> end
> rhs = zeros(n, 1); rhs(m+1) = 1; % unit vector e_m
> w = factorial(m) * (V \ rhs).'; % scaled weights (still need /dx^m)
> end
> ```
> *Look for:* the matrix `V` is exactly the Vandermonde from the derivation.
>
> **Step 2 — Recreate the forward second derivative.** One call with the four forward offsets and $m=2$ reproduces the stencil the problem asks for.
> ```matlab
> disp(fd_weights([0 1 2 3], 2)) % -> [2 -5 4 -1]
> ```
> *Look for:* $(2,-5,4,-1)$, applied as $f''(x_0)\approx(2f_0-5f_1+4f_2-f_3)/\Delta x^2$.
>
> [!example]- R
>
> **Step 1 — Build the general solver.** This is the matrix from the derivation, assembled and solved directly: row $k$ of `V` is $i^k$ across the nodes, the right-hand side is the unit vector $\mathbf e_m$, and the built-in linear solver returns the weights. The $m!$ factor is applied here; the caller still divides by $\Delta x^m$.
> ```r
> fd_weights <- function(offsets, m) {
> s <- offsets
> n <- length(s)
> V <- t(sapply(0:(n - 1), function(k) s^k)) # V[k+1, i] = s_i^k
> rhs <- numeric(n); rhs[m + 1] <- 1 # unit vector e_m
> as.numeric(factorial(m) * solve(V, rhs))
> }
> ```
> *Look for:* the matrix `V` is exactly the Vandermonde from the derivation.
>
> **Step 2 — Recreate the forward second derivative.** One call with the four forward offsets and $m=2$ reproduces the stencil the problem asks for.
> ```r
> print(fd_weights(c(0, 1, 2, 3), 2)) # -> 2 -5 4 -1
> ```
> *Look for:* $(2,-5,4,-1)$, applied as $f''(x_0)\approx(2f_0-5f_1+4f_2-f_3)/\Delta x^2$.
>
> [!example]- Mathematica
>
> **Step 1 — Build the general solver.** This is the matrix from the derivation, assembled and solved directly: row $k$ of `V` is $i^k$ across the nodes, the right-hand side is the unit vector $\mathbf e_m$, and the built-in linear solver returns the weights. The $m!$ factor is applied here; the caller still divides by $\Delta x^m$.
> ```wolfram
> fdWeights[offsets_, m_] := Module[{s = N[offsets], n, V, e},
> n = Length[s];
> V = Table[If[k == 0, ConstantArray[1., n], s^k], {k, 0, n - 1}]; (* V[[k+1, i]] = s_i^k *)
> e = UnitVector[n, m + 1]; (* e_m *)
> m! LinearSolve[V, e]];
> ```
> *Look for:* the matrix `V` is exactly the Vandermonde from the derivation.
>
> **Step 2 — Recreate the forward second derivative.** One call with the four forward offsets and $m=2$ reproduces the stencil the problem asks for.
> ```wolfram
> Print[fdWeights[{0, 1, 2, 3}, 2]]; (* {2., -5., 4., -1.} *)
> ```
> *Look for:* $(2,-5,4,-1)$, applied as $f''(x_0)\approx(2f_0-5f_1+4f_2-f_3)/\Delta x^2$.
>
> [!check]- Checks (analytic truth)
> - forward $f''$ weights $=(2,-5,4,-1)$, i.e. $f''(x_0)\approx(2f_0-5f_1+4f_2-f_3)/\Delta x^2$
> - the $4\times4$ Vandermonde is nonsingular, so the weights are unique
#### Optional — measuring the order (past what the problem asks)
This part is not required; it just confirms the order the derivation predicts and sets up the
[truncation](https://en.wikipedia.org/wiki/Truncation_error)-versus-noise picture of Problem 5.
> [!example]- Python (optional check)
>
> **Measuring the order.** The problem does not ask for this, but since we have the stencil we can *check* how fast its error shrinks. Evaluate it on a function whose second derivative we know ($f=e^x$, so $f''(0)=1$), shrink $\Delta x$, and read the slope of $\log|\text{error}|$ versus $\log\Delta x$.
> ```python
> f = lambda x: np.exp(x) # f'' = e^x, so f''(0) = 1
> def estimate(offsets, m, dx):
> w = fd_weights(offsets, m)
> s = np.asarray(offsets, dtype=float)
> return np.sum(w * f(s*dx)) / dx**m
>
> for offs, name in [([0,1,2], "forward 3-pt"),
> ([0,1,2,3], "forward 4-pt"),
> ([-1,0,1], "centered 3-pt")]:
> dxs = np.array([0.1, 0.05, 0.025, 0.0125])
> errs = np.array([abs(estimate(offs, 2, dx) - 1.0) for dx in dxs])
> p = np.polyfit(np.log(dxs), np.log(errs), 1)[0]
> print(name, "order ~", round(p, 2))
> ```
> *Look for:* forward 3-pt $\approx 1$, forward 4-pt $\approx 2$, centered 3-pt $\approx 2$ — matching $p=n+1-m$, with centering worth a free order.
>
> [!example]- MATLAB
>
> **Measuring the order.** The problem does not ask for this, but since we have the stencil we can *check* how fast its error shrinks. Evaluate it on a function whose second derivative we know ($f=e^x$, so $f''(0)=1$), shrink $\Delta x$, and read the slope of $\log|\text{error}|$ versus $\log\Delta x$.
> ```matlab
> f = @(x) exp(x); % f'' = e^x, so f''(0) = 1
> configs = {[0 1 2], 'forward 3-pt'; [0 1 2 3], 'forward 4-pt'; [-1 0 1], 'centered 3-pt'};
> dxs = [0.1 0.05 0.025 0.0125];
> for r = 1:size(configs,1)
> offs = configs{r,1};
> errs = arrayfun(@(dx) abs(estimate(offs, 2, dx, f) - 1), dxs);
> p = polyfit(log(dxs), log(errs), 1);
> fprintf('%-14s order ~ %.2f\n', configs{r,2}, p(1));
> end
> % estimate() is a local function at the end of the file:
> % function val = estimate(offsets, m, dx, f)
> % val = sum(fd_weights(offsets, m) .* f(offsets * dx)) / dx^m;
> % end
> ```
> *Look for:* forward 3-pt $\approx 1$, forward 4-pt $\approx 2$, centered 3-pt $\approx 2$ — matching $p=n+1-m$, with centering worth a free order.
>
> [!example]- R
>
> **Measuring the order.** The problem does not ask for this, but since we have the stencil we can *check* how fast its error shrinks. Evaluate it on a function whose second derivative we know ($f=e^x$, so $f''(0)=1$), shrink $\Delta x$, and read the slope of $\log|\text{error}|$ versus $\log\Delta x$.
> ```r
> f <- function(x) exp(x) # f'' = e^x, f''(0) = 1
> estimate <- function(offsets, m, dx) sum(fd_weights(offsets, m) * f(offsets * dx)) / dx^m
> configs <- list(list(c(0,1,2), "forward 3-pt"),
> list(c(0,1,2,3), "forward 4-pt"),
> list(c(-1,0,1), "centered 3-pt"))
> dxs <- c(0.1, 0.05, 0.025, 0.0125)
> for (cfg in configs) {
> errs <- sapply(dxs, function(dx) abs(estimate(cfg[[1]], 2, dx) - 1))
> cat(sprintf("%-14s order ~ %.2f\n", cfg[[2]], coef(lm(log(errs) ~ log(dxs)))[2]))
> }
> ```
> *Look for:* forward 3-pt $\approx 1$, forward 4-pt $\approx 2$, centered 3-pt $\approx 2$ — matching $p=n+1-m$, with centering worth a free order.
>
> [!example]- Mathematica
>
> **Measuring the order.** The problem does not ask for this, but since we have the stencil we can *check* how fast its error shrinks. Evaluate it on a function whose second derivative we know ($f=e^x$, so $f''(0)=1$), shrink $\Delta x$, and read the slope of $\log|\text{error}|$ versus $\log\Delta x$.
> ```wolfram
> f[x_] := Exp[x]; (* f'' = e^x, f''(0) = 1 *)
> estimate[offsets_, m_, dx_] := Total[fdWeights[offsets, m] f[offsets dx]]/dx^m;
> order[offs_] := Log2[Abs[estimate[offs, 2, 0.1] - 1]/Abs[estimate[offs, 2, 0.05] - 1]];
> Print["forward 3-pt order ~ ", order[{0, 1, 2}]];
> Print["forward 4-pt order ~ ", order[{0, 1, 2, 3}]];
> Print["centered 3-pt order ~ ", order[{-1, 0, 1}]];
> ```
> *Look for:* forward 3-pt $\approx 1$, forward 4-pt $\approx 2$, centered 3-pt $\approx 2$ — matching $p=n+1-m$, with centering worth a free order.
>
![[hw1_p4_order_plot.png]]
## Problem 5 — Missing data and noisy sampling
**What the problem actually asks.** Sample $f(x)=e^{-x^2}$ on a uniform grid, **lose one interior
reading**, rebuild the stencil on the points that remain to estimate $f''(0)$ and compare to the
exact $-2$; then corrupt the samples with noise, recompute, and explain why the estimate is so
sensitive.
**The mathematics.**
*Why a missing sample is a problem.* Problem 4's tidy formulas assume the samples sit in a fixed
pattern — equally spaced, often symmetric about $x_0$. Here we sample on the uniform grid
$x_i=i\Delta x$ for $i=-1,0,1,2,3$ but **lose the reading at $i=1$**. The points that remain,
$\{-\Delta x,\,0,\,2\Delta x,\,3\Delta x\}$, have a hole in them — they are no longer any textbook
stencil — so there is no memorized formula to reach for.
*The fix is to not reach for a memorized formula.* Look back at Problem 4: the derivation never
required even spacing or symmetry. It only used each node's offset $s_i$. So we hand it the offsets
we *actually* have, $s\in\{-1,0,2,3\}$, and solve for a custom stencil. (Two of those offsets are on
the left and two on the right of $x_0$, so this is a two-sided stencil — which the general system
handles with no change.)
*Building the stencil concretely — where the missing point goes.* In the Vandermonde picture each
**node is a column** and each **derivative-order condition is a row**. If nothing were lost, the full
grid $i\in\{-1,0,1,2,3\}$ would give five nodes and a $5\times5$ system enforcing orders $k=0,\dots,4$:
$\underbrace{\begin{bmatrix}1&1&1&1&1\\ -1&0&1&2&3\\ 1&0&1&4&9\\ -1&0&1&8&27\\ 1&0&1&16&81\end{bmatrix}}_{\textstyle s\,=\,-1,\ 0,\ 1,\ 2,\ 3}
\begin{bmatrix}c_{-1}\\c_{0}\\c_{1}\\c_{2}\\c_{3}\end{bmatrix}=\begin{bmatrix}0\\0\\1\\0\\0\end{bmatrix}.$
Losing the reading at $i=1$ deletes that node — and with it the **middle column** (the one under
$s=1$, namely $(1,1,1,1,1)^{\mathsf T}$) together with its unknown $c_1$. That leaves only four
unknowns, so we can no longer demand five conditions: we also drop the **bottom row**, the
highest-order one ($k=4$) that we can no longer afford. What survives is the $4\times4$ system on the
offsets $\{-1,0,2,3\}$:
$\begin{bmatrix}1&1&1&1\\ -1&0&2&3\\ 1&0&4&9\\ -1&0&8&27\end{bmatrix}
\begin{bmatrix}c_{-1}\\c_{0}\\c_{2}\\c_{3}\end{bmatrix}=\begin{bmatrix}0\\0\\1\\0\end{bmatrix}.$
So a **missing sample is a missing column**, and paying for it with a dropped row costs exactly one
order of accuracy (the $p=n+1-m$ rule, now with a smaller $n+1$). Solving (the same moves as Problem
4's example, one size larger) gives $\mathbf c=\tfrac1{12}(5,-8,4,-1)$, so the scaled weights are
$w=m!\,\mathbf c=\tfrac16(5,-8,4,-1)$ and
$f''(0)\approx\frac{1}{\Delta x^2}\sum_i w_i\,y_i.$
On the clean samples (with $\Delta x=\tfrac1{12}$) this returns $-2.006$ against the exact
$f''(0)=-2$ — about $0.3\%$ off. **So far there is no noise: the missing-data stencil simply works.**
*Now add noise — and where the numbers come from.* The corruption is itself a draw from the
[normal distribution](https://en.wikipedia.org/wiki/Normal_distribution) of Problem 1, recentered at
$0$: $\epsilon_i\sim\mathcal N(0,\sigma_i^2)$. The only thing
to decide is the spread $\sigma_i$, and the problem fixes it indirectly by saying "$99\%$ of the
noise should fall within $3\%$ of $y_i$." To turn that sentence into a number, recall a basic fact
about *any* normal distribution: a fixed fraction of its draws lands within a fixed number of
[standard deviations](https://en.wikipedia.org/wiki/Standard_deviation) of the center — about $68\%$
within $\pm1\sigma$, $95\%$ within $\pm2\sigma$, $99.7\%$ within $\pm3\sigma$ (the
[empirical rule](https://en.wikipedia.org/wiki/68%E2%80%9395%E2%80%9399.7_rule)). We want the
$99\%$ cutoff, and we want the band measured as a
percent of the true value.
"Within $3\%$ of $y_iquot; is the symmetric band $|\epsilon_i|<0.03\,y_i$, of half-width $0.03\,y_i$
around $0$. To trap $99\%$ of a normal distribution's mass in a symmetric band, the band must extend
$z_{0.995}\approx2.576$ standard deviations on each side. The subscript is $0.995$, not $0.99$,
because leaving $1\%$ *outside* a two-sided band puts $0.5\%$ in each tail, so the edge sits at the
$99.5$th [percentile](https://en.wikipedia.org/wiki/Percentile) — in the language of Problem 3, the
point where the standard-normal CDF reaches
$0.995$. Matching the band half-width to that reach,
$0.03\,y_i=z_{0.995}\,\sigma_i\quad\Longrightarrow\quad
\sigma_i=\frac{0.03\,y_i}{2.576}\approx0.0116\,y_i.$
So each sample is given noise whose standard deviation is about $1.2\%$ of its own value — larger
samples get proportionally larger noise. The $\pm3\%$ is the visible "error bar"; $\sigma_i$ is the
dial that produces it.
Now feed the noisy samples through the **same** weights and repeat over many random draws (a
[Monte Carlo](https://en.wikipedia.org/wiki/Monte_Carlo_method) experiment). The $1/\Delta x^2$ out
front magnifies that $\sim1.2\%$ input into roughly **$140\%$** scatter in $f''(0)$: the data barely
moves, yet the second derivative is meaningless. The two-panel figure shows both sides — a tiny
wiggle in, chaos out.
### Code Bases
The build follows the math exactly: first *initialize* the grid, the samples, and the linear system
$V\mathbf c=\mathbf e_m$; then call the **built-in solver** for the weights; then add noise and push
it through.
> [!example]- Python
>
> **Step 1 — Initialize the gapped grid (initialization).** The surviving offsets are $s=(-1,0,2,3)$; the nodes are $x_i=s_i\,\Delta x$ and the samples are $y_i=f(x_i)$ with $f(x)=e^{-x^2}$.
> ```python
> import numpy as np
> from math import factorial
>
> f = lambda x: np.exp(-x**2)
> dx = 1/12
> s = np.array([-1, 0, 2, 3]) # surviving offsets s_i
> x = s * dx # nodes x_i = s_i * dx
> y = f(x) # samples y_i = f(x_i)
> ```
> *Look for:* four nodes with a hole where $i=1$ used to be.
>
> **Step 2 — Initialize the linear system (initialization).** Build the Vandermonde matrix $V_{ki}=s_i^{\,k}$ and the unit right-hand side $\mathbf e_m$ (a $1$ in row $m=2$). These are exactly the matrix and vector written out in the math above.
> ```python
> m = 2
> n = len(s)
> V = np.vstack([s.astype(float)**k for k in range(n)]) # V[k, i] = s_i^k
> e = np.zeros(n); e[m] = 1.0 # e_m
> ```
> *Look for:* `V` has rows $s^0,s^1,s^2,s^3$ — the $4\times4$ system from the derivation.
>
> **Step 3 — Solve with the built-in linear solver.** The built-in linear solver solves $V\mathbf c=\mathbf e_m$ for the moment coefficients $c_i$; scaling by $m!$ gives the stencil weights, and the estimate is $f''(0)\approx\frac{1}{\Delta x^2}\sum_i w_i y_i$.
> ```python
> c = np.linalg.solve(V, e) # V c = e_m
> w = factorial(m) * c # weights w_i
> fpp_clean = np.sum(w * y) / dx**2
> print(w, fpp_clean) # w = (5,-8,4,-1)/6 ; f''(0) ~ -2.006
> ```
> *Look for:* $w=\tfrac16(5,-8,4,-1)$ and a clean estimate $-2.006$ against the exact $-2$.
>
> **Step 4 — Set the noise level.** Choose $\sigma_i$ so $99\%$ of the noise lands within $3\%$ of $y_i$: $0.03\,y_i=z_{0.995}\,\sigma_i$ with $z_{0.995}\approx2.576$, i.e. $\sigma_i\approx0.0116\,y_i$.
> ```python
> z995 = 2.5758293035489004 # 0.995 quantile of N(0,1)
> sigma = 0.03 * np.abs(y) / z995 # sigma_i ~ 1.16% of y_i
> ```
> *Look for:* the noise is only about $1.2\%$ of each sample.
>
> **Step 5 — Push noisy samples through the same weights.** Draw many noisy sample sets $z_i=y_i+\epsilon_i$ and recompute $f''(0)=\frac1{\Delta x^2}\sum_i w_i z_i$ for each; the spread of those estimates is the noise sensitivity. (The remaining lines of the script draw the two-panel figure below.)
> ```python
> rng = np.random.default_rng(0)
> Z = y + rng.normal(0, sigma, size=(20000, len(y))) # noisy sample sets z_i
> est = Z @ w / dx**2 # noisy f''(0) estimates
> print(est.std()) # ~ 2.86 (143% of |exact 2|)
> ```
> *Look for:* a spread of about $143\%$ of $|{-2}|$, even though the inputs moved only $\sim1\%$.
>
> [!example]- MATLAB
>
> **Step 1 — Initialize the gapped grid (initialization).** The surviving offsets are $s=(-1,0,2,3)$; the nodes are $x_i=s_i\,\Delta x$ and the samples are $y_i=f(x_i)$ with $f(x)=e^{-x^2}$.
> ```matlab
> f = @(x) exp(-x.^2);
> fpp = @(x) (4*x.^2 - 2) .* exp(-x.^2);
> dx = 1/12;
> s = [-1 0 2 3]; % surviving offsets s_i
> x = s * dx; % nodes x_i = s_i * dx
> y = f(x); % samples y_i = f(x_i)
> ```
> *Look for:* four nodes with a hole where $i=1$ used to be.
>
> **Step 2 — Initialize the linear system (initialization).** Build the Vandermonde matrix $V_{ki}=s_i^{\,k}$ and the unit right-hand side $\mathbf e_m$ (a $1$ in row $m=2$). These are exactly the matrix and vector written out in the math above.
> ```matlab
> m = 2; n = numel(s);
> V = zeros(n, n);
> for k = 0:n-1
> V(k+1, :) = s.^k; % V(k,i) = s_i^k
> end
> e = zeros(n, 1); e(m+1) = 1; % e_m (1-based)
> ```
> *Look for:* `V` has rows $s^0,s^1,s^2,s^3$ — the $4\times4$ system from the derivation.
>
> **Step 3 — Solve with the built-in linear solver.** The built-in linear solver solves $V\mathbf c=\mathbf e_m$ for the moment coefficients $c_i$; scaling by $m!$ gives the stencil weights, and the estimate is $f''(0)\approx\frac{1}{\Delta x^2}\sum_i w_i y_i$.
> ```matlab
> c = V \ e; % V c = e_m
> w = factorial(m) * c.'; % f''(0) ~ sum w_i y_i / dx^2
> fpp_clean = sum(w .* y) / dx^2;
> fprintf('weights w = [%g %g %g %g]\n', w);
> fprintf('clean f''''(0) = %.5f exact %.5f\n', fpp_clean, fpp(0));
> ```
> *Look for:* $w=\tfrac16(5,-8,4,-1)$ and a clean estimate $-2.006$ against the exact $-2$.
>
> **Step 4 — Set the noise level.** Choose $\sigma_i$ so $99\%$ of the noise lands within $3\%$ of $y_i$: $0.03\,y_i=z_{0.995}\,\sigma_i$ with $z_{0.995}\approx2.576$, i.e. $\sigma_i\approx0.0116\,y_i$.
> ```matlab
> z995 = 2.5758293035489004; % 0.995 quantile of N(0,1)
> sigma = 0.03 * abs(y) / z995; % sigma_i ~ 1.2% of y_i
> ```
> *Look for:* the noise is only about $1.2\%$ of each sample.
>
> **Step 5 — Push noisy samples through the same weights.** Draw many noisy sample sets $z_i=y_i+\epsilon_i$ and recompute $f''(0)=\frac1{\Delta x^2}\sum_i w_i z_i$ for each; the spread of those estimates is the noise sensitivity. (The remaining lines of the script draw the two-panel figure below.)
> ```matlab
> rng(0);
> trials = 20000;
> Z = y + randn(trials, n) .* sigma; % implicit expansion over rows
> est = (Z * w.') / dx^2;
> fprintf('noisy f''''(0): mean %.3f std %.3f (%.0f%% of |exact|)\n', ...
> mean(est), std(est), 100*std(est)/abs(fpp(0)));
> ```
> *Look for:* a spread of about $143\%$ of $|{-2}|$, even though the inputs moved only $\sim1\%$.
>
> [!example]- R
>
> **Step 1 — Initialize the gapped grid (initialization).** The surviving offsets are $s=(-1,0,2,3)$; the nodes are $x_i=s_i\,\Delta x$ and the samples are $y_i=f(x_i)$ with $f(x)=e^{-x^2}$.
> ```r
> f <- function(x) exp(-x^2)
> fpp <- function(x) (4 * x^2 - 2) * exp(-x^2)
> dx <- 1/12
> s <- c(-1, 0, 2, 3) # surviving offsets s_i
> x <- s * dx # nodes x_i = s_i * dx
> y <- f(x) # samples y_i = f(x_i)
> ```
> *Look for:* four nodes with a hole where $i=1$ used to be.
>
> **Step 2 — Initialize the linear system (initialization).** Build the Vandermonde matrix $V_{ki}=s_i^{\,k}$ and the unit right-hand side $\mathbf e_m$ (a $1$ in row $m=2$). These are exactly the matrix and vector written out in the math above.
> ```r
> m <- 2; n <- length(s)
> V <- t(sapply(0:(n - 1), function(k) s^k)) # V[k+1, i] = s_i^k
> e <- numeric(n); e[m + 1] <- 1 # e_m
> ```
> *Look for:* `V` has rows $s^0,s^1,s^2,s^3$ — the $4\times4$ system from the derivation.
>
> **Step 3 — Solve with the built-in linear solver.** The built-in linear solver solves $V\mathbf c=\mathbf e_m$ for the moment coefficients $c_i$; scaling by $m!$ gives the stencil weights, and the estimate is $f''(0)\approx\frac{1}{\Delta x^2}\sum_i w_i y_i$.
> ```r
> cc <- solve(V, e) # V c = e_m
> w <- as.numeric(factorial(m) * cc) # f''(0) ~ sum w_i y_i / dx^2
> fpp_clean <- sum(w * y) / dx^2
> cat("weights w =", w, "\n")
> cat(sprintf("clean f''(0) = %.5f exact %.5f\n", fpp_clean, fpp(0)))
> ```
> *Look for:* $w=\tfrac16(5,-8,4,-1)$ and a clean estimate $-2.006$ against the exact $-2$.
>
> **Step 4 — Set the noise level.** Choose $\sigma_i$ so $99\%$ of the noise lands within $3\%$ of $y_i$: $0.03\,y_i=z_{0.995}\,\sigma_i$ with $z_{0.995}\approx2.576$, i.e. $\sigma_i\approx0.0116\,y_i$.
> ```r
> z995 <- 2.5758293035489004 # 0.995 quantile of N(0,1)
> sigma <- 0.03 * abs(y) / z995 # sigma_i ~ 1.2% of y_i
> ```
> *Look for:* the noise is only about $1.2\%$ of each sample.
>
> **Step 5 — Push noisy samples through the same weights.** Draw many noisy sample sets $z_i=y_i+\epsilon_i$ and recompute $f''(0)=\frac1{\Delta x^2}\sum_i w_i z_i$ for each; the spread of those estimates is the noise sensitivity. (The remaining lines of the script draw the two-panel figure below.)
> ```r
> set.seed(0)
> trials <- 20000
> Z <- matrix(rnorm(trials * n), trials, n)
> Z <- sweep(Z, 2, sigma, "*") # scale columns by sigma_i
> Z <- sweep(Z, 2, y, "+") # add the true samples y_i
> est <- as.numeric(Z %*% w) / dx^2
> cat(sprintf("noisy f''(0): mean %.3f std %.3f (%.0f%% of |exact|)\n",
> mean(est), sd(est), 100 * sd(est) / abs(fpp(0))))
> ```
> *Look for:* a spread of about $143\%$ of $|{-2}|$, even though the inputs moved only $\sim1\%$.
>
> [!example]- Mathematica
>
> **Step 1 — Initialize the gapped grid (initialization).** The surviving offsets are $s=(-1,0,2,3)$; the nodes are $x_i=s_i\,\Delta x$ and the samples are $y_i=f(x_i)$ with $f(x)=e^{-x^2}$.
> ```wolfram
> f[x_] := Exp[-x^2];
> fpp[x_] := (4 x^2 - 2) Exp[-x^2];
> dx = N[1/12];
> s = {-1, 0, 2, 3}; (* surviving offsets *)
> x = s dx; (* nodes x_i = s_i dx *)
> y = f[x]; (* samples y_i = f(x_i) *)
> ```
> *Look for:* four nodes with a hole where $i=1$ used to be.
>
> **Step 2 — Initialize the linear system (initialization).** Build the Vandermonde matrix $V_{ki}=s_i^{\,k}$ and the unit right-hand side $\mathbf e_m$ (a $1$ in row $m=2$). These are exactly the matrix and vector written out in the math above.
> ```wolfram
> m = 2; n = Length[s];
> V = Table[If[k == 0, ConstantArray[1., n], N[s]^k], {k, 0, n - 1}]; (* V[[k+1, i]] = s_i^k *)
> e = UnitVector[n, m + 1]; (* e_m *)
> ```
> *Look for:* `V` has rows $s^0,s^1,s^2,s^3$ — the $4\times4$ system from the derivation.
>
> **Step 3 — Solve with the built-in linear solver.** The built-in linear solver solves $V\mathbf c=\mathbf e_m$ for the moment coefficients $c_i$; scaling by $m!$ gives the stencil weights, and the estimate is $f''(0)\approx\frac{1}{\Delta x^2}\sum_i w_i y_i$.
> ```wolfram
> c = LinearSolve[V, e]; (* V c = e_m *)
> w = m! c; (* f''(0) ~ Sum w_i y_i / dx^2 *)
> fppClean = Total[w y]/dx^2;
> Print["weights w = ", w];
> Print["clean f''(0) = ", fppClean, " exact ", fpp[0]];
> ```
> *Look for:* $w=\tfrac16(5,-8,4,-1)$ and a clean estimate $-2.006$ against the exact $-2$.
>
> **Step 4 — Set the noise level.** Choose $\sigma_i$ so $99\%$ of the noise lands within $3\%$ of $y_i$: $0.03\,y_i=z_{0.995}\,\sigma_i$ with $z_{0.995}\approx2.576$, i.e. $\sigma_i\approx0.0116\,y_i$.
> ```wolfram
> z995 = 2.5758293035489004; (* 0.995 quantile of N(0,1) *)
> sigma = 0.03 Abs[y]/z995; (* sigma_i ~ 1.2% of y_i *)
> ```
> *Look for:* the noise is only about $1.2\%$ of each sample.
>
> **Step 5 — Push noisy samples through the same weights.** Draw many noisy sample sets $z_i=y_i+\epsilon_i$ and recompute $f''(0)=\frac1{\Delta x^2}\sum_i w_i z_i$ for each; the spread of those estimates is the noise sensitivity. (The remaining lines of the script draw the two-panel figure below.)
> ```wolfram
> SeedRandom[0];
> trials = 20000;
> Z = Table[y + RandomVariate[NormalDistribution[0, 1], n] sigma, {trials}];
> est = (Z . w)/dx^2;
> Print["noisy f''(0): mean ", Mean[est], " std ", StandardDeviation[est],
> " (", Round[100 StandardDeviation[est]/Abs[fpp[0]]], "% of |exact|)"];
> ```
> *Look for:* a spread of about $143\%$ of $|{-2}|$, even though the inputs moved only $\sim1\%$.
>
> [!check]- Checks (analytic truth)
> - stencil weights $w=\tfrac16(5,-8,4,-1)$; clean $f''(0)\approx-2.006$ (exact $-2$, ~0.3% off)
> - noisy estimate std $\approx 143\%$ of $|{-2}|$ from only ~1.2% input noise — the $1/\Delta x^2$ effect
![[hw1_p5_noise.png]]