Chapter 3

Finite difference methods for two-point boundary value problems Let us consider a model problem u (x) = f (x),

0 < x < 1,

u(0) = ua ,

u(1) = ub ,

to illustrate the general procedure using a finite difference method as follows. Note that, if f (x) is continuous in [0, 1], then the solution to the ODE with the given boundary condition exists and it is unique. Furthermore, the solution u(x) has continuous up to second order derivatives in [0, 1]. There are many application of one-dimensional second order two-point boundary value problems. • Physical problems such as the pendulum model, the elastic string model, etc. • For high dimensional axi-symmetric problems, the elliptic second partial differential equations can be written as a two-point boundary value problems. • The steady state solution of a boundary value problem of the heat equation equation ∂u ∂2u = c2 2 , ∂t ∂x u(x, 0) = u0 (x),

a 0. If u denotes velocity, then r(x)u (x) may be called the advection term. When the advection is strong (i.e., |r(x)| is large), the approach to a wave equation can lead to solutions with non-physical oscillations, e.g., when ri ∼ 1/h. • Upwinding discretization for the first order derivative and central FD scheme for the diffusion term:

pi

Ui−1 − 2Ui + Ui+1 Ui+1 − Ui − qi Ui = fi , if ri ≥ 0, + ri 2 h h

pi

Ui−1 − 2Ui + Ui+1 Ui − Ui−1 − qi Ui = fi , if ri < 0. + ri h2 h

This scheme increases the diagonal dominance, but it is only first order accurate. If |r(x)| is very large (say |r(x)| ∼ 1/h), it is best to solve the linear system of equations with either a direct or an iterative method for diagonally dominant matrices.

3.6

The ghost point method for boundary conditions involving derivatives

In this section, we discuss how to treat Neumann and mixed (Robin) boundary conditions. Let us first consider the problem u (x) = f (x), u (a) = α,

a < x < b, u(b) = ub ,

where the solution at x = a is unknown. If we use a uniform Cartesian grid xi = a + ih, then U0 is one component of the solution. We can still use Ui−1 − 2Ui + Ui+1 = fi , h2

i = 1, 2, · · · n − 1,

but we need an additional equation at x0 = a given the Neumann boundary condition at a. One approach is to take U 1 − U0 =α h

or

−U0 + U1 α = , 2 h h

(3.1)

54

Z. Li

