A boundary value problem (BVP) is a differential equation with conditions specified at 2 distinct points, called boundary conditions.
y ′ ′ = f ( x , y , y ′ ) , a ≤ x ≤ b , y ( a ) = α , y ( b ) = β y'' = f(x, y, y'), \quad a \leq x \leq b, \quad y(a) = \alpha,\ y(b) = \beta y ′′ = f ( x , y , y ′ ) , a ≤ x ≤ b , y ( a ) = α , y ( b ) = β Linear and Nonlinear BVPs
The BVP is linear if f f f
f ( x , y , y ′ ) = p ( x ) y ′ + q ( x ) y + r ( x ) f(x,y,y') = p(x)y' + q(x)y + r(x) f ( x , y , y ′ ) = p ( x ) y ′ + q ( x ) y + r ( x ) for some functions p p p q q q r r r
y ′ ′ = p ( x ) y ′ + q ( x ) y + r ( x ) , a ≤ x ≤ b , y ( a ) = α , y ( b ) = β y'' = p(x)y' + q(x)y + r(x), \quad a \leq x \leq b, \quad y(a) = \alpha,\ y(b) = \beta y ′′ = p ( x ) y ′ + q ( x ) y + r ( x ) , a ≤ x ≤ b , y ( a ) = α , y ( b ) = β Otherwise, the BVP is nonlinear.
Existence and Uniqueness
The BVP has a unique solution on [ a , b ] [a,b] [ a , b ] iff
f f f f y f_y f y f y ′ f_{y'} f y ′ D = { ( x , y , y ′ ) ∣ a ≤ x ≤ b , y , y ′ ∈ R } D = \{(x,y,y') \mid a \leq x \leq b,\ y, y' \in \mathbb{R}\} D = {( x , y , y ′ ) ∣ a ≤ x ≤ b , y , y ′ ∈ R } f y ( x , y , y ′ ) > 0 f_y(x,y,y') > 0 f y ( x , y , y ′ ) > 0 ( x , y , y ′ ) ∈ D (x,y,y') \in D ( x , y , y ′ ) ∈ D a constant M > 0 M > 0 M > 0 ∣ f y ′ ( x , y , y ′ ) ∣ ≤ M |f_{y'}(x,y,y')| \leq M ∣ f y ′ ( x , y , y ′ ) ∣ ≤ M ( x , y , y ′ ) ∈ D (x,y,y') \in D ( x , y , y ′ ) ∈ D
Examples:
y ′ ′ + e − x y + sin ( y ′ ) = 0 y'' + e^{-xy} + \sin(y') = 0 y ′′ + e − x y + sin ( y ′ ) = 0 0 ≤ x ≤ 1 0 \leq x \leq 1 0 ≤ x ≤ 1 y ( 0 ) = y ( 1 ) = 0 y(0) = y(1) = 0 y ( 0 ) = y ( 1 ) = 0 f y = x e − x y > 0 f_y = xe^{-xy} > 0 f y = x e − x y > 0 x , y x,y x , y ∣ f y ′ ∣ = ∣ cos ( y ′ ) ∣ ≤ 1 |f_{y'}| = |\cos(y')| \leq 1 ∣ f y ′ ∣ = ∣ cos ( y ′ ) ∣ ≤ 1 y ′ y' y ′
A linear BVP is solved by the linear shooting method or a linear finite-difference method . A nonlinear BVP is solved by the nonlinear shooting method or a nonlinear finite-difference method .
Linear Shooting Method
The linear shooting method solves a linear BVP by combining solutions to 2 initial value problems.
Applying the general existence and uniqueness theorem to f ( x , y , y ′ ) = p ( x ) y ′ + q ( x ) y + r ( x ) f(x,y,y') = p(x)y' + q(x)y + r(x) f ( x , y , y ′ ) = p ( x ) y ′ + q ( x ) y + r ( x ) f y = q ( x ) f_y = q(x) f y = q ( x ) f y ′ = p ( x ) f_{y'} = p(x) f y ′ = p ( x ) iff
p p p q q q r r r [ a , b ] [a,b] [ a , b ] q ( x ) > 0 q(x) > 0 q ( x ) > 0 x ∈ [ a , b ] x \in [a,b] x ∈ [ a , b ]
Superposition Idea
The ODE is linear, so if y 1 y_1 y 1
y 1 ′ ′ = p ( x ) y 1 ′ + q ( x ) y 1 + r ( x ) y_1'' = p(x)y_1' + q(x)y_1 + r(x) y 1 ′′ = p ( x ) y 1 ′ + q ( x ) y 1 + r ( x ) and y 2 y_2 y 2
y 2 ′ ′ = p ( x ) y 2 ′ + q ( x ) y 2 y_2'' = p(x)y_2' + q(x)y_2 y 2 ′′ = p ( x ) y 2 ′ + q ( x ) y 2 then y = y 1 + c y 2 y = y_1 + c\,y_2 y = y 1 + c y 2 c c c
Auxiliary IVPs
Fix
y 1 ( a ) = α , y 1 ′ ( a ) = 0 , y 2 ( a ) = 0 , y 2 ′ ( a ) = 1 y_1(a) = \alpha,\ y_1'(a) = 0, \qquad y_2(a) = 0,\ y_2'(a) = 1 y 1 ( a ) = α , y 1 ′ ( a ) = 0 , y 2 ( a ) = 0 , y 2 ′ ( a ) = 1 so that y ( a ) = y 1 ( a ) + c y 2 ( a ) = α y(a) = y_1(a) + c\,y_2(a) = \alpha y ( a ) = y 1 ( a ) + c y 2 ( a ) = α c c c
y 1 ′ ′ = p ( x ) y 1 ′ + q ( x ) y 1 + r ( x ) , a ≤ x ≤ b , y 1 ( a ) = α , y 1 ′ ( a ) = 0 y_1'' = p(x)y_1' + q(x)y_1 + r(x), \quad a \leq x \leq b, \quad y_1(a) = \alpha,\ y_1'(a) = 0 y 1 ′′ = p ( x ) y 1 ′ + q ( x ) y 1 + r ( x ) , a ≤ x ≤ b , y 1 ( a ) = α , y 1 ′ ( a ) = 0 y 2 ′ ′ = p ( x ) y 2 ′ + q ( x ) y 2 , a ≤ x ≤ b , y 2 ( a ) = 0 , y 2 ′ ( a ) = 1 y_2'' = p(x)y_2' + q(x)y_2, \quad a \leq x \leq b, \quad y_2(a) = 0,\ y_2'(a) = 1 y 2 ′′ = p ( x ) y 2 ′ + q ( x ) y 2 , a ≤ x ≤ b , y 2 ( a ) = 0 , y 2 ′ ( a ) = 1 by reducing each to a first-order system and applying order 4 Runge-Kutta . q ( x ) > 0 q(x) > 0 q ( x ) > 0 y 2 ( b ) ≠ 0 y_2(b) \neq 0 y 2 ( b ) = 0 y 2 y_2 y 2
Combining Solutions
c c c y ( b ) = β y(b) = \beta y ( b ) = β
y 1 ( b ) + c y 2 ( b ) = β ⟹ c = β − y 1 ( b ) y 2 ( b ) y_1(b) + c\,y_2(b) = \beta \implies c = \frac{\beta - y_1(b)}{y_2(b)} y 1 ( b ) + c y 2 ( b ) = β ⟹ c = y 2 ( b ) β − y 1 ( b ) y ( x ) = y 1 ( x ) + β − y 1 ( b ) y 2 ( b ) y 2 ( x ) y(x) = y_1(x) + \frac{\beta - y_1(b)}{y_2(b)}\, y_2(x) y ( x ) = y 1 ( x ) + y 2 ( b ) β − y 1 ( b ) y 2 ( x ) is the unique solution to the BVP.
Nonlinear Shooting Method
The nonlinear shooting method solves a nonlinear BVP by approximating the initial slope y ′ ( a ) y'(a) y ′ ( a )
Idea
For a guess t t t y ( x , t ) y(x,t) y ( x , t )
y ′ ′ = f ( x , y , y ′ ) , a ≤ x ≤ b , y ( a ) = α , y ′ ( a ) = t y'' = f(x,y,y'), \quad a \leq x \leq b, \quad y(a) = \alpha,\ y'(a) = t y ′′ = f ( x , y , y ′ ) , a ≤ x ≤ b , y ( a ) = α , y ′ ( a ) = t The method finds t t t y ( b , t ) − β = 0 y(b,t) - \beta = 0 y ( b , t ) − β = 0 Newton’s method .
Newton Iteration
t k = t k − 1 − y ( b , t k − 1 ) − β y ′ ( b , t k − 1 ) t_k = t_{k-1} - \frac{y(b,t_{k-1}) - \beta}{y'(b,t_{k-1})} t k = t k − 1 − y ′ ( b , t k − 1 ) y ( b , t k − 1 ) − β y ′ ( b , t k − 1 ) y'(b,t_{k-1}) y ′ ( b , t k − 1 ) y ( b , t k − 1 ) y(b,t_{k-1}) y ( b , t k − 1 ) z ( x , t ) = ∂ y ( x , t ) ∂ t z(x,t) = \dfrac{\partial y(x,t)}{\partial t} z ( x , t ) = ∂ t ∂ y ( x , t ) t t t
z ′ ′ = ∂ f ∂ y ( x , y , y ′ ) z + ∂ f ∂ y ′ ( x , y , y ′ ) z ′ , a ≤ x ≤ b , z ( a ) = 1 , z ′ ( a ) = 0 z'' = \frac{\partial f}{\partial y}(x,y,y')\,z + \frac{\partial f}{\partial y'}(x,y,y')\,z', \quad a \leq x \leq b, \quad z(a) = 1,\ z'(a) = 0 z ′′ = ∂ y ∂ f ( x , y , y ′ ) z + ∂ y ′ ∂ f ( x , y , y ′ ) z ′ , a ≤ x ≤ b , z ( a ) = 1 , z ′ ( a ) = 0 Substituting z ( b , t k − 1 ) z(b,t_{k-1}) z ( b , t k − 1 ) y ′ ( b , t k − 1 ) y'(b,t_{k-1}) y ′ ( b , t k − 1 )
t k = t k − 1 − y ( b , t k − 1 ) − β z ( b , t k − 1 ) t_k = t_{k-1} - \frac{y(b,t_{k-1}) - \beta}{z(b,t_{k-1})} t k = t k − 1 − z ( b , t k − 1 ) y ( b , t k − 1 ) − β Procedure
Choose an initial guess t 0 t_0 t 0
Repeat for k = 0 , 1 , 2 , … k = 0, 1, 2, \ldots k = 0 , 1 , 2 , … ∣ y ( b , t k ) − β ∣ |y(b,t_k) - \beta| ∣ y ( b , t k ) − β ∣
Approximate y k ( x ) = y ( x , t k ) y_k(x) = y(x,t_k) y k ( x ) = y ( x , t k )
y ′ ′ = f ( x , y , y ′ ) , a ≤ x ≤ b , y ( a ) = α , y ′ ( a ) = t k y'' = f(x,y,y'), \quad a \leq x \leq b, \quad y(a) = \alpha,\ y'(a) = t_k y ′′ = f ( x , y , y ′ ) , a ≤ x ≤ b , y ( a ) = α , y ′ ( a ) = t k
Approximate z k ( x ) = z ( x , t k ) z_k(x) = z(x,t_k) z k ( x ) = z ( x , t k )
z ′ ′ = ∂ f ∂ y ( x , y k , y k ′ ) z + ∂ f ∂ y ′ ( x , y k , y k ′ ) z ′ , a ≤ x ≤ b , z ( a ) = 1 , z ′ ( a ) = 0 z'' = \frac{\partial f}{\partial y}(x,y_k,y_k')\,z + \frac{\partial f}{\partial y'}(x,y_k,y_k')\,z', \quad a \leq x \leq b, \quad z(a) = 1,\ z'(a) = 0 z ′′ = ∂ y ∂ f ( x , y k , y k ′ ) z + ∂ y ′ ∂ f ( x , y k , y k ′ ) z ′ , a ≤ x ≤ b , z ( a ) = 1 , z ′ ( a ) = 0
Update
t k + 1 = t k − y k ( b ) − β z k ( b ) t_{k+1} = t_k - \frac{y_k(b) - \beta}{z_k(b)} t k + 1 = t k − z k ( b ) y k ( b ) − β
Linear Finite-Difference Method
The finite-difference method approximates the solution of a linear BVP by replacing derivatives with centered difference formulas at mesh points and solving the resulting linear system.
Mesh Points
Partition [ a , b ] [a,b] [ a , b ] N + 1 N+1 N + 1 h = b − a N + 1 h = \dfrac{b-a}{N+1} h = N + 1 b − a
x i = a + i h , i = 0 , 1 , … , N + 1 x_i = a + ih, \quad i = 0, 1, \ldots, N+1 x i = a + ih , i = 0 , 1 , … , N + 1 y ′ ′ ( x i ) = 1 h 2 ( y ( x i + 1 ) − 2 y ( x i ) + y ( x i − 1 ) ) − h 2 12 y ( 4 ) ( ξ i ) y''(x_i) = \frac{1}{h^2}\big(y(x_{i+1}) - 2y(x_i) + y(x_{i-1})\big) - \frac{h^2}{12}y^{(4)}(\xi_i) y ′′ ( x i ) = h 2 1 ( y ( x i + 1 ) − 2 y ( x i ) + y ( x i − 1 ) ) − 12 h 2 y ( 4 ) ( ξ i ) y ′ ( x i ) = 1 2 h ( y ( x i + 1 ) − y ( x i − 1 ) ) − h 2 6 y ′ ′ ′ ( η i ) y'(x_i) = \frac{1}{2h}\big(y(x_{i+1}) - y(x_{i-1})\big) - \frac{h^2}{6}y'''(\eta_i) y ′ ( x i ) = 2 h 1 ( y ( x i + 1 ) − y ( x i − 1 ) ) − 6 h 2 y ′′′ ( η i ) for some ξ i , η i ∈ ( x i − 1 , x i + 1 ) \xi_i, \eta_i \in (x_{i-1}, x_{i+1}) ξ i , η i ∈ ( x i − 1 , x i + 1 )
Linear Difference Equation
Substituting the centered difference formulas into the linear BVP and dropping the O ( h 2 ) O(h^2) O ( h 2 ) y ( x i ) ≈ w i y(x_i) \approx w_i y ( x i ) ≈ w i
w 0 = α , w N + 1 = β w_0 = \alpha, \quad w_{N+1} = \beta w 0 = α , w N + 1 = β − ( 1 + h 2 p ( x i ) ) w i − 1 + ( 2 + h 2 q ( x i ) ) w i − ( 1 − h 2 p ( x i ) ) w i + 1 = − h 2 r ( x i ) -\left(1 + \frac{h}{2}p(x_i)\right)w_{i-1} + \big(2 + h^2 q(x_i)\big)w_i - \left(1 - \frac{h}{2}p(x_i)\right)w_{i+1} = -h^2 r(x_i) − ( 1 + 2 h p ( x i ) ) w i − 1 + ( 2 + h 2 q ( x i ) ) w i − ( 1 − 2 h p ( x i ) ) w i + 1 = − h 2 r ( x i ) for i = 1 , … , N i = 1, \ldots, N i = 1 , … , N O ( h 2 ) O(h^2) O ( h 2 )
A w = b Aw = b A w = b with A A A
a i = − 1 − h 2 p ( x i ) , d i = 2 + h 2 q ( x i ) , c i = − 1 + h 2 p ( x i ) a_i = -1 - \frac{h}{2}p(x_i), \quad d_i = 2 + h^2 q(x_i), \quad c_i = -1 + \frac{h}{2}p(x_i) a i = − 1 − 2 h p ( x i ) , d i = 2 + h 2 q ( x i ) , c i = − 1 + 2 h p ( x i ) A = ( d 1 c 1 a 2 d 2 c 2 ⋱ ⋱ ⋱ a N d N ) A = \begin{pmatrix}
d_1 & c_1 & & \\
a_2 & d_2 & c_2 & \\
& \ddots & \ddots & \ddots \\
& & a_N & d_N
\end{pmatrix} A = d 1 a 2 c 1 d 2 ⋱ c 2 ⋱ a N ⋱ d N w = ( w 1 w 2 ⋮ w N − 1 w N ) , b = ( − h 2 r ( x 1 ) − a 1 α − h 2 r ( x 2 ) ⋮ − h 2 r ( x N − 1 ) − h 2 r ( x N ) − c N β ) w = \begin{pmatrix} w_1 \\ w_2 \\ \vdots \\ w_{N-1} \\ w_N \end{pmatrix}, \quad
b = \begin{pmatrix}
-h^2r(x_1) - a_1\alpha \\
-h^2r(x_2) \\
\vdots \\
-h^2r(x_{N-1}) \\
-h^2r(x_N) - c_N\beta
\end{pmatrix} w = w 1 w 2 ⋮ w N − 1 w N , b = − h 2 r ( x 1 ) − a 1 α − h 2 r ( x 2 ) ⋮ − h 2 r ( x N − 1 ) − h 2 r ( x N ) − c N β Solve for w w w LU decomposition or another linear solver.
If p p p q q q r r r [ a , b ] [a,b] [ a , b ] q ( x ) ≥ 0 q(x) \geq 0 q ( x ) ≥ 0 [ a , b ] [a,b] [ a , b ] h < 2 L h < \dfrac{2}{L} h < L 2 L = max a ≤ x ≤ b ∣ p ( x ) ∣ L = \max_{a \leq x \leq b} |p(x)| L = max a ≤ x ≤ b ∣ p ( x ) ∣
Nonlinear Finite-Difference Method
The finite-difference method for a nonlinear BVP replaces derivatives with centered difference formulas and solves the resulting nonlinear system with Newton’s method.
A unique solution exists iff
f f f f y f_y f y f y ′ f_{y'} f y ′ x ∈ [ a , b ] x \in [a,b] x ∈ [ a , b ] y , y ′ ∈ ( − ∞ , ∞ ) y, y' \in (-\infty, \infty) y , y ′ ∈ ( − ∞ , ∞ ) f y ( x , y , y ′ ) ≥ δ f_y(x,y,y') \geq \delta f y ( x , y , y ′ ) ≥ δ δ > 0 \delta > 0 δ > 0 a constant k k k k = max ( ∣ f y ( x , y , y ′ ) ∣ , ∣ f y ′ ( x , y , y ′ ) ∣ ) k = \max\big(|f_y(x,y,y')|,\ |f_{y'}(x,y,y')|\big) k = max ( ∣ f y ( x , y , y ′ ) ∣ , ∣ f y ′ ( x , y , y ′ ) ∣ )
Nonlinear Difference Equation
Using the centered difference formulas at mesh points x 0 , … , x N + 1 x_0, \ldots, x_{N+1} x 0 , … , x N + 1 [ a , b ] [a,b] [ a , b ] y ( x i ) ≈ w i y(x_i) \approx w_i y ( x i ) ≈ w i
w 0 = α , w N + 1 = β w_0 = \alpha, \quad w_{N+1} = \beta w 0 = α , w N + 1 = β − w i + 1 − 2 w i + w i − 1 h 2 + f ( x i , w i , w i + 1 − w i − 1 2 h ) = 0 -\frac{w_{i+1} - 2w_i + w_{i-1}}{h^2} + f\left(x_i, w_i, \frac{w_{i+1}-w_{i-1}}{2h}\right) = 0 − h 2 w i + 1 − 2 w i + w i − 1 + f ( x i , w i , 2 h w i + 1 − w i − 1 ) = 0 for i = 1 , … , N i = 1, \ldots, N i = 1 , … , N O ( h 2 ) O(h^2) O ( h 2 )
Nonlinear System
2 w 1 − w 2 + h 2 f ( x 1 , w 1 , w 2 − α 2 h ) − α = 0 − w 1 + 2 w 2 − w 3 + h 2 f ( x 2 , w 2 , w 3 − w 1 2 h ) = 0 ⋮ − w N − 2 + 2 w N − 1 − w N + h 2 f ( x N − 1 , w N − 1 , w N − w N − 2 2 h ) = 0 − w N − 1 + 2 w N + h 2 f ( x N , w N , β − w N − 1 2 h ) − β = 0 \begin{aligned}
2w_1 - w_2 + h^2 f\left(x_1, w_1, \frac{w_2 - \alpha}{2h}\right) - \alpha &= 0 \\
-w_1 + 2w_2 - w_3 + h^2 f\left(x_2, w_2, \frac{w_3 - w_1}{2h}\right) &= 0 \\
&\ \ \vdots \\
-w_{N-2} + 2w_{N-1} - w_N + h^2 f\left(x_{N-1}, w_{N-1}, \frac{w_N - w_{N-2}}{2h}\right) &= 0 \\
-w_{N-1} + 2w_N + h^2 f\left(x_N, w_N, \frac{\beta - w_{N-1}}{2h}\right) - \beta &= 0
\end{aligned} 2 w 1 − w 2 + h 2 f ( x 1 , w 1 , 2 h w 2 − α ) − α − w 1 + 2 w 2 − w 3 + h 2 f ( x 2 , w 2 , 2 h w 3 − w 1 ) − w N − 2 + 2 w N − 1 − w N + h 2 f ( x N − 1 , w N − 1 , 2 h w N − w N − 2 ) − w N − 1 + 2 w N + h 2 f ( x N , w N , 2 h β − w N − 1 ) − β = 0 = 0 ⋮ = 0 = 0 This is a system of N N N N N N h < 2 L h < \dfrac{2}{L} h < L 2 Newton’s method .
