Solving Boundary Value Problems for Ordinary Di erential ...

1 Solving Boundary Value Problems for OrdinaryDi erential Equations in Matlabwithbvp4cLawrence F. Shampine Jacek KierzenkayMark W. ReicheltzOctober 26, 20001 IntroductionOrdinary di erential equations (ODEs) describe phenomena that change contin-uously. They arise in models throughout mathematics, science, and itself, a system of ODEs has many solutions. Commonly a solution of inter-est is determined by specifying the values of all its components at a single pointx=a. This is an initial Value problem (IVP). However, in many applications asolution is determined in a more complicated way. A Boundary Value problem (BVP) speci es values or equations for solution components at more than onex. Unlike IVPs, a Boundary Value problem may not have a solution, or mayhave a nite number, or may have in nitely many. Because of this, programsfor Solving BVPs require users to provide a guess for the solution desired.

2 Of-ten there are parameters that have to be determined so that the BVP has asolution. Again there might be more than one possibility, so programs requirea guess for the parameters desired. Singularities in coe cients and problemsposed on in nite intervals are not unusual. Simple examples are used inx2toillustrate some of these tutorial shows how to formulate, solve, and plot the solution of a BVPwith theMatlabprogrambvp4c. It aims to make Solving a typical BVP aseasy as possible. BVPs are much harder to solve than IVPs and any solver mightfail, even with good guesses for the solution and unknown e ective solver, but the underlying method and the computing environmentare not appropriate for high accuracies nor for Problems with extremely sharpchanges in their solutions. Section 3 describes briefly the numerical 4 is a collection of examples that illustrate the solution of BVPs withbvp4c.

3 The rst three should be read in order because they introduce suc-cessively features of the solver as it is applied to typical Problems . Although Math. Dept., SMU, Dallas, TX 75275 MathWorks, Inc., 3 Apple Hill Drive, Natick, MA 01760 Coolidge Road, Wayland, MA 01778 quite general BVPs, Problems arise in the most diverse forms andthey may require some preparation for their solution. The remaining examplesillustrate this preparation for common tasks. Some exercises are included forpractice. M- les for the solution of all the examples and exercises accompanythis Boundary Value ProblemsIf the functionfis smooth on [a; b], the initial Value problemy0=f(x;y),y(a)given, has a solution, and only one. Two-point Boundary Value Problems areexempli ed by the equationy00+y=0(1)with Boundary conditionsy(a)=A,y(b)=B. An important way to analyzesuch Problems is to consider a family of solutions of IVPs.

4 Lety(x; s)bethesolution of equation (1) with initial valuesy(a)=A,y0(a)=s. Eachy(x; s)extends tox=band we ask, for what values ofsdoesy(b; s)=B?Ifthereis a solutionsto this algebraic equation, the correspondingy(x; s) provides asolution of the di erential equation that satis es the two Boundary linearity we can sort out the possibilities easily. Letu(x) be the solutionde ned byy(a)=A,y0(a)=0andv(x) be the solution de ned byy(a)=0,y0(a) = 1. Linearity implies thaty(x; s)=u(x)+sv(x), and the boundaryconditionB=y(b; s)=u(b)+sv(b) amounts to a linear algebraic equation forthe unknown initial slopes. The familiar facts of existence and uniqueness ofsolutions of linear algebraic equations then tell us that there is either exactlyone solution to the BVP, or there are Boundary valuesBfor which there is nosolution and others for which there are in nitely many Problems , more speci cally Sturm- liouville Problems , are exem-pli ed byy00+ y=0withy(0) = 0,y( ) = 0.

5 Such a problem obviously has the trivial solutiony(x) 0, but for some values of , there are non-trivial solutions. Such arecalled eigenvalues and the corresponding solutions are called eigenfunctions . Ify(x) is a solution of this BVP, it is obvious that y(x) is, too. Accordingly,we need a normalizing condition to specify a solution of interest, for instancey0(0) = 1. For >0, the solution of the IVP withy(0) = 0,y0(0) = 1 isy(x)=sin xp =p . The Boundary conditiony( ) = 0 amounts to a non-linear algebraic equation for . Generally existence and uniqueness of solutionsof nonlinear algebraic equations are di cult matters. For this example the al-gebraic equation is solved easily to nd that the BVP has a non-trivial solutionif, and only if, =k2fork=1;2;:::. This example shows that when Solving aSturm- liouville problem , we have to specify not only a normalizing condition,but also which eigenvalue interests 2 1 1: Two solutions fory00+jyj= introduces other complications illustrated by the problem [3]y00+jyj=0withy(0) = 0,y(b)=B.