and the resulting linear system of equations is again tri-diagonal and symmetric negative definite: ⎤ ⎡ ⎤⎡ ⎤ ⎡ α 1 U − h12 0 2 h h ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ 1 ⎥ ⎥ 1 ⎢ ⎥ ⎢ ⎢ h2 − h22 ⎥ ⎥ U ) f (x 1 1 2 ⎢ ⎥ h ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ 1 2 1 ⎢ ⎢ ⎥ ⎥ ⎢ ⎥ − U ) f (x 2 ⎥ 2 ⎢ ⎥⎢ h2 h2 h2 ⎥ ⎢ ⎥⎢ ⎥=⎢ (3.2) ⎢ ⎥ ⎢ ⎥⎢ . ⎥ ⎢ . . . . ⎥ .. .. .. .. ⎢ ⎥ ⎢ .. ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ 1 2 1 ⎥ − h2 f (xn−2 ) ⎢ ⎥ ⎢ Un−2 ⎥ ⎢ 2 2 h h ⎥ ⎣ ⎦⎣ ⎦ ⎣ ⎦ u b 1 2 − U f (x ) − n−1 2 2 n−1 h h h2 However, this approach is only first order accurate if α = 0. To maintain second order accuracy, the ghost point method is recommended, where a ghost grid point x−1 = x0 − h = a − h is added and the solution is extended to the interval [a − h, a]. Then the central FD scheme can be used at all grid points (where the solution is unknown), i.e., for i = 0, 1, · · · n. However, we now have n equations and n + 1 unknowns including U−1 , so one more equation is needed to close the system. The additional equation is the central finite difference equation for the Neumann boundary condition U1 − U−1 = α, (3.3) 2h which yields U−1 = U1 − 2hα. Inserting this into the FD equation for the differential equation at x = a (i.e. at x0 ), we have U−1 − 2U0 + U1 = f0 , h2 U1 − 2hα − 2U0 + U1 = f0 , h2 f0 α −U0 + U1 + , = h2 2 h where the coefficient matrix is precisely the same as for (3.2) and the only difference in this second order method is the component f0 /2 + α/h in the vector on the right-hand side, rather than α/2 in the previous first order method. We compare the two methods in Figure 3.2, where the differential equation u (x) = f (x) is subject to a Dirichlet boundary condition at x = 0 and a Neumann boundary condition at x = 0.5 (the exact solution is of course u(x) = cos πx). Figure 3.2 (a) shows the grid refinement analysis. The error in the second order method is evidently much smaller than in the first order method. In Figure 3.2 (b), we show the error plot of the ghost point method, and note the error at x = b is no longer zero.

The ghost point method

55

(a)

(b) 0

−5

The error plot in log−log scale for 1−st and 2−nd order methods

10

0

−1

10

Error

x 10

−0.5

The 1−st order method, slope = 1

−2

10

−1

−3

10

−1.5

−4

10

−2 The 2−nd order method, slope = 2

−4

10

−3

10

−2

10

−1

10

0

10

−2.5

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

Figure 3.2: (a) A grid refinement analysis of the ghost point method and the first order method. The slopes of the curves are the order of convergence. (b) The error plot of the computed solution from the ghost point method.

3.6.1

A Matlab code of the ghost point method

function [x,U] = two_pointa(a,b,ua,uxb,f,n) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This matlab function two_point solves the following two-point % % boundary value problem: u’’(x) = f(x) using the center finite % % difference scheme. % % Input: % % a, b: Two end points. % % ua, uxb: Dirichlet and Neumann boundary conditions at a and b % % f: external function f(x). % % n: number of grid points. % % Output: % % x: x(1),x(2),...x(n) are grid points % % U: U(1),U(2),...U(n) are approximate solution at grid points. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% h = (b-a)/n; h1=h*h; A = sparse(n,n); F = zeros(n,1); for i=1:n-1, A(i,i) = -2/h1; A(i+1,i) = 1/h1; A(i,i+1)= 1/h1;

56

Z. Li

end A(n,n) = -2/h1; A(n,n-1) = 2/h1; for i=1:n, x(i) = a+i*h; F(i) = feval(f,x(i)); end F(1) = F(1) - ua/h1; F(n) = F(n) - 2*uxb/h; U = A\F; return

3.6.2

The Matlab driver program

Below is a Matlab driver code to solve the two-point boundary value problem u (x) = f (x),

a < x < b,

u (a) = ua,

u(b) = uxb.

%%%%%%%% Clear all unwanted variable and graphs. clear;

close all

%%%%%%% Input a =0; b=1/2; ua = 1; uxb = -pi; %%%%%% Call solver: U is the FD solution n=10; k=1; for k=1:5 [x,U] = two_point(a,b,ua,uxb,’f’,n); u=zeros(n,1);

%ghost-point method.

An example of a non-linear BVP

57

for i=1:n, u(i) = cos(pi*x(i)); end h(k) = 1/n; e(k) = norm(U-u,inf); %%% Print out the maximum error. k = k+1; n=2*n; end loglog(h,e,h,e,’o’); axis(’equal’); axis(’square’), title(’The error plot in log-log scale, the slope = 2’); % Ghost point figure(2); plot(x,U-u); title(’Error’)

3.6.3

Dealing with mixed boundary conditions

The ghost point method can be used to discretize a mixed boundary condition. Suppose that αu (a) + βu(b) = γ at x = a, where α = 0. Then we can discretize the boundary condition by α or

U1 − U−1 + βU0 = γ, 2h

U−1 = U1 +

2βh 2hγ U0 − α α

and substitute this into the central FD equation at x = x0 :  2β 2 2γ 2 , − 2+ U0 + 2 U 1 = f 0 + h αh h αh  1 γ β 1 f0 or − 2+ U0 + 2 U1 = + , h αh h 2 αh

(3.4) (3.5)

yielding a symmetric coefficient matrix.

3.7

An example of a non-linear boundary value problem

Discretrizing a non-linear differential equation generally produces a non-linear algebraic system. Furthermore, if we can solve this system then we can get an approximate solution. We now present an example in this section to illustrate the procedure. Consider the simple pendulum model θ¨ + K sin θ = 0

58

Z. Li

˙ we are given the angle at the with θ(0) = θ0 . But instead of given the initial velocity θ(0), final time T , θ(T ) = θ1 , then we have a two-point boundary value problem. If we apply the central FD scheme, then we obtain the system of non-linear equations θi−1 − 2θi + θi+1 + K sin θi = 0, h2

i = 1, 2, · · · , n − 1,

(3.1)

where h = T /n is the step size. The non-linearity comes from sin θi . The non-linear system of equations above can be solved in several ways: • Approximate linearized model. For example, we can approximate θ + K sin θ = 0 using the linear equation θ + Kθ = 0, justifying this by noting that sin θ = θ + O(θ3 ) for θ small. Unfortunately, this liberalization loses a key feature of the non-linear solution, which depends on the pendulum oscillation amplitude in addition to K. • Substitution method, where the non-linear term depends upon the previous approximation. Thus given an initial guess θ(0) (x), form the iteration (θk+1 ) + K sin θk = 0,

k = 0, 1, · · · ,

(3.2)

involving a two-point boundary value problem at each iteration. The main concerns are then whether or not the method converges, and the rate of the convergence if it does. • Solve the non-linear problem directly, as now explained below. In general, a non-linear system of equations F(U) = 0 is obtained if we discretize a non-linear ODE or PDE, i.e., ⎧ ⎪ F1 (U1 , U2 , · · · , Un ) = 0, ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ F2 (U1 , U2 , · · · , Un ) = 0, (3.3) .. .. .. .. ⎪ ⎪ . . . . ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ F (U , U , · · · , U ) = 0, n

1

2

n

which generally is solved by Newton’s method. Given an initial guess U(0) , the Newton iteration is U(k+1) = U(k) − (J(U(k) ))−1 F(U(k) ) or

⎧ ⎨ J(U(k) )ΔU(k) = −F(U(k) ), ⎩ U(k+1) = U(k) + ΔU(k) ,

k = 0, 1, · · ·

(3.4)

An example of a non-linear BVP

59

where J(U(k) ) is the Jacobi matrix defined ⎡ ∂F1 ∂F1 ⎢ ∂U1 ∂U2 ⎢ ⎢ ⎢ ∂F2 ∂F2 ⎢ ⎢ ∂U1 ∂U2 ⎢ ⎢ . .. ⎢ .. . ⎢ ⎢ ⎣ ∂Fn ∂Fn ∂U1 ∂U2

as ··· ··· .. . ···

⎤

∂F1 ∂Un ∂F2 ∂Un .. .

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥. ⎥ ⎥ ⎥ ⎥ ⎦

∂Fn ∂Un

For the pendulum problem, we have ⎡ −2 + h2 K cos θ1 1 ⎢ ⎢ 1 −2 + h2 K cos θ2 1 ⎢ ⎢ J(θ) = 2 ⎢ .. h ⎢ . ⎢ ⎣

⎤ 1 ..

..

.

.

−2 + h2 K cos θn−1

1

⎥ ⎥ ⎥ ⎥ ⎥. ⎥ ⎥ ⎦

In Figure 3.3, we show the numerical solution of the pendulum problem using Newton’s method with θ(0) = 0.7, θ(2π) = 0.5. Several approximate solutions of the Newton method are shown in the left plot. In the right plot, approximate solutions using the linear and non-linear differential equations are shown. (a)

(b) nu = 10

Final computed solutions of linear and non−linear ODE

1

0.8

0.8

0.6 initial

0.6 0.4 0.4 maximum error = 3.6459e−02

0.2

0.2

0

0

−0.2

−0.2

−0.4 −0.4 −0.6 −0.6

−0.8

−1

0

1

2

3

4

5

6

7

−0.8

0

1

2

3

4

5

6

7

Figure 3.3: (a) Plot of some intermediate solution to the non-linear system of equations. (b) Comparison of the solution of the linear and non-linear solution to the pendulum motion. It is not always easy to find J(U) and it can be computationally expensive. Furthermore, Newton’s method is only locally convergent, i.e., it requires a close initial guess to guarantee its convergence. Quasi-Newton methods, such as the Broyden and BFGS rank-one and rank-two update methods, and also conjugate gradient methods, can avoid evaluating the Jacobian matrix, The main issues remain global convergence, the convergence order

60

Z. Li

(Newton’s method is quadratically convergent locally), and computational issues (storage, etc.). A well known software package called MINPACK is available through the netlib, cf., [?] for more complete discussions.

Index autonomous, 31 backward Euler’s method, 24 consistency, 23 convergency, 24 Crank-Nicholson method, 25 equilibrium, 31 Euler’s method, 21 explicit finite difference method, 21 initial condition, 15 local truncation error, 23 Lotka-Volterra model, 18 Matlab ODE-Suite, 29 Newton cooling model, 15 pendulum model, 17 Runge-Kutta Method, 26 stability condition, 23 steady state solution, 31

61