43.10.06 · numerical-analysis / 10-numerical-odes

Finite-difference methods for two-point boundary-value problems

shipped3 tiersLean: none

Anchor (Master): LeVeque 2007 *Finite Difference Methods for Ordinary and Partial Differential Equations* (SIAM) Ch. 2 (the steady-state BVP in full: variable coefficients $(\kappa u')'$, the discrete Green's function and comparison-function stability bound, the half-order Neumann pitfall, and Newton on the discrete nonlinear system); Keller 1992 *Numerical Methods for Two-Point Boundary-Value Problems* (Dover) Ch. 1-3 (shooting, multiple shooting, and the finite-difference method with its convergence theory); Stoer-Bulirsch 2002 *Introduction to Numerical Analysis* 3e (Springer) §7.3 (shooting, multiple shooting, and difference methods for boundary-value problems)

Intuition Beginner

An initial-value problem hands you everything at one end: a starting height and a starting slope, and you march forward. A two-point boundary-value problem is stingier. It tells you the height at the left end of an interval and the height at the right end, and asks for the curve in between that obeys a differential equation and threads both pins. You know where the rope is nailed at both walls; you must find the sag.

This change of question changes the method. You can no longer just walk forward, because you are not given the starting slope. Two honest tricks get around this. The first is to guess the missing slope, march forward as if it were an initial-value problem, and see where you land at the far wall. If you overshoot, lower the guess; if you undershoot, raise it. Correcting the guess until you hit the right end on the nose is the shooting method, named for aiming a cannon.

The second trick refuses to march at all. Lay down a row of grid points across the interval and, at each one, replace the second derivative by a simple comparison of a point with its two neighbours: how far it sits below their average. Writing that relation at every interior point gives a tidy system of linear equations, one per unknown, each touching only a point and its two neighbours. That banded system is the finite-difference method, and you solve it all at once instead of stepping through it.

Two questions then matter, the same two that matter for every approximation. Does the grid answer approach the true answer as the grid gets finer? And can the method invent a value the true solution would never take? For the steady problems here the answers are yes-it-converges, at a rate set by the spacing squared, and no-it-cannot, because of a discrete echo of a fact about the smooth problem: with no source pushing, the solution sits between its boundary values and never overshoots them.

Visual Beginner

The picture is a single interior grid point with its two neighbours, and the weights the second-difference stencil attaches to each. The centre point is compared against the average of the point to its left and the point to its right.

The stencil places a weight of two on the centre and a weight of minus one on each neighbour, all divided by the square of the spacing. When the source is zero, this forces the centre value to be exactly the average of its left and right neighbours — a straight line through the two.

          -1        +2        -1
   (j-1) ------- ( j ) ------- (j+1)
     |             |             |
   U_{j-1}        U_j         U_{j+1}

position	weight in the stencil	meaning
centre point	$+ 2$	the value being solved for
each of the two neighbours	$- 1$	pulls the centre toward the neighbour average
anything else	$0$	the stencil is local, so the system is tridiagonal

The takeaway: the entire finite-difference method is this one small picture, copied to every interior grid point in a single line of points. The two endpoint dots are not unknowns — their values are handed to you by the boundary data — so the unknowns are only the interior dots, and each one is tied to just two neighbours. That is why the linear system is tridiagonal and cheap to solve.

Worked example Beginner

Let us solve the smallest interesting case by hand. Take the interval from $0$ to $1$ and put down a grid with spacing $h = 1/4$ , so there are three interior points at $x_{1} = 1/4$ , $x_{2} = 1/2$ , $x_{3} = 3/4$ . We solve $- u^{''} = f$ with a constant source $f = 8$ , and boundary values $u (0) = 0$ and $u (1) = 0$ .

The true solution of $- u^{''} = 8$ with both ends pinned at $0$ is the downward parabola $u (x) = 4 x (1 - x)$ , so we have an exact answer to check against: $u (1/4) = 3/4$ , $u (1/2) = 1$ , $u (3/4) = 3/4$ .

Step 1. Write the stencil at each interior point. The relation $- u^{''} = f$ becomes $(- U_{j - 1} + 2 U_{j} - U_{j + 1}) / h^{2} = 8$ . With $h = 1/4$ , $h^{2} = 1/16$ , so $- U_{j - 1} + 2 U_{j} - U_{j + 1} = 8/16 = 1/2$ at every interior point.

Step 2. Plug in the three points, using $U_{0} = 0$ and $U_{4} = 0$ at the ends. At $x_{1}$ : $2 U_{1} - U_{2} = 1/2$ . At $x_{2}$ : $- U_{1} + 2 U_{2} - U_{3} = 1/2$ . At $x_{3}$ : $- U_{2} + 2 U_{3} = 1/2$ .

Step 3. Use the left-right symmetry of the data: $U_{1} = U_{3}$ . The middle equation becomes $- 2 U_{1} + 2 U_{2} = 1/2$ , so $U_{2} = U_{1} + 1/4$ . The first equation is $2 U_{1} - U_{2} = 1/2$ ; substitute $U_{2}$ : $2 U_{1} - U_{1} - 1/4 = 1/2$ , giving $U_{1} = 3/4$ .

Step 4. Read off the rest: $U_{1} = U_{3} = 3/4$ and $U_{2} = 3/4 + 1/4 = 1$ .

What this tells us: the three grid values are $3/4$ , $1$ , $3/4$ — and they match the exact answer at these points exactly. That is special to this example: the source is constant, so the true solution is a quadratic, and the centred second difference reproduces the second derivative of a quadratic with no error at all. For a general source the grid answer would be close but not exact, with a gap that shrinks like the spacing squared.

Check your understanding Beginner

Formal definition Intermediate+

Consider the linear two-point boundary-value problem on a bounded interval $[a, b]$ , $$ -u''(x) + q(x),u(x) = f(x), \quad x \in (a,b), \qquad u(a) = \alpha, \ \ u(b) = \beta, $$ with $q \in C [a, b]$ , $q (x) \geq 0$ , $f \in C [a, b]$ , and Dirichlet data $α, β$ . The continuous existence and uniqueness of a classical solution $u \in C^{2} [a, b]$ under $q \geq 0$ is taken as input from the boundary-value theory of the underlying ODE; this unit is the numerical treatment. The more general variable-coefficient (divergence/self-adjoint) form is $$ -\big(\kappa(x),u'(x)\big)' + q(x),u(x) = f(x), \qquad \kappa(x) \ge \kappa_0 > 0, $$ which reduces to the first form when $κ \equiv 1$ .

The grid and the second-difference operator. Fix $N \geq 1$ , set the mesh width $h = (b - a) / (N + 1)$ , and place grid points $x_{j} = a + j h$ for $0 \leq j \leq N + 1$ . The interior indices are $1 \leq j \leq N$ , numbering $N$ unknowns $U_{j} \approx u (x_{j})$ ; the boundary values $U_{0} = α$ , $U_{N + 1} = β$ are fixed by the data. For a grid function $V$ the second-difference operator is $$ (D^2 V)j = \frac{V{j-1} - 2V_j + V_{j+1}}{h^2}, $$ and the discretised problem replaces $- u^{''}$ by $- D^{2}$ : $$ -\frac{U_{j-1} - 2U_j + U_{j+1}}{h^2} + q_j U_j = f_j, \qquad q_j = q(x_j), \ f_j = f(x_j), $$ at every interior $j$ , with $U_{0}, U_{N + 1}$ substituted from the data.

The tridiagonal system. Collecting the interior unknowns into $U \in R^{N}$ gives the linear system $$ A_h,U = F, \qquad A_h = \frac{1}{h^2},\mathrm{tridiag}(-1,,2,,-1) + \mathrm{diag}(q_1,\dots,q_N), $$ where $F_{j} = f_{j}$ except in the two boundary-adjacent rows, which carry the moved data: $F_{1} = f_{1} + α / h^{2}$ and $F_{N} = f_{N} + β / h^{2}$ . The matrix $A_{h} \in R^{N \times N}$ is symmetric, tridiagonal, has diagonal entries $2/ h^{2} + q_{j} > 0$ and off-diagonal entries $- 1/ h^{2} < 0$ , and is solved by the Thomas algorithm — Gaussian elimination 43.03.01 specialised to a tridiagonal matrix, costing $O (N)$ operations with no fill-in. For the variable-coefficient form, the conservative discretisation evaluates $κ$ at the cell midpoints $x_{j \pm 1/2} = x_{j} \pm h /2$ : $$ -\big(\kappa(x)u'\big)'\big|{x_j} \approx -\frac{\kappa{j+1/2}(U_{j+1}-U_j) - \kappa_{j-1/2}(U_j - U_{j-1})}{h^2}, $$ which keeps $A_{h}$ symmetric tridiagonal with off-diagonals $- κ_{j \pm 1/2} / h^{2}$ and preserves the sign pattern.

Local truncation error. The local truncation error is the residual obtained by inserting the exact solution into the discrete operator, $$ \tau_j = \Big[-D^2 u(x_j) + q_j u(x_j)\Big] - f_j = -D^2 u(x_j) + u''(x_j), $$ measuring the failure of the exact solution to satisfy the discrete equation. The scheme is consistent of order $p$ if $∥ τ ∥_{\infty} = O (h^{p})$ . The symbols $D^{2}$ , $A_{h}$ , the tridiagonal generator $tridiag (- 1, 2, - 1)$ , the discrete maximum norm $∥ \cdot ∥_{\infty}$ , and the Landau symbol $O (h^{2})$ are introduced here with these definitions.

M-matrix and monotonicity. A real square matrix is an M-matrix if it has nonpositive off-diagonal entries and a nonnegative inverse ( $A_{h}^{- 1} \geq 0$ entrywise); equivalently here, $A_{h}$ is monotone — $A_{h} V \geq 0$ implies $V \geq 0$ entrywise. The sign pattern of $A_{h}$ (positive diagonal, nonpositive off-diagonals) together with weak diagonal dominance ( $2/ h^{2} + q_{j} \geq 2/ h^{2}$ , with strict dominance in the boundary-adjacent rows once the data is moved to $F$ ) and irreducibility delivers this structure, the algebraic carrier of the discrete maximum principle.

Counterexamples to common slips

Consistency is not convergence. The truncation error $τ$ measures how well the exact solution satisfies the discrete equation; alone it says nothing about how far the computed $U$ lies from $u$ . Bridging the two needs the stability bound on $A_{h}^{- 1}$ . A scheme can be consistent yet diverge if its discrete operator is not stable.
The naive one-sided Neumann difference is only first order. Imposing a Neumann condition $u^{'} (a) = γ$ by the one-sided difference $(U_{1} - U_{0}) / h = γ$ has truncation error $O (h)$ , which pollutes the global error down to first order even though the interior stencil is $O (h^{2})$ . The centred ghost-point treatment restores second order; conflating the two loses an order of accuracy.
The shooting residual can be ill-conditioned even when the BVP is well-posed. For problems with growing modes, the initial-value map $s \mapsto u (b; s)$ amplifies the slope guess exponentially, so the residual $ϕ (s)$ is steep and a small error in $s$ swings $u (b; s)$ wildly. The BVP is well-conditioned; the IVP reformulation need not be. Multiple shooting is the remedy.
The boundary values are data, not unknowns. Only the $N$ interior values are solved for; the boundary grid values enter $A_{h} U = F$ through the right-hand side $F$ . Treating boundary points as unknowns mis-sizes the system and breaks the M-matrix structure.

Key theorem with proof Intermediate+

The signature result is that the centred finite-difference scheme for the two-point BVP converges at second order in the maximum norm, and that this follows from two separately checkable facts: a consistency estimate (the truncation error is $O (h^{2})$ ) and a stability estimate (the discrete solution operator $A_{h}^{- 1}$ is bounded in the maximum norm uniformly in $h$ ). The stability estimate is itself a consequence of the discrete maximum principle, and the sharp constant is geometric.

Theorem (consistency + stability $\Rightarrow$ $O (h^{2})$ convergence). Let $u \in C^{4} [a, b]$ solve $- u^{''} + q u = f$ on $(a, b)$ with $q \geq 0$ and Dirichlet data $u (a) = α$ , $u (b) = β$ , and let $U$ solve the tridiagonal system $A_{h} U = F$ . Then there is a constant $C$ , independent of $h$ , with $$ |U - u|\infty := \max{1\le j\le N} |U_j - u(x_j)| \le C,h^2 \max_{[a,b]} |u''''|, \qquad C = \frac{(b-a)^2}{96}. $$ ^{[LeVeque, R. J. — Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems; Süli, E. & Mayers, D. F. — An Introduction to Numerical Analysis]}

Proof. Let $E_{j} = U_{j} - u (x_{j})$ be the global error grid function, vanishing at $j = 0, N + 1$ . Since $A_{h} U = F$ encodes $- D^{2} U_{j} + q_{j} U_{j} = f_{j}$ , and the exact solution satisfies $- D^{2} u (x_{j}) + q_{j} u (x_{j}) = f_{j} + τ_{j}$ by the definition of the truncation error, subtraction gives $- D^{2} E_{j} + q_{j} E_{j} = - τ_{j}$ at every interior point, that is $A_{h} E = - τ$ with $E$ vanishing on the boundary.

Consistency. Expand $u$ in a Taylor series about $x_{j}$ . The centred second difference gives $$ \frac{u(x_{j-1}) - 2u(x_j) + u(x_{j+1})}{h^2} = u''(x_j) + \frac{h^2}{12}u''''(\xi_j), $$ for some $ξ_{j} \in (x_{j - 1}, x_{j + 1})$ , the odd-order terms cancelling by symmetry. Hence $τ_{j} = - D^{2} u (x_{j}) + u^{''} (x_{j}) = - \frac{h ^{2}}{12} u^{''''} (ξ_{j})$ , so $∥ τ ∥_{\infty} \leq \frac{h ^{2}}{12} max_{[a, b]} ∣ u^{''''} ∣$ .

Stability. The claim is the uniform bound $∥ A_{h}^{- 1} ∥_{\infty} \leq (b - a)^{2} /8$ . Introduce the comparison (barrier) grid function $ϕ_{j} = \frac{1}{2} (x_{j} - a) (b - x_{j})$ , the restriction of a quadratic with $- ϕ^{''} = 1$ , hence $- D^{2} ϕ_{j} = 1$ exactly at every interior point (the centred second difference of a quadratic is its exact second derivative). On $[a, b]$ , $0 \leq ϕ_{j} \leq (b - a)^{2} /8$ , the maximum attained at the midpoint. Given any interior grid function $W$ with $A_{h} W = R$ and $W = 0$ on the boundary, set $W^{\pm} = \pm W - ∥ R ∥_{\infty} ϕ$ . Using $q \geq 0$ , at every interior point $$ (A_h W^{\pm})j = \pm R_j - |R|\infty(-D^2\phi_j + q_j\phi_j) = \pm R_j - |R|\infty(1 + q_j\phi_j) \le \pm R_j - |R|\infty \le 0, $$ since $q_{j} ϕ_{j} \geq 0$ and each component of $\pm R$ is at most $∥ R ∥_{\infty}$ ; on the boundary $W^{\pm} = - ∥ R ∥_{\infty} ϕ \leq 0$ . The discrete maximum principle (Proposition 1 below) gives $W^{\pm} \leq 0$ throughout, that is $\pm W_{j} \leq ∥ R ∥_{\infty} ϕ_{j} \leq \frac{( b - a ) ^{2}}{8} ∥ R ∥_{\infty}$ . Hence $∥ W ∥_{\infty} \leq \frac{( b - a ) ^{2}}{8} ∥ A_{h} W ∥_{\infty}$ , which is exactly $∥ A_{h}^{- 1} ∥_{\infty} \leq (b - a)^{2} /8$ .

Synthesis. Apply the stability bound to $E$ , which satisfies $A_{h} E = - τ$ and vanishes on the boundary: $$ |E|\infty \le \frac{(b-a)^2}{8}|A_h E|\infty = \frac{(b-a)^2}{8}|\tau|\infty \le \frac{(b-a)^2}{8}\cdot\frac{h^2}{12}\max{[a,b]}|u''''| = \frac{(b-a)^2}{96}h^2\max_{[a,b]}|u''''|. $$ This is the stated estimate with $C = (b - a)^{2} /96$ . $□$

Bridge. This proof is the foundational reason the difference method is trusted: it factors the error into a consistency piece, controlled by Taylor expansion of the stencil, and a stability piece, controlled by the maximum-norm bound on $A_{h}^{- 1}$ , and their product is the convergence rate. This is exactly the consistency-plus-stability template that builds toward the two-dimensional elliptic theory of 43.11.01, where the same comparison-function argument bounds the inverse of the 5-point Laplacian and the very tridiagonal generator $tridiag (- 1, 2, - 1)$ here reappears as the Kronecker factor $T$ assembled into $A_{h} = I \otimes T + T \otimes I$ ; it generalises the discrete-Gronwall stability idea of the IVP theory 43.10.01, where the role here played by $A_{h}^{- 1}$ was played by the accumulated one-step amplification factor. The stability half is dual to the consistency half — one bounds how the scheme amplifies residuals, the other bounds the residual the exact solution leaves behind — and putting these together is the central insight that a consistent, stable discretisation converges at the consistency rate. The same pattern appears again in the shooting method, where the stability question migrates from $A_{h}^{- 1}$ to the conditioning of the initial-value map $s \mapsto u (b; s)$ that the Newton root-solve inverts.

Exercises Intermediate+

Exercise 3 (medium, symbolic).

Show by Taylor expansion that $D^{2} u (x_{j}) = u^{''} (x_{j}) + \frac{h ^{2}}{12} u^{''''} (x_{j}) + O (h^{4})$ , so the second-difference operator is second-order accurate.

Hint

Expand $u (x_{j} \pm h)$ to fourth order and add; the odd-order terms cancel.

Answer

Expand $u (x_{j} + h) = u + h u^{'} + \frac{h ^{2}}{2} u^{''} + \frac{h ^{3}}{6} u^{'''} + \frac{h ^{4}}{24} u^{''''} + O (h^{5})$ and $u (x_{j} - h) = u - h u^{'} + \frac{h ^{2}}{2} u^{''} - \frac{h ^{3}}{6} u^{'''} + \frac{h ^{4}}{24} u^{''''} + O (h^{5})$ , all derivatives at $x_{j}$ . Adding, the odd-order terms cancel: $u (x_{j} + h) + u (x_{j} - h) = 2 u + h^{2} u^{''} + \frac{h ^{4}}{12} u^{''''} + O (h^{6})$ . Subtract $2 u$ and divide by $h^{2}$ : $D^{2} u (x_{j}) = u^{''} + \frac{h ^{2}}{12} u^{''''} + O (h^{4})$ . The leading error is the $\frac{h ^{2}}{12} u^{''''}$ term, so the operator is second-order. Rubric: full credit for both expansions, the cancellation, and the coefficient $\frac{1}{12}$ .

Exercise 4 (medium, symbolic).

Show that the linear shooting method is exact in one step. Let $- u^{''} + q u = f$ on $[a, b]$ with $u (a) = α$ , $u (b) = β$ , and let $u (\cdot; s)$ be the IVP solution with $u (a; s) = α$ , $u^{'} (a; s) = s$ . Prove that the shooting residual $ϕ (s) = u (b; s) - β$ is affine in $s$ , so a single secant/Newton step from any two guesses lands the exact slope.

Hint

Write $u (\cdot; s) = u_{p} + s w$ where $u_{p}$ solves the IVP with slope $0$ and $w$ solves the homogeneous IVP with slope $1$ , using linearity of the ODE.

Answer

By linearity, $u (\cdot; s) = u_{p} + s w$ , where $u_{p}$ solves $- u_{p}^{''} + q u_{p} = f$ with $u_{p} (a) = α$ , $u_{p}^{'} (a) = 0$ , and $w$ solves the homogeneous equation $- w^{''} + q w = 0$ with $w (a) = 0$ , $w^{'} (a) = 1$ . Then $u (b; s) = u_{p} (b) + s w (b)$ , so $ϕ (s) = u_{p} (b) + s w (b) - β$ , which is affine in $s$ with slope $w (b)$ . Solving $ϕ (s) = 0$ gives $s^{⋆} = (β - u_{p} (b)) / w (b)$ exactly, provided $w (b) \neq = 0$ (which holds when the BVP is uniquely solvable). A secant step through two guesses $s_{0}, s_{1}$ reproduces $s^{⋆}$ because the secant of an affine function is the function itself. Rubric: full credit for the superposition decomposition, the affine residual, and the exact one-step conclusion.

Exercise 5 (medium, numeric).

On $[0, 1]$ with $h = 1/4$ (three interior points), solve $- u^{''} = 8$ with $u (0) = u (1) = 0$ by the finite-difference method, and compare the centre value $U_{2}$ with the exact $u (1/2)$ .

Hint

The discrete equations are $- U_{j - 1} + 2 U_{j} - U_{j + 1} = 8 h^{2} = 1/2$ ; use the symmetry $U_{1} = U_{3}$ .

Answer

With $h^{2} = 1/16$ , each interior equation is $- U_{j - 1} + 2 U_{j} - U_{j + 1} = 1/2$ . By symmetry $U_{1} = U_{3}$ ; the middle equation gives $- 2 U_{1} + 2 U_{2} = 1/2$ so $U_{2} = U_{1} + 1/4$ , and the first gives $2 U_{1} - U_{2} = 1/2$ , hence $U_{1} = 3/4$ and $U_{2} = 1$ . The exact solution $u (x) = 4 x (1 - x)$ has $u (1/2) = 1$ , so $U_{2} = 1$ matches exactly. The match is exact here because a constant source makes $u$ quadratic, and the centred difference is exact on quadratics — the truncation error $\frac{h ^{2}}{12} u^{''''}$ vanishes when $u^{''''} \equiv 0$ .

Exercise 6 (medium, numeric).

For $- u^{''} + q u = f$ with constant $q = π^{2}$ on $[0, 1]$ , the diagonal entries of $A_{h}$ are $2/ h^{2} + π^{2}$ . With $h = 1/4$ , compute the diagonal entry and confirm $A_{h}$ is strictly diagonally dominant in the interior rows.

Hint

Interior rows have one $+ (2/ h^{2} + q)$ diagonal and two $- 1/ h^{2}$ off-diagonals. Compare the diagonal magnitude with the sum of the off-diagonal magnitudes.

Answer

With $h = 1/4$ , $1/ h^{2} = 16$ , so the diagonal entry is $2 \times 16 + π^{2} = 32 + 9.8696 \dots = 41.87$ . The two off-diagonals each have magnitude $16$ , summing to $32$ . Since $41.87 > 32$ , the interior rows are strictly diagonally dominant; the dominance margin is exactly $q = π^{2}$ . When $q = 0$ the interior rows are only weakly dominant ( $32 = 32$ ), and strict dominance comes solely from the boundary-adjacent rows. Either way $A_{h}$ is a nonsingular M-matrix.

Exercise 7 (hard, symbolic).

Prove the discrete maximum principle for $A_{h}$ with $q \geq 0$ : if a grid function $V$ satisfies $(- D^{2} V_{j} + q_{j} V_{j}) \leq 0$ at every interior point, then $max_{0 \leq j \leq N + 1} V_{j}$ is attained on the boundary (or $V \leq 0$ throughout).

Hint

At an interior maximum with $V_{j} > 0$ , use $q_{j} V_{j} \geq 0$ to deduce $V_{j}$ is at most the average of its two neighbours.

Answer

The hypothesis $- D^{2} V_{j} + q_{j} V_{j} \leq 0$ rearranges to $2 V_{j} \leq V_{j - 1} + V_{j + 1} - h^{2} q_{j} V_{j}$ . Let $M = max_{j} V_{j}$ and suppose $M$ is attained at an interior point $j_{0}$ with $M > 0$ . Then $q_{j_{0}} V_{j_{0}} = q_{j_{0}} M \geq 0$ , so $2 V_{j_{0}} \leq V_{j_{0} - 1} + V_{j_{0} + 1} - h^{2} q_{j_{0}} M \leq V_{j_{0} - 1} + V_{j_{0} + 1}$ , i.e. $V_{j_{0}}$ is at most the average of its neighbours. Since each neighbour is $\leq M = V_{j_{0}}$ , the average is $\leq M$ , forcing equality: both neighbours equal $M$ and (if $q_{j_{0}} > 0$ ) also $h^{2} q_{j_{0}} M = 0$ , so $M = 0$ , contradicting $M > 0$ unless $q_{j_{0}} = 0$ . Where $q \equiv 0$ , equality propagates along the connected interior to a boundary-adjacent point, whose boundary neighbour must then equal $M$ ; hence $M$ is attained on the boundary. If $M \leq 0$ then $V \leq 0$ throughout. Rubric: full credit for the sub-mean rearrangement using $q \geq 0$ , the equality-propagation argument, and the boundary conclusion.

Exercise 8 (hard, symbolic).

Set up Newton's method for the nonlinear two-point BVP $- u^{''} = g (x, u)$ , $u (a) = α$ , $u (b) = β$ , discretised as $G_{j} (U) := - D^{2} U_{j} - g (x_{j}, U_{j}) = 0$ . Compute the Jacobian $\partial G / \partial U$ and show it is tridiagonal, so each Newton step is an $O (N)$ tridiagonal solve.

Hint

Differentiate $G_{j}$ with respect to each $U_{k}$ . The $- D^{2}$ part contributes the constant tridiagonal stencil; the $- g$ part contributes only a diagonal term.

Answer

Write $G_{j} (U) = (- U_{j - 1} + 2 U_{j} - U_{j + 1}) / h^{2} - g (x_{j}, U_{j})$ . Differentiating: $\partial G_{j} / \partial U_{j - 1} = - 1/ h^{2}$ , $\partial G_{j} / \partial U_{j + 1} = - 1/ h^{2}$ , and $\partial G_{j} / \partial U_{j} = 2/ h^{2} - g_{u} (x_{j}, U_{j})$ , with all other partials zero. The Jacobian is therefore tridiagonal: $$ J(U) = \frac{1}{h^2}\mathrm{tridiag}(-1,2,-1) - \mathrm{diag}\big(g_u(x_j, U_j)\big). $$ Newton's method 43.02.03 iterates $J (U^{(k)}) δ = - G (U^{(k)})$ , $U^{(k + 1)} = U^{(k)} + δ$ , each step a tridiagonal solve costing $O (N)$ by the Thomas algorithm 43.03.01. When $g_{u} \leq 0$ (so $- g$ is monotone increasing in $u$ ), $J$ keeps the M-matrix sign pattern and the linear convergence theory carries over; Newton then converges quadratically near a solution. Rubric: full credit for the correct tridiagonal Jacobian, the diagonal $g_{u}$ term, and the $O (N)$ -per-step observation.

Advanced results Master

Theorem 1 ( $A_{h}$ is a symmetric positive definite M-matrix). For $q \geq 0$ the matrix $A_{h} = \frac{1}{h ^{2}} tridiag (- 1, 2, - 1) + diag (q_{j})$ is symmetric, has positive diagonal entries $2/ h^{2} + q_{j}$ and nonpositive off-diagonal entries $- 1/ h^{2}$ , is irreducible (the grid path is connected) and weakly diagonally dominant with strict dominance in the boundary-adjacent rows, and is therefore a nonsingular M-matrix with $A_{h}^{- 1} \geq 0$ entrywise. It is moreover positive definite: the tridiagonal part $\frac{1}{h ^{2}} tridiag (- 1, 2, - 1)$ has eigenvalues $\frac{4}{h ^{2}} sin^{2} (\frac{k π h}{2 ( b - a ) ^{- 1}}) > 0$ and $diag (q_{j}) \geq 0$ , so the sum is positive definite. Monotonicity ( $A_{h} V \geq 0 \Rightarrow V \geq 0$ ) is equivalent to $A_{h}^{- 1} \geq 0$ and is the matrix form of the discrete maximum principle ^{[LeVeque, R. J. — Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems]}.

Theorem 2 (the discrete Green's function and the sharp stability constant). For the Dirichlet problem with $q \geq 0$ the inverse $A_{h}^{- 1}$ has nonnegative entries $G_{h} (x_{j}, x_{k}) \geq 0$ , the discrete Green's function, satisfying the row-sum identity $\sum_{k} h^{2} G_{h} (x_{j}, x_{k}) \leq ϕ (x_{j})$ for the barrier $ϕ (x) = \frac{1}{2} (x - a) (b - x)$ with $- ϕ^{''} = 1$ ; consequently $∥ A_{h}^{- 1} ∥_{\infty} \leq max_{j} ϕ (x_{j}) \leq (b - a)^{2} /8$ , a mesh-independent constant. For the pure problem $q \equiv 0$ this is sharp, and $G_{h}$ converges entrywise to the continuous Green's function $G (x, ξ) = (x_{<} - a) (b - x_{>}) / (b - a)$ of $- d^{2} / d x^{2}$ on $[a, b]$ . The maximum-norm convergence rate $O (h^{2})$ is the product of this $O (1)$ stability constant with the $O (h^{2})$ truncation error ^{[Süli, E. & Mayers, D. F. — An Introduction to Numerical Analysis]}.

Theorem 3 (Neumann/Robin conditions, ghost points, and the half-order pitfall). A Neumann condition $u^{'} (a) = γ$ imposed by the one-sided difference $(U_{1} - U_{0}) / h = γ$ has truncation error $\frac{h}{2} u^{''} (a) + O (h^{2}) = O (h)$ , and this $O (h)$ boundary residual degrades the global convergence to first order despite the $O (h^{2})$ interior stencil. The second-order remedy introduces a ghost point $x_{- 1} = a - h$ and a fictitious value $U_{- 1}$ , imposes the centred condition $(U_{1} - U_{- 1}) / (2 h) = γ$ , and applies the interior stencil at the boundary node $x_{0}$ itself; eliminating $U_{- 1}$ yields a modified first row that restores $∥ τ ∥_{\infty} = O (h^{2})$ , at the cost of breaking the symmetry of $A_{h}$ unless the equation is rescaled by $\frac{1}{2}$ in the Neumann row. A Robin condition $κ u^{'} (a) + σ u (a) = γ$ ( $σ \geq 0$ ) is handled identically by the ghost-point centred difference, and the resulting first row preserves the M-matrix sign pattern. The discrete operator on a pure-Neumann or periodic interval is singular (the constant grid function is a null vector), reflecting the continuous solvability constraint $\int f = 0$ ^{[LeVeque, R. J. — Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems]}.

Theorem 4 (nonlinear BVPs by Newton, and the contrast with shooting). For the nonlinear two-point BVP $- u^{''} = g (x, u, u^{'})$ with Dirichlet data, the centred discretisation $G_{j} (U) = - D^{2} U_{j} - g (x_{j}, U_{j}, D_{0} U_{j}) = 0$ (with $D_{0} U_{j} = (U_{j + 1} - U_{j - 1}) / (2 h)$ ) is a nonlinear system whose Newton iteration 43.02.03 has tridiagonal Jacobian $$ J(U) = \tfrac1{h^2}\mathrm{tridiag}(-1,2,-1) - \mathrm{diag}(g_u(x_j, U_j, D_0 U_j)) - \tfrac{1}{2h}\mathrm{antidiag}\big(g_{u'}\big), $$ solved in $O (N)$ per step; when $g_{u} \leq 0$ and $h ∣ g_{u^{'}} ∣ \leq 2$ the Jacobian remains an M-matrix and Newton converges quadratically from a good initial guess. The shooting alternative recasts the BVP as the IVP $u^{''} = - g (x, u, u^{'})$ , $u (a) = α$ , $u^{'} (a) = s$ , integrated by a one-step method 43.10.01, and solves the residual $ϕ (s) = u (b; s) - β = 0$ by Newton, with $ϕ^{'} (s)$ obtained from the variational (sensitivity) IVP $v^{''} = - (g_{u} v + g_{u^{'}} v^{'})$ , $v (a) = 0$ , $v^{'} (a) = 1$ . Shooting inherits the conditioning of the IVP: for problems with exponentially growing modes the map $s \mapsto u (b; s)$ is exponentially sensitive, and multiple shooting — partitioning $[a, b]$ and matching at interface nodes — is the stabilised variant that trades the ill-conditioned single trajectory for a larger but well-conditioned coupled system ^{[Keller, H. B. — Numerical Methods for Two-Point Boundary-Value Problems]}.

Synthesis. The one-dimensional two-point BVP is the model case in which consistency and stability are each visible as a separate, checkable structure, and convergence is their product: the foundational reason the difference method works is that $A_{h}$ is an M-matrix, and this single algebraic fact carries the discrete maximum principle, the inverse-monotonicity $A_{h}^{- 1} \geq 0$ , the sharp mesh-independent stability constant $(b - a)^{2} /8$ , and the nonnegativity of the discrete Green's function. This is exactly the discrete shadow of the continuous BVP, whose Green's function $G_{h}$ approximates and whose comparison/barrier argument the stability proof restricts to the grid. The central insight is that the maximum norm, not the energy norm, is the natural setting for the finite-difference method — the M-matrix monotonicity is an $L^{\infty}$ phenomenon — and this same apparatus is dual to the two-dimensional elliptic theory of 43.11.01, where the tridiagonal generator here becomes the Kronecker factor of the 5-point Laplacian and the quadratic barrier becomes the paraboloid comparison function. Putting these together, the consistency-plus-stability proof here generalises into the Lax-Richtmyer framework for evolution problems and specialises, in the shooting method, into a statement about the conditioning of an initial-value flow composed with a Newton root-solve; the bridge throughout is that local accuracy fixed by Taylor matching, propagated under a stability bound, yields global convergence at the consistency rate, whether the stability lives in $A_{h}^{- 1}$ or in the derivative of the shooting map.

Full proof set Master

Proposition 1 (discrete maximum principle, $q \geq 0$ ). If a grid function $V$ on ${x_{0}, \dots, x_{N + 1}}$ satisfies $(- D^{2} V_{j} + q_{j} V_{j}) \leq 0$ at every interior point with $q_{j} \geq 0$ , then $max_{0 \leq j \leq N + 1} V_{j} \leq max (0, max_{j \in {0, N + 1}} V_{j})$ ; in particular if $V \leq 0$ on the boundary then $V \leq 0$ throughout.

Proof. The hypothesis is $2 V_{j} \leq V_{j - 1} + V_{j + 1} - h^{2} q_{j} V_{j}$ . Let $M = max_{j} V_{j}$ and suppose $M > 0$ is attained at an interior index $j_{0}$ . Then $h^{2} q_{j_{0}} V_{j_{0}} = h^{2} q_{j_{0}} M \geq 0$ , so $2 V_{j_{0}} \leq V_{j_{0} - 1} + V_{j_{0} + 1} - h^{2} q_{j_{0}} M \leq V_{j_{0} - 1} + V_{j_{0} + 1}$ , hence $V_{j_{0}} \leq \frac{1}{2} (V_{j_{0} - 1} + V_{j_{0} + 1}) \leq M$ . Equality forces both neighbours to equal $M$ and $h^{2} q_{j_{0}} M = 0$ . If $q_{j_{0}} > 0$ this gives $M = 0$ , contradicting $M > 0$ . Where $q \equiv 0$ on a stretch of interior points, equality propagates along the connected path until it reaches an index adjacent to the boundary; that boundary value then also equals $M$ , so $M = max_{j \in {0, N + 1}} V_{j}$ . In all cases $M \leq max (0, max_{\partial} V)$ . The final clause is the case $max_{\partial} V \leq 0$ . $□$

Proposition 2 ( $A_{h}$ is monotone: $A_{h}^{- 1} \geq 0$ ). For $q \geq 0$ the matrix $A_{h}$ representing $- D^{2} + q$ on the interior with Dirichlet boundary is nonsingular and $A_{h}^{- 1} \geq 0$ entrywise; equivalently, $A_{h} V \geq 0$ implies $V \geq 0$ .

Proof. Nonsingularity: if $A_{h} V = 0$ with $V = 0$ on the boundary, then $V$ satisfies both $- D^{2} V + q V \leq 0$ and $\geq 0$ , so by Proposition 1 applied to $V$ and to $- V$ both $max V$ and $min V$ are controlled by the boundary values $0$ ; hence $V \equiv 0$ and $A_{h}$ is injective, thus invertible. For monotonicity, suppose $A_{h} V \geq 0$ , i.e. $- D^{2} V_{j} + q_{j} V_{j} \geq 0$ , with $V = 0$ on the boundary. Apply Proposition 1 to $- V$ : $- V$ satisfies $- D^{2} (- V) + q (- V) \leq 0$ , so $max (- V) \leq max (0, max_{\partial} (- V)) = 0$ , giving $V \geq 0$ everywhere. Taking $V = A_{h}^{- 1} e_{k}$ for each basis vector $e_{k} \geq 0$ shows every column of $A_{h}^{- 1}$ is nonnegative, hence $A_{h}^{- 1} \geq 0$ . The sign pattern (positive diagonal, nonpositive off-diagonals) together with this nonnegative inverse is the definition of an M-matrix. $□$

Proposition 3 (uniform stability bound $∥ A_{h}^{- 1} ∥_{\infty} \leq (b - a)^{2} /8$ ). For every $h = (b - a) / (N + 1)$ and $q \geq 0$ , the discrete solution operator satisfies $∥ A_{h}^{- 1} ∥_{\infty} \leq (b - a)^{2} /8$ , independent of $h$ .

Proof. Let $ϕ$ be the grid restriction of $Φ (x) = \frac{1}{2} (x - a) (b - x)$ , for which $Φ^{''} = - 1$ and, since $Φ$ is quadratic and the centred second difference is exact on quadratics, $- D^{2} ϕ_{j} = - Φ^{''} (x_{j}) = 1$ at every interior point. On $[a, b]$ , $0 \leq Φ (x) \leq (b - a)^{2} /8$ with the maximum at the midpoint. Fix a right-hand side $R$ and let $W = A_{h}^{- 1} R$ , so $W = 0$ on the boundary. Set $Z^{\pm} = \pm W - ∥ R ∥_{\infty} ϕ$ . At interior points, using $q_{j} \geq 0$ and $q_{j} ϕ_{j} \geq 0$ , $$ (-D^2 Z^\pm_j + q_j Z^\pm_j) = \pm R_j - |R|\infty(-D^2\phi_j + q_j\phi_j) = \pm R_j - |R|\infty(1 + q_j\phi_j) \le \pm R_j - |R|\infty \le 0, $$ since $|R_j| \le |R|\infty $. O n t h e b o u n d a r y$ Z^\pm = -|R|\infty\phi \le 0 $. B y P r o p os i t i o n 1,$ Z^\pm \le 0 $t h r o ug h o u t, so$ \pm W_j \le |R|\infty\phi_j \le \tfrac{(b-a)^2}{8}|R|\infty $f or a l l$ j $. H e n ce$ |W|\infty \le \tfrac{(b-a)^2}{8}|R|\infty $f or e v er y$ R $, w hi c hi s$ |A_h^{-1}|\infty \le (b-a)^2/8 $.$ \square$

Proposition 4 (exact eigensystem of the Dirichlet second-difference matrix). On $[a, b]$ with $h = (b - a) / (N + 1)$ , the matrix $T_{h} = \frac{1}{h ^{2}} tridiag (- 1, 2, - 1) \in R^{N \times N}$ has the orthogonal eigenvectors $v_{j}^{(k)} = sin (\frac{k π ( x _{j} - a )}{b - a})$ , $1 \leq k \leq N$ , with eigenvalues $μ_{k} = \frac{4}{h ^{2}} sin^{2} (\frac{k π h}{2 ( b - a )})$ .

Proof. Write $θ_{k} = k π h / (b - a)$ and $s_{j}^{(k)} = sin (j θ_{k})$ , noting $x_{j} - a = j h$ so $v_{j}^{(k)} = s_{j}^{(k)}$ . Apply $T_{h}$ at interior index $j$ : $$ (T_h s^{(k)})_j = \frac{-\sin((j-1)\theta_k) + 2\sin(j\theta_k) - \sin((j+1)\theta_k)}{h^2}. $$ The identity $sin ((j - 1) θ) + sin ((j + 1) θ) = 2 sin (j θ) cos θ$ gives $$ (T_h s^{(k)})_j = \frac{2\sin(j\theta_k)(1 - \cos\theta_k)}{h^2} = \frac{2(1-\cos\theta_k)}{h^2},s^{(k)}_j = \frac{4}{h^2}\sin^2!\tfrac{\theta_k}{2},s^{(k)}_j, $$ the last step by the half-angle identity. The boundary entries vanish: $s_{0}^{(k)} = sin 0 = 0$ and $s_{N + 1}^{(k)} = sin ((N + 1) θ_{k}) = sin (k π) = 0$ , so $s^{(k)}$ is a genuine interior eigenvector with $μ_{k} = \frac{4}{h ^{2}} sin^{2} (θ_{k} /2)$ . The $N$ vectors $s^{(1)}, \dots, s^{(N)}$ are orthogonal (discrete sine transform) and complete in $R^{N}$ . Adding $diag (q_{j})$ shifts the quadratic form by a nonnegative amount, so $A_{h} = T_{h} + diag (q_{j})$ is positive definite with smallest eigenvalue $\geq μ_{1} = \frac{4}{h ^{2}} sin^{2} (\frac{π h}{2 ( b - a )}) \to π^{2} / (b - a)^{2}$ as $h \to 0$ . $□$

Connections Master

The one-step ODE integrators of 43.10.01 are the engine inside the shooting method: shooting recasts the two-point BVP as the initial-value problem $u^{''} = - g (x, u, u^{'})$ with an unknown initial slope $s$ , marches it to the far endpoint with an Euler/Runge-Kutta scheme of known order, and reads the boundary residual $ϕ (s) = u (b; s) - β$ . The accuracy of the integrator sets the accuracy of $ϕ$ , and the IVP convergence theory of that unit is what guarantees the shooting trajectory itself is computed to order $h^{p}$ ; this unit supplies the boundary-coupling and root-solve layer that wraps around the IVP march.
The Newton and secant root-finders of 43.02.03 solve the two nonlinear equations this unit produces: the scalar shooting residual $ϕ (s) = 0$ (with $ϕ^{'}$ from the variational IVP, or one exact secant step in the linear case), and the discrete nonlinear system $G (U) = 0$ from a nonlinear BVP, whose tridiagonal Jacobian makes each Newton step an $O (N)$ solve. The superlinear/quadratic convergence theory there is what certifies that the outer iteration of both methods terminates rapidly near a solution.
The Gaussian-elimination and LU theory of 43.03.01 is what actually solves the tridiagonal system $A_{h} U = F$ at the heart of the finite-difference method: specialised to a tridiagonal matrix, elimination becomes the Thomas algorithm, costing $O (N)$ operations with no fill-in, and the symmetric positive definiteness of $A_{h}$ established here means the factorisation needs no pivoting. Each Newton step of a nonlinear BVP is one such tridiagonal solve, so the linear-solver unit is invoked once per outer iteration.
The two-dimensional elliptic finite-difference theory of 43.11.01 is the direct multidimensional successor: the tridiagonal generator $tridiag (- 1, 2, - 1)$ that is the whole operator here becomes the Kronecker factor $T$ assembled into the 5-point Laplacian $A_{h} = I \otimes T + T \otimes I$ there, the quadratic barrier $\frac{1}{2} (x - a) (b - x)$ becomes the paraboloid comparison function on the square, and the one-dimensional discrete maximum principle becomes the two-dimensional one. This unit is the genuinely one-dimensional ODE-BVP precursor of that elliptic PDE unit; it owns the scalar two-point problem, shooting, variable coefficients, and the ghost-point boundary conditions, while that unit owns the Kronecker assembly and the explicit two-dimensional spectrum.
The continuous nth-order linear ODE theory of 02.06.02 is the object this unit discretises and the source of the superposition principle that makes linear shooting exact in one step: the decomposition $u (\cdot; s) = u_{p} + s w$ into a particular and a homogeneous solution is exactly the structure of the general solution of a linear ODE, and the well-posedness of the BVP (uniqueness when $w (b) \neq = 0$ ) rests on the continuous boundary-value theory the numerical method approximates without reproving.

Historical & philosophical context Master

The finite-difference treatment of boundary-value problems for ordinary differential equations descends from the same 1928 work of Richard Courant, Kurt Friedrichs, and Hans Lewy, Über die partiellen Differenzengleichungen der mathematischen Physik (Mathematische Annalen 100, 32–74) ^{[Courant, Friedrichs & Lewy 1928]}, that founded the convergence theory of difference schemes; their use of the discrete maximum principle and a comparison function to bound the discrete inverse is the argument that, specialised to one dimension, gives the stability constant $(b - a)^{2} /8$ used here. The recasting of convergence into the algebraic language of monotone and M-matrices — the inverse-positivity property $A^{- 1} \geq 0$ that carries the discrete maximum principle — is due to the matrix-theoretic tradition of Lothar Collatz, whose The Numerical Treatment of Differential Equations (Springer, 1960) developed the discrete maximum principle for two-point problems, and Richard Varga, whose Matrix Iterative Analysis (Prentice-Hall, 1962) connected diagonal dominance and irreducibility to monotonicity.

The shooting method and its multiple-shooting stabilisation, together with the systematic finite-difference theory for linear and nonlinear two-point boundary-value problems, were consolidated by Herbert B. Keller in Numerical Methods for Two-Point Boundary-Value Problems (Blaisdell, 1968; reprinted Dover, 1992) ^{[Keller, H. B. — Numerical Methods for Two-Point Boundary-Value Problems]}, which gave the conditioning analysis explaining why single shooting fails for problems with growing modes and how matching across subintervals repairs it. The variable-coefficient conservative discretisation and the ghost-point treatment of Neumann and Robin conditions, with the attendant half-order accuracy caveat, are developed in Randall LeVeque's Finite Difference Methods for Ordinary and Partial Differential Equations (SIAM, 2007) ^{[LeVeque, R. J. — Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems]}, whose Chapter 2 is the steady-state one-dimensional core that this unit follows.

Bibliography Master

@book{leveque2007fdm,
  author    = {LeVeque, Randall J.},
  title     = {Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems},
  publisher = {Society for Industrial and Applied Mathematics (SIAM)},
  year      = {2007}
}

@book{sulimayers2003,
  author    = {S\"{u}li, Endre and Mayers, David F.},
  title     = {An Introduction to Numerical Analysis},
  publisher = {Cambridge University Press},
  year      = {2003}
}

@book{keller1992twopoint,
  author    = {Keller, Herbert B.},
  title     = {Numerical Methods for Two-Point Boundary-Value Problems},
  publisher = {Dover Publications},
  year      = {1992},
  note      = {Reprint of the Blaisdell edition, 1968}
}

@book{stoerbulirsch2002,
  author    = {Stoer, Josef and Bulirsch, Roland},
  title     = {Introduction to Numerical Analysis},
  edition   = {3},
  series    = {Texts in Applied Mathematics},
  volume    = {12},
  publisher = {Springer-Verlag},
  year      = {2002}
}

@book{collatz1960differential,
  author    = {Collatz, Lothar},
  title     = {The Numerical Treatment of Differential Equations},
  edition   = {3},
  publisher = {Springer-Verlag},
  year      = {1960}
}

@book{varga2000matrix,
  author    = {Varga, Richard S.},
  title     = {Matrix Iterative Analysis},
  edition   = {2},
  publisher = {Springer-Verlag},
  year      = {2000}
}

@article{courantfriedrichslewy1928,
  author  = {Courant, Richard and Friedrichs, Kurt and Lewy, Hans},
  title   = {{\"U}ber die partiellen Differenzengleichungen der mathematischen Physik},
  journal = {Mathematische Annalen},
  volume  = {100},
  number  = {1},
  year    = {1928},
  pages   = {32--74}
}

Prerequisites

43.10.01
43.02.03
43.03.01

Tier anchors

beginner: LeVeque 2007 *Finite Difference Methods for Ordinary and Partial Differential Equations* (SIAM) §2.1-2.4 (the steady-state two-point problem $u''=f$ on an interval, the second-difference stencil, and the tridiagonal system as a picture); Süli-Mayers 2003 *An Introduction to Numerical Analysis* (Cambridge) §13.1-13.2 (shooting as guess-the-slope-and-correct, and the difference replacement of the BVP)
intermediate: LeVeque 2007 *Finite Difference Methods for Ordinary and Partial Differential Equations* (SIAM) §2.4-2.12 (the centred discretisation, the $O(h^2)$ truncation error, the tridiagonal $A_h$, stability via $\|A_h^{-1}\|_\infty$, Neumann/ghost-point handling, and nonlinear BVPs by Newton); Süli-Mayers 2003 (Cambridge) §13.1-13.6 (the shooting method, the finite-difference method, the discrete maximum principle and the convergence theorem)
master: LeVeque 2007 *Finite Difference Methods for Ordinary and Partial Differential Equations* (SIAM) Ch. 2 (the steady-state BVP in full: variable coefficients $(\kappa u')'$, the discrete Green's function and comparison-function stability bound, the half-order Neumann pitfall, and Newton on the discrete nonlinear system); Keller 1992 *Numerical Methods for Two-Point Boundary-Value Problems* (Dover) Ch. 1-3 (shooting, multiple shooting, and the finite-difference method with its convergence theory); Stoer-Bulirsch 2002 *Introduction to Numerical Analysis* 3e (Springer) §7.3 (shooting, multiple shooting, and difference methods for boundary-value problems)

References

LeVeque, R. J. — Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems · SIAM 2007. Chapter 2. §2.1-2.3 pose the one-dimensional steady-state two-point BVP $u''=f$ (and the variable-coefficient $(\kappa u')'=f$) on $[a,b]$ with Dirichlet data; §2.4 builds the centred second-difference discretisation, writes the symmetric tridiagonal system $A_h U=F$, and gives the local truncation error $\tau=\tfrac{h^2}{12}u''''+O(h^4)$. §2.5-2.7 prove the $O(h^2)$ global convergence $\|E\|_\infty=O(h^2)$ via consistency plus the stability bound $\|A_h^{-1}\|_\infty\le(b-a)^2/8$ from a discrete comparison/Green's-function argument and the discrete maximum principle. §2.11-2.12 treat Neumann/Robin boundary conditions by one-sided and ghost-point stencils, flag the first-order accuracy of the naive one-sided Neumann difference, and develop nonlinear two-point BVPs $u''=g(x,u,u')$ solved by Newton's method on the discrete system with a tridiagonal Jacobian.
Süli, E. & Mayers, D. F. — An Introduction to Numerical Analysis · Cambridge University Press 2003. Chapter 13. §13.1-13.2 introduce the two-point BVP $-u''+q(x)u=f$ with $q\ge0$ and the shooting method: recast the BVP as an IVP with unknown initial slope $s$, integrate with a Chapter-12 one-step method, and solve the residual $\phi(s)=u(b;s)-\beta=0$ for $s$ by Newton/secant, noting that for the linear BVP a single secant/Newton step is exact by superposition. §13.3-13.5 develop the finite-difference method: replace $u''$ by $(U_{j-1}-2U_j+U_{j+1})/h^2$ to obtain the symmetric tridiagonal $A_h U=F$, prove the discrete maximum principle and the monotonicity of $A_h$, establish $\|A_h^{-1}\|_\infty\le(b-a)^2/8$, and conclude $O(h^2)$ convergence from consistency plus stability.
Keller, H. B. — Numerical Methods for Two-Point Boundary-Value Problems · Dover 1992 (reprint of Blaisdell 1968). Chapters 1-3. The simple and multiple shooting methods and their sensitivity to the conditioning of the underlying IVP; the finite-difference method for linear and nonlinear two-point BVPs; the discrete maximum principle; and the convergence theorems relating the consistency order of the difference stencil to the global error under a stability (inverse-boundedness) hypothesis on the discrete operator.

Estimated time

beginner: 20m
intermediate: 50m
master: 90m