6 Proceeding as with the linear examples, it is foundthat for anyb> , there are exactlytwosolutions for anyB<0. One solutionhas the formy(x; s)=ssinhx; it starts o with a negative slopesand decreasesmonotonely toB. The other starts o with a positive slope where it has the formy(x; s)=ssinx. This solution crosses the axis atx= , where its form changesand it decreases thereafter monotonely toB. Figure 1 shows an example ofthis withb=4andB= 2. Much as with eigenvalue Problems , when solvingnonlinear BVPs we have to specify which solution is the one that interests inx4 show that BVPs modelling physical situations do not nec-essarily have unique solutions. Other examples show that Problems involvingphysical parameters might have solutions only for parameter values in certainranges. The examples make it clear that in practice, Solving BVPs may wellinvolve an exploration of the existence and uniqueness of solutions of a is quite di erent from Solving Numerical MethodsThe theoretical approach to BVPs ofx2 is based on the solution of IVPs forODEs and the solution of nonlinear algebraic equations.

7 Because there aree ective programs for both tasks, it is natural to combine them in a programfor the solution of BVPs. The approach is called a shooting method. Because3it appears so straightforward to use quality numerical tools for the solution ofBVPs by shooting, it is perhaps surprising thatbvp4cisnota shooting basic di culty with shooting is that a perfectly nice BVP can require theintegration of IVPs that are unstable. That is, the solution of a BVP canbe insensitive to changes in Boundary values, yet the solutions of the IVPs ofshooting are sensitive to changes in initial values. The simple exampley00 100y=0withy(0) = 1,y(1) =Bmakes the point. Shooting involves the solutiony(x; s)=cosh10x+0:1ssinh 10xof the IVP with initial valuesy(0) = 1,y0(0) =s. Obviously@y=@s=0:1 sinh 10x, which can be as large as 0:1 sinh 10 1101. A little calculation shows that the slope that results in satisfaction ofthe Boundary condition atx=1iss= 10(B cosh 10)=sinh 10 and then thatfor the solution of the BVP,j@y=@Bj=jsinh 10x=sinh 10j 1.

8 Evidently thesolutions of the IVPs are considerably more sensitive to changes in the initialslopesthan the solution of the BVP is to changes in the Boundary IVPs are not too unstable, shooting can be quite e ective. Unstable IVPscan cause a shooting code to fail because the integration \blows up" beforereaching the end of the interval. More often, though, the IVP solver reachesthe end, but is unable to compute an accurate result there and because of this,the nonlinear equation solver is unable to nd accurate initial values. A varietyof techniques are employed to improve shooting, but when the IVPs are veryunstable, shooting is just not a natural approach to Solving a collocation method for the solution of BVPs of the formy0=f(x;y;p);a x bsubject to general nonlinear, two-point Boundary conditionsg(y(a);y(b);p)=0 Herepis a vector of unknown parameters. For simplicity it is suppressed inthe expressions that follow.

9 The approximate solutionS(x) is a continuousfunction that is a cubic polynomial on each subinterval [xn;xn+1]ofamesha=x0<x1< ::: < xN=b. It satis es the Boundary conditionsg(S(a);S(b)) =0and it satis es the di erential equations (collocates) at both ends and the mid-point of each subintervalS0(xn)=f(xn;S(xn))S0((xn+xn+1 )=2) =f((xn+xn+1)=2;S((xn+xn+1)=2))S0(xn+1)=f (xn+1;S(xn+1))These conditions result in a system of nonlinear algebraic equations for thecoe cients de ningS(x). In contrast to shooting, the solutiony(x) is approxi-matedoverthewholeinterval[a; b] and the Boundary conditions are taken into4account at all times. The nonlinear algebraic equations are solved iteratively bylinearization, so this approach relies upon the linear equation solvers ofMatlabrather than its IVP codes. The basic method ofbvp4c, which we call Simpson'smethod, is well-known and is found in a number of codes. It can be shown[8] that with modest assumptions,S(x) is a fourth order approximation to anisolated solutiony(x), ,ky(x) S(x)k the maximum ofthe step sizeshn=xn+1 xnandCis a constant.

10 Because it is not true ofsome popular collocation methods, we stress the important fact that this boundholds for allxin [a; b]. AfterS(x) is computed on a mesh withbvp4c,itcanbeevaluated inexpensively at anyx,orsetofx,in[a; b] with BVPs can have more than one solution, BVP codes require usersto supply a guess for the solution desired. The guess includes a guess for aninitial mesh that reveals the behavior of the desired solution. The codes thenadapt the mesh so as to obtain an accurate numerical solution with a modestnumber of mesh points. Coming up with a su ciently good guess is often thehardest part of Solving a an unusual approach to the controlof error that helps it deal with poor guesses. The continuity ofS(x)on[a; b]andcollocation at the ends of each subinterval imply thatS(x) also has a continuousderivative on [a; b]. For such an approximation, the residualr(x)intheODEsis de ned byr(x)=S0(x) f(x;S(x))Put di erently, this says thatS(x) is the exact solution of the perturbed ODEsS0(x)=f(x;S(x)) +r(x)Similarly, the residual in the Boundary conditions isg(S(a);S(b)).