Boundary value problems

Lessons 7 and 8 solved initial value problems, where everything is known at the start and the solver marches forward. Many engineering problems instead fix conditions at two ends: the temperature at the base of a fin and the heat flow at its tip, the supports at both ends of a beam, the launch point and the target of a projectile, the concentrations at the two faces of a membrane. These boundary value problems (BVPs) cannot be marched directly, because the missing starting values are unknown, and MATH-201 showed that they may have one solution, several or none. This lesson develops the two main numerical approaches. The shooting method turns a BVP into a root-finding problem (Lesson 2) wrapped around an IVP solver. Finite difference methods replace derivatives by differences on a grid and solve a system of equations (Lesson 3), using Newton's method when the problem is nonlinear. The lesson also shows how a BVP with an unknown parameter becomes a matrix eigenvalue problem (Lesson 4), the numerical route to buckling loads and mode shapes.

The shooting method

For a second-order BVP such as \(y'' = f(x, y, y')\) with \(y(a) = \alpha\) and \(y(b) = \beta\), guess the missing initial slope s, solve the IVP from a to b with \(y(a) = \alpha\), \(y'(a) = s\), and see how far \(y(b)\) misses β. The miss distance \(F(s) = y(b; s) - \beta\) is a function of one variable, and the BVP is solved when \(F(s) = 0\). Any root finder will do: bisection if two guesses bracket the target, the secant method otherwise. For a linear equation, F is linear in s, so two shots and one interpolation give the exact answer.

Shooting is intuitive and uses accurate adaptive IVP solvers, and it handles nonlinear problems naturally. It has two weaknesses. If the IVP is sensitive, with solutions that grow rapidly (as in some problems with exponential behavior over long intervals), small changes in s produce enormous changes at b and the root finder struggles; multiple shooting, which breaks the interval into pieces joined by matching conditions, cures this. And it does not extend easily to problems in two or three dimensions.

import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import brentq

# Ball launched at 30 m/s with quadratic drag; find the launch angle that lands 45 m away
g, k, v0, target = 9.81, 0.01, 30.0, 45.0     # k = drag constant per unit mass (1/m)

def rhs(t, s):
    x, y, vx, vy = s
    v = np.hypot(vx, vy)
    return [vx, vy, -k * v * vx, -g - k * v * vy]

def landing(t, s):
    return s[1]
landing.terminal, landing.direction = True, -1

def miss(angle_deg):
    a = np.radians(angle_deg)
    sol = solve_ivp(rhs, [0, 20], [0, 0, v0 * np.cos(a), v0 * np.sin(a)], events=landing,
                    rtol=1e-9, atol=1e-9)
    return sol.y_events[0][0][0] - target

for a in [10, 20, 30, 45, 60, 75]:
    print(f"angle {a:2d} deg: lands {miss(a) + target:6.2f} m")
low = brentq(miss, 5, 38)
high = brentq(miss, 42, 85)
print(f"Shooting solutions: {low:.3f} deg (flat) and {high:.3f} deg (lob)")

angle 10 deg: lands  26.13 m
angle 20 deg: lands  42.96 m
angle 30 deg: lands  52.61 m
angle 45 deg: lands  55.78 m
angle 60 deg: lands  46.90 m
angle 75 deg: lands  27.36 m
Shooting solutions: 21.645 deg (flat) and 61.882 deg (lob)

Finite difference methods

The finite difference method treats the whole interval at once. Divide \([a, b]\) into N intervals of width h, with unknown values \(y_i\) at the grid points, and replace each derivative by a centered difference (Lesson 1):

Both approximations have error proportional to \(h^2\). Writing the differential equation at every interior point gives one equation per unknown. For a linear second-order equation, each equation involves only three neighboring unknowns, so the system is tridiagonal and can be solved in \(O(N)\) operations (the Thomas algorithm, or scipy.linalg.solve_banded), even for millions of points.

Boundary conditions enter in two ways. A fixed value (a Dirichlet condition, such as a given base temperature) simply replaces an unknown with a known number, moved to the right side. A derivative condition (a Neumann condition, such as an insulated tip with zero heat flow) is handled with a ghost point beyond the boundary: the centered condition \((y_{N+1} - y_{N-1})/(2h) = 0\) gives \(y_{N+1} = y_{N-1}\), which is substituted into the equation at the boundary point, keeping second-order accuracy.

import numpy as np
from scipy.linalg import solve_banded

# Fin: theta double-prime = m^2 theta, theta(0) = 1, zero slope at x = L (insulated tip); exact cosh solution
m, L = 10.0, 0.05
exact_tip = 1 / np.cosh(m * L)

def fin_fd(N):
    h = L / N
    # unknowns theta_1 .. theta_N; tip uses a ghost point theta_{N+1} = theta_{N-1}
    main = np.full(N, -(2 + (m * h)**2))
    upper = np.ones(N); lower = np.ones(N)
    lower[N - 2] = 2.0                       # tip row: 2 theta_{N-1} - (2 + m^2 h^2) theta_N = 0
    rhs = np.zeros(N); rhs[0] = -1.0         # theta_0 = 1 moved to the right side
    ab = np.vstack([np.r_[0, upper[:-1]], main, np.r_[lower[:-1], 0]])
    return solve_banded((1, 1), ab, rhs)[-1]

prev = None
for N in [5, 10, 20, 40]:
    err = abs(fin_fd(N) - exact_tip)
    ratio = f"{prev / err:5.2f}" if prev else "   --"
    print(f"N = {N:3d}: tip value {fin_fd(N):.7f}, error {err:.2e}, ratio {ratio}")
    prev = err
print(f"Exact tip value: {exact_tip:.7f}")

N =   5: tip value 0.8869042, error 8.53e-05, ratio    --
N =  10: tip value 0.8868402, error 2.13e-05, ratio  4.00
N =  20: tip value 0.8868242, error 5.34e-06, ratio  4.00
N =  40: tip value 0.8868202, error 1.33e-06, ratio  4.00
Exact tip value: 0.8868189

Nonlinear BVPs and Newton's method

When the equation is nonlinear, the finite difference equations form a nonlinear system \(\mathbf{F}(\mathbf{y}) = \mathbf{0}\), solved by Newton's method for systems (Lesson 2): at each iteration, solve \(J\,\Delta\mathbf{y} = -\mathbf{F}\) with the Jacobian matrix \(J = \partial\mathbf{F}/\partial\mathbf{y}\), which is tridiagonal here, and update. Newton converges quadratically from a good initial guess; a poor guess may fail or, as the next example shows, may converge to a different solution.

import numpy as np

# Bratu problem: y double-prime + lam * exp(y) = 0, y(0) = y(1) = 0 (thermal ignition in a slab)
N, lam = 50, 2.0
h = 1.0 / N
x = np.linspace(0, 1, N + 1)[1:-1]

def newton(y):
    for it in range(1, 50):
        F = np.zeros(N - 1)
        F[:] = -2 * y + lam * h**2 * np.exp(y)
        F[1:] += y[:-1]
        F[:-1] += y[1:]
        J = np.diag(-2 + lam * h**2 * np.exp(y)) + np.diag(np.ones(N - 2), 1) + np.diag(np.ones(N - 2), -1)
        dy = np.linalg.solve(J, -F)
        y = y + dy
        if np.max(np.abs(dy)) < 1e-12:
            return y, it
    return y, it

for name, guess in [("zero guess", np.zeros(N - 1)), ("large guess", 8 * x * (1 - x) * 4)]:
    y, its = newton(guess)
    print(f"{name:11s}: converged in {its} Newton steps, peak value y(0.5) = {y.max():.4f}")

zero guess : converged in 5 Newton steps, peak value y(0.5) = 0.3290
large guess: converged in 12 Newton steps, peak value y(0.5) = 2.8951

Eigenvalue problems from BVPs

A homogeneous BVP containing a parameter, such as the buckling equation \(EI\,y'' + Py = 0\) with pinned ends, has non-zero solutions only for particular parameter values (MATH-201). Discretizing it gives a matrix eigenvalue problem. With the centered second difference and \(y_0 = y_N = 0\):

so the buckling loads are the eigenvalues of the tridiagonal matrix on the left, and the eigenvectors are the buckled shapes sampled at the grid points. For a uniform column the smallest eigenvalue converges to the Euler load \(\pi^2EI/L^2\) with error proportional to \(h^2\). The real value of the method is for non-uniform problems: a tapered column, a pile in soil, a shaft with varying diameter, where EI varies along the length and no formula exists. Then \(EI_i\) simply varies from row to row, and the shifted inverse iteration of Lesson 4 finds the lowest loads. Vibration modes of strings, bars and beams are found the same way.

Library solvers

In practice, scipy.integrate.solve_bvp solves nonlinear BVPs by collocation: it represents the solution by piecewise cubic polynomials, requires the equation to hold at chosen points in each interval, and adapts the mesh until a tolerance is met. It needs the problem written as a first-order system, the boundary conditions as a residual function, and an initial guess for the solution on a mesh. A good guess, often from a simplified linear problem or physical intuition, matters as much here as for Newton's method, and the warnings about multiple solutions apply equally.