Class 17.

ODEs, Part 4 (2-pt BVPs)

Two-point Boundary Value Problems
• NRiC §17.

• BCs specified at two or more points, e.g., start and end.

• For IVP, just integrate away.

• For 2-pt BVP, must make a free choice of unknown BVs at initial point, then integrate
away. But solution will almost certainly not satisfy other BCs at end.

• Strategy: Use information about how much the other BVs “missed” to iteratively
improve initial guess.

=⇒ Techniques are all iterative (and expensive).

Notation
• Denote standard system as:

dyi (x)
= gi (x, y1 , ..., yN ) i = 1, ..., N.
dx

• At x1 , the solution is supposed to satisfy:

B1j (x1 , y1 , ..., yN ) = 0 j = 1, ..., n1 .

• At x2 , it is supposed to satisfy:

B2k (x2 , y1 , ..., yN ) = 0 k = 1, ..., n2 ,

where n2 = N − n1 .

Two Basic Techniques

Shooting method
1. Begin at x1 .

2. Guess values for free BCs (n2 values).

3. Integrate as IVP to x2 .

4. Adjust n2 guesses to get closer to BVs at x2 .

1
 1
y
3 00111100Desired BV
2
0011
Required BV x

• Heart of technique: system of iteratively improving guesses.

=⇒ Multi-D root finding.

Relaxation method
1. Replace ODEs by finite-difference equations on mesh from x1 to x2 .

2. Guess solution on this mesh.

3. Mathematically, FDEs are just algebraic relations between unknowns. Use iterative
technique to relax this solution to true solution.

2
y 3
true solution 00111100Required BV

00111100
Required BV x

• Relaxation very powerful for smooth solutions, or ODEs that must be solved many
times with different parameter values. Also good when ODEs have extraneous solu-
tions, i.e., stiff equations.

• NRiC : “Shoot first, relax later.”

2
2-pt BVP: Shooting Method
Procedure (NRiC §17.1):

1. At x1 , must specify N starting values for yi , i = 1, ..., N .

• n1 values given by BC at x1 .
• ∴ n2 = N − n1 values can be chosen freely.

2. Represent the free values as a vector V of dimension n2 (actually, V represents schemat-

ically any parameter values that specify unknown BVs).

3. Now integrate to x2 .

4. Define “discrepancy vector” F of dimension n2 , where

Fk = B2k (x2 , y) k = 1, ..., n2 .

• We want to find V that zeroes F.

5. Solve n2 linear equations:

J δV = −F,
where Jij ≡ ∂Fi /∂Vj is the Jacobian matrix.

• This is the globally convergent Newton’s method (NRiC §9.7).

6. Then Vnew = Vold + δV.

7. Use Vnew to solve ODEs again as IVP, recompute F, and iterate again until |F| < ε,
the convergence criterion.

• Infeasible to compute Jacobian analytically. Instead evaluate differences numerically,

i.e.,
∂Fi Fi (V1 , ..., Vj + ∆Vj , ..., Vn2 ) − Fi (V1 , ..., Vj , ..., Vn2 )
≃ ,
∂Vj ∆Vj
i.e., solve IVP n 2 times, varying each component of V by ±∆V each time to build
up Jacobian (recall Fi (V1 , ..., Vn2 ) already computed in step 4).

• Overall procedure requires n2 + 1 solutions to ODEs per iteration.

• For linear systems, one iteration is enough.

• For nonlinear systems, many (say M ) iterations may be required to converge =⇒

M × (n2 + 1) solutions of ODEs!

• ∴ need efficient integrator... (NRiC routine shoot() uses odeint()).

3
 • NOTE: Can also shoot to a fitting point xf between x1 and x2 (NRiC §17.2; specify
known points at x1 and x2 , choose the rest, integrate in both directions, and require
that y(xf ) match for both integrations). Useful for singular BC(s) and/or domain
point(s); integrate away from these.

y 00111100

11
00
x1 xf x2

2-pt BVP: Relaxation Method

Procedure (NRiC §17.3):

1. Replace ODEs
dy(x)
= g(x, y)
dx
with finite difference equations (FDEs) on a grid:1
 
yk − yk−1 xk + xk−1 yk + yk−1
=g , ,
xk − xk−1 2 2
where yk refers to the entire set of dependent variables y1 , y2 , . . . , yN at point xk .
• Here xk and the components of yk are discrete values of independent and depen-
dent variables at “mesh points.”

1 2 3 ... M

• For M mesh points and N coupled equations, have M × N yk components to

solve for. Approximate the set of N first-order ODEs by
0 = Ek ≡ yk − yk−1 − (xk − xk−1 )gk (xk , xk−1 , yk , yk−1 ), k = 2, 3, . . . , M.
Here gk can be evaluated using information from both points k, k − 1.
1
This is not a unique representation. We could, for example, evaluate g at (xk , yk ) and (xk−1 , yk−1 ),
then take the average.

4
 • This is (M − 1)N equations; get remaining equations from the boundary condi-
tions:
0 = E1 ≡ B(x1 , y1 )
0 = EM +1 ≡ C(xM , yM )

