The Dirichlet problem asks for a function that satisfies a differential equation inside a region and takes prescribed values on its boundary. For Laplace's equation the interior condition is $\Delta u = 0$, and in one dimension that reads $u'' = 0$, whose solution is the straight line through the two boundary values. There is nothing left to compute. The one-dimensional model problem worth solving numerically is therefore Poisson's equation, $$-u''(x) = f(x), \qquad 0 < x < 1, \qquad u(0) = \alpha, \quad u(1) = \beta.$$ It has the same second-difference structure, the same conditioning, and the same convergence theory as the two-dimensional problem, at a hundredth of the cost, so it is the right place to develop the tools. From Chapter 5 on we set $f = 0$ and let the second dimension supply the difficulty. ### Knowing the answer in advance A convergence claim needs an error, and an error needs a truth to measure against. Reference solutions computed by the same family of methods at higher resolution are the usual substitute, and they are unreliable in exactly the situations one wants to study: near a singularity, at high frequency, or wherever the method being tested is systematically wrong. Here we use manufactured solutions throughout. Choose $u$ first, differentiate it twice by hand, and *define* the data by $$f := -u'', \qquad \alpha := u(0), \qquad \beta := u(1).$$ Then $u$ is the exact solution of the resulting boundary value problem, by construction. The error of any computed approximation is known at every point, to full precision, at every resolution, with no reference computation and no extrapolation. The obvious limitation is that we are testing a method against a solution we chose. A method can be flattered by a convenient choice, so the five manufactured solutions in these labs differ deliberately in how hard they are to resolve: an entire function of moderate amplitude, one with eight wavelengths across the interval, a narrow Gaussian, a boundary layer of width $1/50$, and one that is $C^2$ but not $C^3$. Anything that is claimed below should be checked against all five, and the last two are where the claims get interesting. ### The three-point scheme Put $N$ equal intervals of length $h = 1/N$ on $[0,1]$, write $u_i$ for the value at $x_i = ih$, and replace the second derivative by the second difference, $$u''(x_i) \;\approx\; \frac{u_{i-1} - 2u_i + u_{i+1}}{h^2}.$$ A Taylor expansion of the right-hand side gives $u''(x_i) + \frac{h^2}{12}u''''(x_i) + O(h^4)$, so the substitution is exact for quadratics and wrong by a term proportional to $h^2$ for anything else. Imposing the resulting equation at the $N-1$ interior nodes, and moving the two known boundary values to the right-hand side, gives $$\frac{1}{h^2}\begin{pmatrix} 2 & -1 & & \\ -1 & 2 & \ddots & \\ & \ddots & \ddots & -1 \\ & & -1 & 2\end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \\ \vdots \\ u_{N-1}\end{pmatrix} = \begin{pmatrix} f_1 + \alpha/h^2 \\ f_2 \\ \vdots \\ f_{N-1} + \beta/h^2 \end{pmatrix}.$$ The matrix is symmetric positive definite and tridiagonal, so Gaussian elimination without pivoting needs one forward and one backward sweep, $O(N)$ operations and $O(N)$ storage. Nothing about one dimension is expensive. Two things are worth noticing before moving on. The error is not noise: it is a smooth function of $x$, of one sign over long stretches, with a shape that does not change as $N$ grows. Only its amplitude changes. And the amplitude is a strong function of which solution you ask for. At $N = 32$ the smooth case is already accurate to a few parts in $10^4$, while the boundary layer is not resolved at all, because a mesh that never places a point inside the layer cannot know it is there. ### Why the error behaves that way Let $e_i = u_i - u(x_i)$ be the nodal error. Applying the discrete operator to $e$ and using the Taylor expansion above, $$\frac{e_{i-1} - 2e_i + e_{i+1}}{h^2} \;=\; \frac{h^2}{12}\,u''''(x_i) + O(h^4),$$ with $e_0 = e_N = 0$. The error therefore solves the *same* discrete problem, with the truncation error as its source. Since the discrete solution operator is bounded independently of $h$, the error inherits the size of that source: $$e \;\approx\; \frac{h^2}{12}\, w, \qquad -w'' = u'''', \quad w(0) = w(1) = 0.$$ That is the whole of the classical theory in one line, and it predicts two things the lab above shows. The error is a fixed smooth function, $w/12$, scaled by $h^2$, which is why its shape is stable under refinement. And its size is set by $u''''$, which is why the narrow Gaussian and the boundary layer are hard: their fourth derivatives are enormous. The observed order settles at exactly 2 for the first four solutions. For `lowreg`, where $u''''$ does not exist at $x = 1/2$, the argument above does not apply, and the observed order is lower. This is the first appearance of a theme that runs through the rest of the tutorial: convergence rates are statements about the smoothness of the solution, not about the method alone. Second order is a modest rate. Each extra digit of accuracy costs a factor of $\sqrt{10} \approx 3.2$ in $N$, so three more digits costs a factor of 32 in work. The next chapter asks how far that can be pushed before something other than the discretization gets in the way.