Prev
Unit 10 · Boundary Value Problems and Partial Differential Equations
Boundary Value Problems: Shooting and Finite Differences
60m100 XP
Next
🌡️Boundary Value Problems and Partial Differential Equations

Boundary Value Problems: Shooting and Finite Differences

Two-point BVPs, the shooting method, finite-difference discretisation, convergence and solve_bvp

Leiden, 1638
In his last book, Two New Sciences, Galileo asks how much load a beam can carry. He imagines a cantilever breaking at the wall and reasons about the resisting fibres of its cross-section — and gets the answer wrong by a factor of three, because he assumes every fibre pulls with the same force. It took more than a century, and the work of Jacob Bernoulli, Leonhard Euler and others, to arrive at the beam equation engineers still use: the bending moment is proportional to the curvature, and for a uniform beam under a distributed load qq, EI w′′′′(x)=q(x).EI\,w''''(x) = q(x). Four derivatives, so four conditions are needed. But they are not all given at one end. A beam clamped at the wall and resting on a support at the other has w=w′=0w = w' = 0 at the wall and w=w′′=0w = w'' = 0 (no deflection, no moment) at the support. A computer cannot simply start at one end and march — it does not know the slope and curvature it should leave the wall with. This is a boundary value problem: the conditions are split between two boundaries, and the solution must satisfy all of them at once.
A two-point boundary value problem is an ODE on [a,b][a, b] with conditions at both ends — fin temperatures with the base and the tip both specified, a cable sagging between two towers, the deflection of a beam, the steady concentration in a reactor with inlet and outlet conditions. Unlike initial value problems, BVPs may have no solution or many; for the linear problems of engineering there is usually exactly one. Two families of methods solve them.
Derivation
Where the fin and beam equations come from
The fin. Cut a thin slice of length dxdx from a fin of cross-section AA, perimeter PP and conductivity kk. Heat conducted in at one face minus heat conducted out at the other (Fourier's law, flux =−kA dT/dx= -kA\,dT/dx) must equal heat lost through the slice's surface to the air (Newton's law of cooling, hP dx (T−T∞)hP\,dx\,(T - T_\infty)) in steady state:
kA T′(x+dx)−kA T′(x)=hP dx (T−T∞)⟹T′′=m2(T−T∞),m2=hPkA.kA\,T'(x + dx) - kA\,T'(x) = hP\,dx\,(T - T_\infty) \quad\Longrightarrow\quad T'' = m^2(T - T_\infty), \qquad m^2 = \frac{hP}{kA}.
Its general solution is T−T∞=C1cosh⁡mx+C2sinh⁡mxT - T_\infty = C_1\cosh mx + C_2\sinh mx — the analytic answer the methods below are checked against, with C1,C2C_1, C_2 fixed by the two end temperatures.
The beam. Bernoulli and Euler's law says the bending moment is proportional to the curvature, M=EI w′′M = EI\,w'' (with ww the deflection). Equilibrium of a short slice relates the moment, the shear force VV and the distributed load qq: V=dM/dxV = dM/dx and q=dV/dxq = dV/dx (with the usual sign conventions). Chain the three:
q=d2Mdx2=d2dx2(EI w′′)=EI w′′′′for uniform EI.q = \frac{d^2M}{dx^2} = \frac{d^2}{dx^2}\big(EI\,w''\big) = EI\,w'''' \quad\text{for uniform } EI.
The boundary conditions come from the same quantities: a clamped end fixes w=0w = 0 and w′=0w' = 0; a pinned end w=0w = 0 and M=0M = 0, i.e. w′′=0w'' = 0; a free end M=0M = 0 and V=0V = 0, i.e. w′′=w′′′=0w'' = w''' = 0.
💡The shooting method
Turn the BVP into an IVP by guessing the missing initial conditions. For a fin, T′′=m2(T−T∞)T'' = m^2(T - T_\infty) with T(0)=TbT(0) = T_b and T(L)=TtipT(L) = T_{tip}: guess the base slope ss, integrate to x=Lx = L with an IVP solver (Unit 9), and measure the miss F(s)=T(L;s)−TtipF(s) = T(L; s) - T_{tip}. Then solve F(s)=0F(s) = 0 with a root finder (Unit 2). For the fin with m=20m = 20/m, slopes of −1000-1000 and −5000-5000 °C/m miss by +442+442 and −283-283 °C; brentq finds s=−3437.597127s = -3437.597127 °C/m, matching the analytic cosh⁡\cosh/sinh⁡\sinh solution to all digits shown. Shooting is natural and reuses good IVP solvers. Its weakness: if the ODE has growing solutions, small changes in ss produce enormous changes at the far end and the root finding becomes ill-conditioned — at large mm the shots diverge violently. Multiple shooting splits the interval and matches at interior points.
Derivation
Why the fin's miss is a straight line in ss — and why shooting can blow up
The fin equation is linear, so the solution started with slope ss is a fixed combination of two solutions: one with T(0)=TbT(0) = T_b, T′(0)=0T'(0) = 0, and one with T(0)=0T(0) = 0, T′(0)=1T'(0) = 1, scaled by ss. In the exact solution the second is sinh⁡(mx)/m\sinh(mx)/m. So
F(s)=T(L;s)−Ttip=F(0)+s sinh⁡mLm:F(s) = T(L; s) - T_{tip} = F(0) + s\,\frac{\sinh mL}{m}:
a straight line in ss. The secant method finds its root in a single step from any two shots — which is why two misses are enough in the lab, and brentq agrees to all digits.
The slope of that line is also the warning. A change δs\delta s in the starting slope changes the far end by δs sinh⁡(mL)/m\delta s\,\sinh(mL)/m, which grows like emL/(2m)e^{mL}/(2m): every extra unit of mLmL multiplies the sensitivity by e≈2.7e \approx 2.7, so going from mL=20mL = 20 to mL=40mL = 40 multiplies it by about e20≈5×108e^{20} \approx 5 \times 10^8. The root-finding problem becomes ill-conditioned (lesson 1-4) — not because of the root finder, but because the IVP amplifies everything. Finite differences impose both ends at once and never amplify anything.
💡Finite differences for BVPs
Replace derivatives by difference quotients (lesson 8-1) at the interior grid points. For the fin: Ti−1−2Ti+Ti+1h2−m2(Ti−T∞)=0,i=1,…,n−1,\dfrac{T_{i-1} - 2T_i + T_{i+1}}{h^2} - m^2(T_i - T_\infty) = 0, \qquad i = 1, \dots, n-1, with T0=TbT_0 = T_b and Tn=TtipT_n = T_{tip} moved to the right-hand side. That is a tridiagonal linear system — solved in O(n)O(n) by the Thomas algorithm (lesson 3-3). All the boundary conditions are imposed at once; nothing is guessed. The error is O(h2)O(h^2): at the fin's midpoint, n=10,20,…,160n = 10, 20, \dots, 160 intervals give errors that fall by 3.99, 4.00, 4.00, 4.00 per halving. For nonlinear BVPs the algebraic system is nonlinear and is solved by Newton's method — whose Jacobian is again tridiagonal or banded.
📐A propped cantilever with solve_bvp
A 6 m beam, clamped at x=0x = 0 and pinned at x=Lx = L, EI=40EI = 40 MN m², uniform load 20 kN/m. Write EIw′′′′=qEI w'''' = q as four first-order equations for (w,w′,w′′,w′′′)(w, w', w'', w''') and give scipy.integrate.solve_bvp the four conditions w(0)=w′(0)=w(L)=w′′(L)=0w(0) = w'(0) = w(L) = w''(L) = 0. It returns a maximum deflection of 3.5096 mm at x=3.471x = 3.471 m — exactly the textbook 0.005416 qL4/EI0.005416\,qL^4/EI at 0.5785L0.5785L from the clamp — and a fixed-end moment of 90.000 kN m =qL2/8= qL^2/8. solve_bvp uses collocation: a piecewise polynomial that satisfies the ODE at chosen points in each interval, with Newton's method on the resulting system and automatic mesh refinement.
Calculus systematically evades a great deal of numerical calculation.
– Richard Hamming, Methods of Mathematics Applied to Calculus, Probability, and Statistics (1985)
💡Where you will use this
Structures: beam deflections, columns on elastic foundations, plate strips — and the stiffness method is the finite-difference idea in matrix form. Heat transfer: fins, walls with heat generation, insulation design. Geotechnics: piles under lateral load (beams on Winkler springs) and one-dimensional consolidation profiles. Optimal Guidance: optimal control problems lead to two-point BVPs for the state and costate, solved by shooting or collocation.
📚Further reading
  • Uri Ascher, Robert Mattheij and Robert Russell, Numerical Solution of Boundary Value Problems for Ordinary Differential Equations (Prentice-Hall, 1988; SIAM Classics 1995) — shooting, multiple shooting, finite differences and collocation.
  • Randall LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations (SIAM, 2007), chapter 2, 'Steady States and Boundary Value Problems'.
  • Jacek Kierzenka and Lawrence Shampine, 'A BVP solver based on residual control and the MATLAB PSE', ACM Transactions on Mathematical Software 27(3), 299–316 (2001) — the algorithm inside scipy.integrate.solve_bvp.
↓But there is a problem
A beam or a fin has one spatial coordinate and no time. Most engineering fields have more: temperature through a heat shield that changes from second to second, pore pressure in a clay layer squeezed under a new building over decades, the concentration of a pollutant spreading through an aquifer. In 1947 John Crank and Phyllis Nicolson at Manchester published a careful method for exactly this kind of problem — heat conduction that changes in time — after finding that the obvious explicit approach forced absurdly small time steps. Their scheme became the textbook method for diffusion. It also turned out to hide a subtle flaw.