2. The “solution” of the FDE problem consists of a set of variables yj,k . Need to guess
starting values for all yj,k . We then determine increments ∆yj,k such that yj,k + ∆yj,k
is an improved approximation to the solution.
3. To get increments, expand FDEs in first-order Taylor series about yk (N-R method):
N N
X ∂Ek X ∂Ek
Ek (yk + ∆yk , yk−1 + ∆yk−1 ) ≃ Ek (yk , yk−1 ) + ∆yn,k−1 + ∆yn,k .
n=1
∂y n,k−1 n=1
∂y n,k

Similar relations can be obtained for the boundary conditions (see NRiC §17.3 for
details; the partial derivatives can be computed analytically—it’s just tedious! NRiC
§17.4 gives a worked example).
4. Want E(y + ∆y) = 0. Result is a large (M × N ) × (M × N ) block-diagonal matrix2
that can be solved using optimized Gaussian elimination for the ∆y’s.
5. After applying the new increments, iteratively improve (“relax”) solution until the
boundary conditions are satisfied and the difference equations between grid points are
zeroed to the desired accuracy.
• Need to solve a matrix equation each iteration.
• Choice of grid points is an important issue and leads to adaptive mesh strategies
in modern solvers.

Example: Stellar Structure

• Numerical methods for 2-pt BVPs largely developed by astronomers seeking to solve
equations of stellar structure.
– Form a system of four coupled ODEs.
1. Consider spherical shell, thickness dr, distance r from origin. Then dMr = 4πr2 drρ,
or,
dMr
= 4πr2 ρ.
dr
It is convenient to transform this equation (and the rest) so that Mr is the independent
variable. In this case, just take the reciprocal:
dr 1
= .
dMr 4πr2 ρ
2
Each interior point supplies a block of N equations coupling 2N corrections to the solution variables at
points k, k − 1. The boundary conditions supply smaller blocks, n1 × N and n2 × N .

5
2. Hydrostatic equilibrium =⇒ net force on shell is zero. ∴ −∇r P − ρg = 0, where g =
gravitational acceleration per unit mass = GMr /r2 , or,
dP GMr
= − 2 ρ.
dr r
(To derive, note upward force on shell per unit area = P (r) − P (r + ∆r) = −∆P
must equal downward force per unit area = [GMr /r2 ](Mshell /4πr2 ) = (GMr /r2 )ρ dr.)
Transforming,
dP GMr
=− .
dMr 4πr4
3. Let ε = energy generation rate/unit mass. Then energy transport rate through shell
∆L = L(r + ∆r) − L(r) must equal energy generation rate 4πr2 dr ρε, where L =
luminosity, or,
dLr
= 4πr2 ρε.
dr
Transforming,
dLr
= ε.
dMr
4. Finally, there is an equation describing energy transport. For radiative (and conduc-
tive) transport,
dT κρLr
∝ 2 3,
dr r T
where κ is the mean absorption coefficient (opacity; so higher opacity =⇒ higher T
gradient). Transforming,
dT κLr
∝ 4 3.
dMr r T
• This equation harder to derive since it depends on radiation transport mechanism
and convective stability.

• We have 4 ODEs in 7 unknowns (r, ρ, P, Lr , T, ε, κ).

• Need 3 constitutive relations:

1. P = P (ρ, T ) — equation of state.

2. ε = ε(ρ, T ) — nuclear energy generation rate.
3. κ = κ(ρ, T ) — opacity (for radiative transport; otherwise need an equivalent
relation for convection).

• Boundary conditions (need 4):

– At center (Mr = 0): r = 0, Lr = 0. But P = Pc = ?, T = Tc = ?

– At surface (Mr = M ): P ≃ 0, T ≃ 0 (also ρ ≃ 0). But r = R = ?, Lr = L = ?

• This is a classic 2-pt BVP.

6
 • First solving techniques based on shooting method.

• Singularity at center =⇒ must fit at an intermediate point: Schwarzschild Scheme

(e.g., Schwarzchild, Structure and Evolution of the Stars).

• Modern stellar structure codes use relaxation method with adaptive mesh (e.g., P.
Eggleton, MNRAS, 151, 351 (1971)).

Polytropes
• Can illustrate technique with calculation of structure of polytropes.

• Assume there is no energy generation anywhere inside (ε ≡ 0), e.g., white dwarf or
neutron star.

• Assume EOS of form P = kρ(n+1)/n (no T dependence).

– E.g., EOS for monatomic gas (such as degenerate electron gas) is:
P ∝ ρ5/3 if non-relativistic (n = 3/2);
P ∝ ρ4/3 if relativistic (n = 3).

• This form is called a polytropic EOS. n is polytropic index.

• It is convenient to rewrite ρ = θn . Then the first stellar structure equation becomes

dM
= R2 θ n ,
dR
and the second becomes
dθ M
= − 2,
dR R
3 2
where R ≡ r/r0 , M = 4πr0 M, and r0 = (n + 1)k/4πG (E.F.T.S.).

• These are the “Lane-Emden Equations.”

• This is a system of 2 ODEs with BC M = 0 at R = 0 and θ = 0 at R = R⋆ /r0 .

– If we have a desired R⋆ known in advance (or, equivalently, M⋆ ), then we can

set θ = θc (say) at R = 0 and iterate over different starting θc until the outer
boundary condition is satisfied.