Solution of Two-Point Boundary Value Problems: A Root-Solving Approach

Robert D. Tzou
Professor and Chair
Department of Mechanical and Aerospace Engineering
University of Missouri-Columbia
Columbia, MO 65211
http://web.missouri.edu/~mae/; TzouR@missouri.edu

INTRODUCTION

NDSolve is a powerful embedded function in Mathematica that can be used to solve both linear and nonlinear differential equations.  For problems involving linear differential equations, NDSolve is effective for both initial and boundary value problems.  For problems involving nonlinear differential equations, however, NDSolve can only be applied to solve initial value problems.  This notebook aims to extend the use of NDSolve to obtain the numerical solutions for the two-point boundary value problems along with the other embedded functions in Mathematica.  Without attributing to additional numerical algorithms such as the shooting method, this approach is self-contained within the context of Mathematica, which is valuable in engineering education because the students can better focus on the physics of the solution rather than the numerical methods.

	The sample problem in this notebook employs a nonlinear ordinary differential equation encountered in the boundary layer theory of fluid flow.  More intensive applications for nonlinear boundary value problems in forced and natural convective heat transfer can be found in the web page http://web.missouri.edu/~mae/.

ROOT-SOLVING APPROACH

The following third-order system is encountered in the boundary layer theory of fluid flow over a flat plate:

				f'''+(1/2)ff'' = 0, where f' = (df/(dη)),  f'' = (((d^2) f)/d)(η^2), etc.				(1)

The boundary conditions to be satisfied by f are

				f = 0; f' = 0 at η = 0, and f' = 1 at η = 10				(2)

For an explicit example, the boundary condition for f' at η = 10 replaces the remote boundary condition as η -> ∞.

	Since equation (1) is nonlinear, direct use of NDSolve necessitates recasting of equation (2) into  three initial conditions:

				f = 0; f' = 0; f'' = A at η = 0.						(3)

where the initial slope A is an unknown to be determined.  Obviously, the solution curve for f, and hence the end value of f' at η = 10, depends on the value of A.  The appropriate value of A must render f' = 1 at η = 10, in accordance with the original boundary condition stated in equation (2).  This behavior is further illustrated in the following figure.



Along with f(0) = 0 and f'(0) = 0 specified at η = 0 and f''(0) = A presumed at η = 0, NDSolve can determine the curves of f, f', and f'' in the entire domain of η from 0 to xend.  The curve of f' is focused in the figure because f' is involved in the boundary condition at η = xend; f'(xend) - 1 = 0.  The end value f'(xend), and hence fpend = f'(xend) - 1, varies with the presumed initial slope A.  Suppose a series of A values is selected in the use of NDSolve; in a window from Al to Au with an increment Ainc, an interpolation function can be established for fpend(A).  The task is then the determination of the root A that satisfies fpend(A) = 0 or f'(xend) - 1 = 0.  In Mathematica, the initial value problem is solved by NDSolve, the interpolation function for fpend(A) is accomplished by Interpolation, and the root finding procedure can be supported by findRoot.  To avoid extrapolation in determining the root, Al and Au (defining the window size of A) should be chosen so that fpend(A) changes sign across f'(xend) = 1.

IMPLEMENTATIONS IN Mathematica

1.  Define the size of the window for A - Lower bound of A = Al = 0; Upper bound  of A = Au = 1; Increment of A = Ainc, and xend = 10:

xend = 10; Al = 0; Au = 1; Ainc = (Au - Al)/50  // N;

2.  Define fpend as a function of A, which is f'[xend]  obtained by NDSolve:

fpend[A_] := f'[xend]  /.  
        NDSolve[{f'''[x] + 0.5 f[x] f''[x] == 0, 
            f[0] == 0, f'[0] == 0, f''[0] == A}, 
          f, {x, 0, xend}]

3.  Tabulate the results of fpend for A ∈ [Al, Au] with an increment Ainc; note that the cross-over (fpend = 1) does occur for A between 0.32 and 0.33:

Table[{A, fpend[A]}, {A, Al, Au, Ainc}]

4.  Interpolate the tabulated results and store in fpA:

fpA = Interpolation[%]

5.  Make the interpolated results a function of A:

fpA[A]

6.  The appropriate value of A satisfies fpA[A] = 1, with A in the same domain from Al to Au:

FindRoot[fpA[A] == 1, {A, {Al, Au}}]

7.  Check the result by using the value of A obtained above:

NDSolve[{f'''[x] + 0.5 f[x] f''[x] == 0, 
        f[0] == 0, f'[0] == 0, f''[0] == A  /. %}, 
      f, {x, 0, xend}]

8.  Plot the function f'[x] in the entire domain of x from 0 to xend:

Plot[Evaluate[f'[x] /. %], {x, 0, xend}, 
      AxesLabel -> {"eta", "f(eta)"}, 
      PlotRange -> {0, 1}]

CONCLUSIONS

(1) Comparing to the result of A = 0.332 given in most textbooks of Heat Transfer, the result of A = 0.332058 obtained by the root-solving approach is accurate.  The plot shown above is exactly the same as that obtained by the shooting method, where f'[x] is the dominating flow velocity in the direction parallel to the plate.  If intended, f[x] can be used in (8.0) to plot the flow velocity in the direction perpendicular to the plate.  The boundary condition f' = 1 is achieved for η > 6, implying that the use of xend = 10 is sufficient to model the remote distance (η -> ∞).  The use of xend < 5 will not give good result for the flow field with an infinite extent.

(2) The root-solving approach employs the embedded functions of NDSolve, Interpolation, and FindRoot in Mathematica, totally eliminating the need for decomposing the high-order ODE into numerous first-order ODE's and the use of the shooting method.  The root-solving approach facilitates a better focus on the physics of the problem rather than the numerical method and its related procedure.

(3) Defining the window for A, Al and Au in particular, is a task parallel to the initially assumed slope in the shooting method.  Specifying a bad range for A will lead to failure of the numerical procedure, which is also the case for all the other numerical methods (including the shooting method) should bad initial values be assumed.  In the present root-solving approach, checking the cross-over behavior specified in (3.0) is essential to ensure success.  The window for A sensitively varies with the change in the nonlinear differential equations, unfortunately.  For solving the nonlinear two-point boundary value problem, a general code as straightforward as the embedded function of NDSolve is possible if and only if a systematic way of searching for the appropriate window of A exists.

(4) Direct extension of the root-solving approach has been made to the two-point boundary value problems with coupled nonlinear differential equations.  For the boundary value problems involving two coupled nonlinear differential equations, in large, two initial slopes need to be assumed in NDSolve and there exists two boundary conditions to be satisfied at η = xend.  Interpolation needs to be replaced by ListInterpolation due to the expansion of the number of unknowns.  FindRoot, in addition, involves the solutions for two simultaneous equations from the two boundary conditions.  Various boundary value problems in forced and natural convective heat transfer have been studied, posted, and constantly updated in the author's web page (http://web.missouri.edu/~mae/).