The Shooting Method for the Solutionof Two-Point Boundary Value ProblemsDouglas B. MeadeDepartment of MathematicsUniversity of South CarolinaColumbia, SC 29208, USA;E-mail: meade@math.sc.edu WWW URL: http://www.math.sc.edu/~meade/ Copyright 2004: Douglas B. Meade

One of the strengths of Maple is its ability to provide a wide variety of information about solutions to differential equations. Explicit, implicit, parametric, series, Laplace transform, numerical, and graphical solutions can all be obtained via the dsolve command. Numerical solutions are of particular interest due to the fact that exact solutions do not exist, in closed form, for most engineering and scientific applications. The numerical solution methods available within dsolve are applicable only to initial value problems. Thus, at first glance, Maple appears to be very limited in its ability to analyze the multitude of two-point boundary value problems that occur frequently in engineering analysis.A commonly used numerical method for the solution of two-point boundary value problems is the shooting method. This well-known technique is an iterative algorithm which attempts to identify appropriate initial conditions for a related initial value problem (IVP) that provides the solution to the original boundary value problem (BVP).The first objective of this paper is to describe the shooting method and its Maple implementation, shoot. Then, shoot is used to analyze some common two-point BVPs from chemical engineering:the Blasius solution for laminar boundary-layer flow past a flat plate,the reactivity behavior of porous catalyst particles subject to both internal mass concentration gradients and temperature gradients, andthe steady-state flow near an infinite rotating disk.

This worksheet was originally published in the Maple Technical Newsletter. The citation for that article is: Douglas B. Meade, Bala S. Haran, and Ralph E. White, The shooting technique for the solution of two-point boundary value problems, MapleTech, 3(1) 1996, pp. 85-93.While the applications are description of the problems is nearly the same as in the original paper, the Shoot package has been almost completely rewritten and updated for Maple 9.5.This worksheet, and shoot, were originally created in Maple V Release 3. The update to Release 5, and Maple 6, has required a complete reimplementation of shoot. The worksheet, for Maple 9.5, is also significantly modified in some places. However, the content and functionality are as close as possible to that of the original.Please let me know if you have any problems or questions. I also would like to hear about any successes that you have using shoot.Thank you!Douglas B. Meademeade@math.sc.edu December 2004

restart;ShootLib := "C:\\Documents and Settings/DMeade/Desktop/Shoot9/":libname := ShootLib, libname:with( Shoot );with( plots ):

The basic idea of the shooting method for two-point boundary value problems is to reformulate the problem as a nonlinear parameter estimation problem. The new problem requires the solution of a related initial value problem (IVP) with initial conditions chosen to approximate the boundary conditions at the other endpoint. If these boundary conditions are not satisfied to the desired accuracy, the process is repeated with a new set of initial conditions until the desired accuracy is achieved or an iteration limit is reached.To be more specific, consider the two-point BVP for a coupled system of NiMlIm5H first-order ODEs NiMvLSUlZGlmZkc2JC0lInlHNiMlInRHRiotJSJmRzYkRipGJw== NiMvLSYlInlHNiMlImlHNiMlImFHJiUmYWxwaGFHRic= , NiMvJSJpRzsiIiImJSJtRzYjRiY= (1) NiMvLSYlInlHNiMsJiYlIm1HNiMiIiJGLCUiakdGLDYjJSJiRyYlJWJldGFHNiNGLQ== , NiMvJSJqRzsiIiImJSJtRzYjIiIj The vector NiMlInlH contains the NiMlIm5H unknown functions of the independent variable NiMlInRH. The unknown functions are ordered so that the first NiMmJSJtRzYjIiIi (0<NiMmJSJtRzYjIiIi<NiMlIm5H) components of NiMlInlH have boundary conditions at NiMvJSJ0RyUiYUc=. The remaining NiMmJSJtRzYjIiIj := NiMsJiUibkciIiImJSJtRzYjRiUhIiI= components of the solution have boundary conditions specified at a second point, NiMvJSJ0RyUiYkc=. The Maple procedure shoot supports nonlinear boundary conditions at NiMvJSJ0RyUiYkc=. In particular, NiMmSSQoMSlHNiI2IyIiIw== and NiMmSSQoMSlHNiI2IyIiJA== can be replaced with NiMvLSZJInJHNiI2I0kiaUdGJzYjLUkieUdGJzYjSSJhR0YnIiIh and NiMvLSYlInJHNiMlImpHNiMtJSJ5RzYjJSJiRyIiIQ==, where each NiNJInJHNiI= is a differentiable function of NiNJIm5HNiI= variables. The examples discussed in this paper all use NiMtJiUickc2IyUiakc2Iy0lInlHNiMlImJH := NiMsJi0mJSJ5RzYjLCYmJSJtRzYjIiIiRiwlImpHRiw2IyUiYkdGLCYlJWJldGFHNiNGLSEiIg==. (Note that if NiMvJiUibUc2IyIiIyIiIQ==, then NiNJJCgxKUc2Ig== is an initial value problem.)The shooting method seeks to identify a vector of parameters NiMlInNH in NiMpJSJSRyYlIm1HNiMiIiM=so that the solution, denoted by NiMtJSJ5RzYkJSJ0RyUic0c=, to the initial value problem NiMvLSUlZGlmZkc2JCUieUclInRHLSUiZkc2JEYoLUYnNiRGKCUic0c= NiMvLSYlInlHNiMlImlHNiQlImFHJSJzRyYlJmFscGhhR0Yn , NiMvJSJpRzsiIiImJSJtRzYjRiY= (2) NiMvLSYlInlHNiMsJiYlIm1HNiMiIiJGLCUiakdGLDYkJSJhRyUic0cmRjA2I0Yt , NiMvJSJqRzsiIiImJSJtRzYjIiIjagrees with the solution to (1). Note that (2) is simply (1) with the boundary conditions at NiMvJSJ0RyUiYkc= replaced with unknown initial conditions at NiMvJSJ0RyUiYUc=. To determine the correct initial values, consider the ``objective function'' NiMlIkZH with components NiMtJiUiRkc2IyUiakc2IyUic0c= := NiMsJi0mJSJ5RzYjLCYmJSJtRzYjIiIiRiwlImpHRiw2JCUiYkclInNHRiwmJSViZXRhRzYjRi0hIiI=, NiMvJSJqRzsiIiImJSJtRzYjIiIj .Then (1) is solvable if and only if there exists NiMlInNH in NiMpJSJSRyYlIm1HNiMiIiM= so that NiMvLSUiRkc2IyUic0ciIiE=.The success of this process depends primarily on the iterative procedure used to construct a sequence of parameter vectors that converges to a zero of NiMlIkZH. While any numerical root-finding algorithm could be employed for this step, one step of the Newton-Raphson method is most commonly used. That is, given an initial guess NiMqJCUic0ciIiE= in NiMpJSJSRyYlIm1HNiMiIiM=, define a sequence of initial conditions {NiMpJSJzRyUia0c=} by NiMpJSJzRywmJSJrRyIiIkYnRic= := NiMsJiklInNHJSJrRyIiIiomKS0lI0pGRzYjRiQsJEYnISIiRictJSJGR0YsRidGLg==for all NiMlImtH >= 0. To implement this, note that the vector NiMtJSJGRzYjKSUic0clImtH is directly available from the solution of (2), but the Jacobian matrix NiMtJSNKRkc2IyklInNHJSJrRw== requires the values of NiMtJSVkaWZmRzYkLSYlInlHNiMsJiYlIm1HNiMiIiJGLiUiaUdGLjYkJSJiRyklInNHJSJrRyZGMzYjJSJqRw== for all NiQlImlHLyUiakc7IiIiJiUibUc2IyIiIw==. These values can be obtained by solving the NiMlIm5H IVPs in (2) together with the NiMqJiUibkciIiImJSJtRzYjIiIjRiU= sensitivity equations ([2, p. 226], [3, pp. 54--58]):NiMvLSUlZGlmZkc2JC1GJTYkLSYlInlHNiMlImlHNiMlInRHJiUic0c2IyUiakdGLyomLUYlNiQlImZHRioiIiItRiU2JEYqRjBGOA==, NiMvJSJpRzsiIiIlIm5H, NiMvJSJqRzsiIiImJSJtRzYjIiIj with corresponding initial conditions NiMvLSUlZGlmZkc2JCYlInlHNiMlImlHJiUic0c2IyUiakciIiE= , for all NiMvJSJpRzsiIiImJSJtRzYjRiY=, NiMvJSJqRzsiIiImJSJtRzYjIiIj NiMvLSUlZGlmZkc2JCYlInlHNiMsJiYlIm1HNiMiIiJGLiUiakdGLiYlInNHNiNGLyYlJmRlbHRhRzYjKiYlImlHRi5GL0Yu, for all NiQlImlHLyUiakc7IiIiJiUibUc2IyIiIw== The above algorithm is known as the simple, or single, shooting method. While this method is effective for many problems, there are some problems that are encountered in practice [4]. Assume (1) has a unique solution. There is no guarantee that the initial value problem (2) has a solution on the interval NiM3JEkiYUc2IkkiYkdGJQ== for all NiNJInNHNiI= in NiMpJSJSRyYlIm1HNiMiIiM= . Even if (2) does have a solution on NiM3JEkiYUc2IkkiYkdGJQ== the problem may be stiff. In such a case the solution at NiMvJSJ0RyUiYkc= may be so inaccurate as to make the results of the Newton-Raphson step meaningless. When the solution at NiMvJSJ0RyUiYkc= is known with sufficient accuracy, the local convergence of the Newton-Raphson step may prevent the iterations from converging to a solution of the original boundary value problem (1). This difficulty can be addressed by replacing the Newton-Raphson step with another iterative solver with improved convergence properties (e.g., modified Newton's method and Broyden's method). The multiple, or parallel, shooting method addresses the other difficulties. As these problems are generally more pronounced as NiMsJiUiYkciIiIlImFHISIi increases, it seems appropriate to consider partitioning the interval into NiMlIk5H subintervals. Then, using compatibility conditions between the subintervals, a well-posed IVP is obtained on each subinterval.The combination of Maple's symbolic and numerical facilities for manipulating and solving IVPs provide an excellent environment for the translation of the simple shooting method into the Maple programming language. The syntax for this procedure, called shoot, closely follows the syntax of dsolve. In particular, the data structures returned by shoot are identical to the ones returned by dsolve/numeric. Thus, all the techniques and tools for manipulating Maple's numerical solutions of differential equations, e.g. plots[odeplot], can be used to interpret the results from shoot. (For convenience, shoot is stored as a Maple repository (i.e., package). This repository can be obtained from the author's website at http://www.math.sc.edu/~meade/maple/Shoot9/Shoot9.zip . The source code for shoot, including on-line help documentation, is available upon request from the author.)Other numerical methods for the solution of two-point boundary value problems have also been developed on other platforms. One of the most well-known is the fortran subroutine COLSYS [4] (available on the WWW from Netlib). The numerical method used in this routine is collocation of B-splines at Gaussian points [5]. Finite difference methods can be easily implemented using Matlab's fsolve command.The examples discussed in this paper provide samples of how simple shooting, in particular shoot, can be applied to the analysis of boundary value problems that are encountered in chemical engineering.

Prior to showing the use of shoot for realistic problems in chemical engineering, let's start with a relatively standard problem from elementary ODEs: NiMvLCYtSSVkaWZmR0kqcHJvdGVjdGVkR0YnNiUtSSJZRzYiNiNJInRHRitGLUYtIiIiRilGLi1JJGNvc0c2JEYnSShfc3lzbGliR0YrRiw= NiMvLUkiWUc2IjYjIiIiIiIk (3) NiMvLUkiWUc2IjYjIiIkIiIi This problem is solvable directly with Maple's built-in dsolve command:dsolve( {diff( Y(t),t,t ) + Y(t) = cos(t), Y(1)=3, Y(3)=1}, Y(t) );Sexact := unapply( rhs(%), t ):This was not the case when shoot was first written (in 1996 for Maple V, Release 3).To solve the same problem with the shooting method, first restate the problem as a first-order systemODE:={diff(Y(t),t) = Yp(t), diff(Yp(t),t) = -Y(t)+cos(t)}:for the function, Y, and its first derivative, Yp = Y' FNS:={ Y(t), Yp(t) }:with initial conditionsIC:={ Y(1)=3, Yp(1)=alpha }:and boundary conditionsBC:={ Y(3)=1 }:Note that the initial condition for the auxiliary function Yp reflects that alpha is the shooting parameter for this problem. This can be handled by shoot without modification:infolevel[shoot]:=1:S:=shoot( ODE, IC, BC, FNS, [alpha=0.0] ):S(3);evalf( Sexact(3) );Pshoot := odeplot( S, [t,Y(t)], t=1..3, legend=["Y(t) [shoot]"] ):Pexact := plot( Sexact(t), t=1..3, style=point, color=blue, legend=["Y(t) [exact]"] ):display( [Pshoot,Pexact], view=[1..3,0..5], title="Comparison of shooting method and exact solutions" );

Consider a fluid stream with velocity NiMmJSJ1RzYjIiIh and kinematic viscosity NiMlI251Rw== in which a thin plate is inserted parallel with the fluid flow. Determining the velocity of the fluid in the region close to the plate is the Blasius problem [6, p. 233]. Assuming the flow is steady, incompressible, and Newtonian, the plate is infinitely wide, and neglecting buoyancy, the equations of motion and continuity are:alias( U=u(x,y), V=v(x,y) ):PDE:={ U*diff(U,x)+V*diff(U,y)-nu*diff(U,y$2)=0, diff(U,x)+diff(V,y)=0 };where NiMlInVH and NiMlInZH are the NiMlInhH- and NiMlInlH-components of the fluid velocity. The boundary conditions consist of the ``no-slip'' conditions on the boundary of the plate: NiMtJSJ1RzYkJSJ4RyIiIQ===NiMtJSJ2RzYkJSJ4RyIiIQ===0, and the free stream-merge condition NiMvLSUmbGltaXRHNiQtJSJ1RzYkJSJ4RyUieUcvRislKWluZmluaXR5RyZGKDYjIiIh , for all NiMlInhH>0.A similarity transformation can be used to reduce this parabolic system of PDEs to a single ODE. This can be done by choosing the dimensionless similarity variable to be:simsubs:=eta(x,y)=y*sqrt(u[0]/nu/x/2);The corresponding nondimensional stream function for the flow is:stream:=psi(x,y)=sqrt(2*nu*x*u[0])*f(eta(x,y));Thus, the velocities can be expressed as:Usubs:=U=diff(rhs(stream),y);Vsubs:=V=-diff(rhs(stream),x);Substituting the stream function representations of the velocities into the PDEs is tedious to complete by hand. Fortunately, this is exactly one of Maple's strengths:ODE:=simplify(subs(Usubs,Vsubs,simsubs,PDE));The trivial fulfillment of the continuity equation is evident in this result. However, the first equation is not so readily identified. To simplify this further, note that each argument to NiMlImZH, and its derivatives, is simply NiMlJGV0YUc=. To force this simplification,simsubs2:=solve(subs(eta(x,y)=eta,simsubs),{y});ODE:=simplify(subs(simsubs2,ODE),symbolic);It is now easy to identify the Blasius equation NiMvLCYtJSVkaWZmRzYkLSUiZkc2IyUkZXRhRy0lIiRHNiRGKyIiJCIiIiomRihGMC1GJjYkRigtRi02JEYrIiIjRjBGMCIiIQ==, NiMlJGV0YUc=>0.The conversion of the boundary conditions can be done by inspection. The resulting conditions are NiMtJSJmRzYjIiIh = NiMqJiUiZkciIiItJSInRzYjIiIhRiU= = NiMiIiE= and NiMvLSUmbGltaXRHNiQtJSVkaWZmRzYkLSUiZkc2IyUkZXRhR0YtL0YtJSlpbmZpbml0eUciIiI=.The ``boundary condition'' at NiMvJSRldGFHJSlpbmZpbml0eUc= presents a problem. However, a simple asymptotic analysis shows that solutions at NiMlJGV0YUc==10 are safely in the far-field. The problem now takes the form of a third-order two-point boundary value problem for NiMlImZH on [0,10]. The reformulation as a first-order system is straightforward; let NiM+JSJnRyomJSJmRyIiIiUiJ0dGJw== and NiM+JSJoRyomJSJmRyIiIiUjJydHRic=, then NiMvKiYlImZHIiIiJSInR0YmJSJnRw== , NiMvLSUiZkc2IyIiIUYn NiMvKiYlImdHIiIiJSInR0YmJSJoRw==, NiMvLSUiZ0c2IyIiIUYn (4) NiMvKiYlImhHIiIiJSInR0YmLCQqJiUiZkdGJkYlRiYhIiI=, NiMvLSUiaEc2IyIiISUlYmV0YUc= where NiMlJWJldGFH is the control parameter. The objective is to find NiMlJWJldGFH so that the solution to (4) satisfies the boundary condition NiMvLSUiZ0c2IyIjNSIiIg==. Since there is only one boundary condition at NiMvJSRldGFHIiM1 we have NiMvJiUibUc2IyIiIyIiIg== and NiMvJSJzRyUlYmV0YUc=. The sensitivity equations that are obtained from (3), together with their accompanying initial values, are: NiMvKiYmJSJmRzYjJSViZXRhRyIiIiUiJ0dGKSYlImdHRic= , NiMvLSYlImZHNiMlJWJldGFHNiMiIiFGKg== NiMvKiYmJSJnRzYjJSViZXRhRyIiIiUiJ0dGKSYlImhHRic= , NiMvLSYlImdHNiMlJWJldGFHNiMiIiFGKg== (5) NiMvKiYmJSJoRzYjJSViZXRhRyIiIiUiJ0dGKSwmKiYmJSJmR0YnRilGJkYpISIiKiZGLkYpRiVGKUYv , NiMvLSYlImhHNiMlJWJldGFHNiMiIiEiIiI= where, for example, NiM+LSYlImZHNiMlJWJldGFHNiQlInRHRigtJSVkaWZmRzYkLUYmRilGKA==. Now, assume that an approximate solution to (4) and (5) has been computed with NiMvJSJzRylGJCUia0c=, i.e., NiMvJSViZXRhRylGJCUia0c=. Then, assuming NiMtJSRhYnNHNiMsJi0lImdHNiQiIzUpJSViZXRhRyUia0ciIiJGLiEiIg== is not sufficiently small, the next iteration will be made with NiM+KSUlYmV0YUcsJiUia0ciIiJGKEYoLCYpRiVGJ0YoKiYsJi0lImdHNiQiIzVGKkYoRighIiJGKC0mRi42I0YlRi9GMUYx. While it is important to understand how shoot obtains its approximations, it is also a tremendous advantage to use Maple to both compute automatically, and symbolically, the sensitivity equations and iteratively solve the combined system of IVPs and take one step of the Newton-Raphson method until the boundary conditions are satisfied within a prescribed tolerance.Now, define the dependent variables for the first-order system:FNS:={ f(eta), g(eta), h(eta) }:and the differential equations which they satisfy:ODE:={ diff(f(eta),eta)=g(eta), diff(g(eta),eta)=h(eta), diff(h(eta),eta)=-f(eta)*h(eta) }:The initial condition for NiMqJiUiZkciIiIlIycnR0Yl is unknown, and the boundary condition at NiMvJSRldGFHIiM1 is the control for our iterations:IC:={ f(0)=0, g(0)=0, h(0)=beta }:BC:=g(10)=1:We take the first shot with NiMvJSViZXRhRyIiIQ==, and iterate until the boundary condition is satisfied to six decimal places:infolevel[shoot]:=1;S:=shoot( ODE, IC, BC, FNS, beta=0, abserr=Float(5,-7), output=listprocedure, method=taylorseries ):This solution can be analyzed using any of the standard Maple tools for manipulating the numerical solution of a differential equation. A common place to begin is a plot of the solution, including labels to distinguish the different curves:odeplot( S, [ [eta,f(eta)], [eta,g(eta)], [eta,h(eta)] ],0..4, labels=[`eta`,``], legend=["f","f '","f ''"], title="Blasius solution for flat-plate boundary layer" );With increasing distance from the leading edge of the plate in the downstream direction the thickness NiMlJmRlbHRhRw== of the retarded boundary layer increases continuously as increasing quantities of fluid become affected. In the boundary layer the velocity of the fluid increases from zero at the wall (no slip) to its full value which corresponds to external frictionless flow. The velocity reaches 99% of its bulk value when NiMvKiYlImZHIiIiLSUiJ0c2IyUkZXRhR0YmJCIjKiohIiM=:fp0:=subs( S, g(eta) ):fp := proc(x) if not type(evalf(x),`numeric`) then 'procname'(x); else fp0(x); end if;end proc:eta[`99%`]=fsolve( fp(eta)=0.99, eta=3..4 );The solution returned by shoot is obtained with method=taylorseries and the function used in fsolve must be written to avoid evaluation in cases when the argument is not numeric. But, in the end, the result is equivalent to the one reported in the original paper.The corresponding boundary layer thickness is:delta[`99%`]=subs(eta=rhs(%), combine(rhs(op(simsubs2))) );As is now evident, the thickness of the boundary layer decreases with decreasing viscosity.The shear stress,NiMvJSR0YXVHKiYlI211RyIiIi0lJWRpZmZHNiQlInVHJSJ5R0Yn, is:tau:=mu*simplify(diff(subs(simsubs,rhs(Usubs)),y));Note that a large velocity gradient across the flow creates considerable shear stress in the boundary layer. The wall shear stress is:tau[w]:=subs( (D@@2)(f)(0)=subs(S,h(eta))(0), subs( y=0, tau ) );The analysis can be continued to obtain an understanding of parameters such as the displacement thickness and drag coefficient.

Consider the steady fluid flow generated when the infinite plane NiMvJSJ6RyIiIQ==, immersed in a Newtonian viscous fluid, rotates about the axis NiMvJSJyRyIiIQ== with a constant angular velocity NiMlJm9tZWdhRw==. The viscous drag of the rotating surface creates a swirling flow toward the disk. The motion is characterized in terms of the pressure, NiMlInBH, and the three components of the velocity, NiMpJSJ2RyUickc=, NiMpJSJ2RyUmdGhldGFH, NiMpJSJ2RyUiekc=, in cylindrical coordinates. The radial symmetry of this problem eliminates NiMlJnRoZXRhRw== as a independent variable. Thus, with NiMlJHJob0c= and NiMlI251Rw== denoting the density and kinematic viscosity of the fluid and writing partial derivatives as subscripts, the equations of continuity and conservation of momentum reduce to [6]: NiMqJiIiIkYkJSJyRyEiIg== NiMmKiYlInJHIiIiKSUidkdGJUYmNiNGJQ== + NiMmKSUidkclInpHNiNGJg== = NiMiIiE=, NiMsJiomKSUidkclInJHIiIiJkYlNiNGJ0YoRigmLSZGJjYjJSJ6RzYjRiVGLkYo - NiMqJiIiIkYkJSJyRyEiIg== NiMqJCklInZHJSZ0aGV0YUciIiM= = NiMsJComIiIiRiUlJHJob0chIiJGJw== NiMmJSJwRzYjJSJyRw== + NiMlI251Rw== ( NiMmKSUidkclInJHNiMlI3JyRw== + NiMqJiIiIkYkJSJyRyEiIg== NiMmKSUidkclInJHNiNGJg== + NiMmKSUidkclInJHNiMlI3p6Rw== - NiMqJiklInZHJSJyRyIiIiokRiYiIiMhIiI= ) (6) NiMsJiomKSUidkclInJHIiIiJilGJiUmdGhldGFHNiNGJ0YoRigqJilGJiUiekdGKCZGKjYjRi9GKEYo + NiMqJiIiIkYkJSJyRyEiIg== NiMqJiklInZHJSJyRyIiIilGJSUmdGhldGFHRic= = NiMlI251Rw== ( NiMmKSUidkclJnRoZXRhRzYjJSNyckc= + NiMqJiIiIkYkJSJyRyEiIg== NiMmKSUidkclJnRoZXRhRzYjJSJyRw== + NiMmKSUidkclJnRoZXRhRzYjJSN6ekc= - NiMqJiklInZHJSZ0aGV0YUciIiIqJCUickciIiMhIiI= ) NiMsJiomKSUidkclInJHIiIiJilGJiUiekc2I0YnRihGKComRipGKCZGKjYjRitGKEYo = NiMsJComIiIiRiUlJHJob0chIiJGJw== NiMmJSJwRzYjJSJ6Rw== + NiMlI251Rw== ( NiMmKSUidkclInpHNiMlI3JyRw== + NiMqJiIiIkYkJSJyRyEiIg== NiMmKSUidkclInpHNiMlInJH + NiMmKSUidkclInpHNiMlI3p6Rw== ) The boundary conditions are chosen (for all NiMlInJH>0) to enforce no slippage at the interface with the disk: NiMtKSUidkclInJHNiRGJiIiIQ== = NiMtKSUidkclInpHNiQlInJHIiIh = 0, NiMvLSklInZHJSZ0aGV0YUc2JCUickciIiEqJkYpIiIiJSZvbWVnYUdGLA==, (7) NiMvLSUicEc2JCUickciIiFGKA== and no non-axial viscous effect in the far-field: NiMvLSUmbGltaXRHNiQtKSUidkclInJHNiRGKiUiekcvRiwlKWluZmluaXR5RyIiIQ==, NiMvLSUmbGltaXRHNiQtKSUidkclJnRoZXRhRzYkJSJyRyUiekcvRi0lKWluZmluaXR5RyIiIQ==, (8) NiMvLSUmbGltaXRHNiQtJiklInZHJSJ6RzYjRis2JCUickdGKy9GKyUpaW5maW5pdHlHIiIh. Similarity solutions to this equation were found by Karman [9], who noted that each of NiMqJiklInZHJSJyRyIiIkYmISIi, NiMqJiklInZHJSZ0aGV0YUciIiIlInJHISIi, NiMqJiklInZHJSJ6RyIiIiUickchIiI=, and NiMlInBH depends only on the distance from the disk, NiMlInpH. The system of PDEs reduces to a system of ODEs with the introduction of NiMlInpH*=NiMqJiUiekciIiItJSVzcXJ0RzYjKiYlJm9tZWdhR0YlJSNudUchIiJGJQ== as the dimensionless independent variable for the dimensionless functions NiMlIkZH, NiMlIkdH, NiMlIkhH, and NiMlIlBH defined according to NiMvKSUidkclInJHKihGJiIiIiUmb21lZ2FHRigtJSJGRzYjKiYlInpHRiglIipHRihGKA==, NiMvKSUidkclJnRoZXRhRyooJSJyRyIiIiUmb21lZ2FHRiktJSJHRzYjKiYlInpHRiklIipHRilGKQ==, (9) NiMvKSUidkclInpHKiYtJSVzcXJ0RzYjKiYlJm9tZWdhRyIiIiUjbnVHRi1GLS0lIkhHNiMqJkYmRi0lIipHRi1GLQ==, NiMvJSJwRyoqJSRyaG9HIiIiJSZvbWVnYUdGJyUjbnVHRictJSJQRzYjKiYlInpHRiclIipHRidGJw==, The substitution of (9) into (6), (7), and (8) can be completed in a manner analogous to that used in the Blasius problem. These steps are omitted here in the interest of space; the resulting BVP is: NiMvKiYlIkhHIiIiJSInR0YmLCQqJiIiI0YmJSJGR0YmISIi, NiMvKiYlIkZHIiIiJSMnJ0dGJiwoKiRGJSIiI0YmKiQlIkdHRiohIiIqKEYlRiYlIidHRiYlIkhHRiZGJg==, NiMvKiYlIkdHIiIiJSMnJ0dGJiwmKigiIiNGJiUiRkdGJkYlRiZGJiooJSJIR0YmRiVGJiUiJ0dGJkYm NiMvKiYlIlBHIiIiJSInR0YmLCYqKCUiSEdGJkYqRiZGJ0YmISIiKiZGKkYmJSMnJ0dGJkYm, (10) NiMvLSUiRkc2IyIiIS0lIkhHRiY= = NiMvLSUiUEc2IyIiIUYn, NiMvLSUiR0c2IyIiISIiIg==, NiMtJSZsaW1pdEc2JC0lIkZHNiMqJiUiekciIiIlIipHRisvRiklKWluZmluaXR5Rw== = NiMvLSUmbGltaXRHNiQtJSJHRzYjKiYlInpHIiIiJSIqR0YsL0YqJSlpbmZpbml0eUciIiE=. Note that (10) can be solved as a two-point BVP for NiMlIkZH, NiMlIkdH, NiMlIkhH. Once this solution is known, the pressure equation can be integrated to yield NiMlIlBH=-NiMqJiIiIkYkIiIjISIi NiMqJCUiSEciIiM= + NiMlIkhH'. In addition, note that any physically realistic solution must have NiMlIkZH, NiMlIkdH, NiMlIlBH > 0 and NiMlIkhH < 0 for all NiMlInpH* [6, p. 164].Assuming the boundary condition at z* = NiMlKWluZmluaXR5Rw== can be applied at a finite distance from the disk, shoot can be used to obtain an approximate solution to this system. First, introduce the first derivatives of NiMlIkZH and NiMlIkdH as new dependent variables:FNS:={ F(Z), G(Z), H(Z), Fp(Z), Gp(Z) }:and reformulate the system as a first-order system of ODEs:ODE:={diff(H(Z),Z) = -2*F(Z), diff(F(Z),Z) = Fp(Z), diff(Fp(Z),Z) = -G(Z)^2+F(Z)^2+Fp(Z)*H(Z), diff(G(Z),Z) = Gp(Z), diff(Gp(Z),Z) = 2*F(Z)*G(Z)+H(Z)*Gp(Z)}:with initial conditions:IC:={ F(0)=0, G(0)=1, H(0)=0, Fp(0)=alpha, Gp(0)=beta }:and boundary conditions (again, NiMlInpH*=10 can be shown to be in the far-field):BC:={ F(10)=0, G(10)=0 }:Note that, as in the first example, the boundary condition at infinity has been moved to a finite position. The value NiMlInpH*=10 is chosen, as before, on the basis that this is already in the far field. (In fact, truncating the computations at NiMlInpH*=7 yields essentially the same solution.) Note also that since there are two boundary conditions at the second boundary point, there are two parameters to be determined in the shooting method (i.e., NiMvJiUibUc2IyIiI0Yn). This can be handled by shoot without modification:infolevel[shoot]:=1:S:=shoot( ODE, IC, BC, FNS, [alpha=0.51, beta=-0.62] ):P:=Z->-H(Z)^2/2-2*F(Z):odeplot(S,[ [Z,F(Z)], [Z,G(Z)], [Z,H(Z)], [Z,-P(Z)] ], 0..10, labels=[`z*`,``], title=`Infinite Disk: alpha=0.51, beta=-0.62`, legend=["F","G","H","-1/2*H^2-2*F"] );This plot confirms that the solution satisfies the original boundary conditions at NiMlInpH*=NiMlKWluZmluaXR5Rw==. (In fact, the solution obtained when the computational domain is truncated at NiMlInpH*=7.) Compare the solution in the preceding plot with the solution obtained when the starting values of the control parameters are NiMvJSZhbHBoYUckIiNdISIj and NiMvJSViZXRhRywkJCIjaCEiIyEiIg== -- each just 0.01 less than the first attempt:S2:=shoot( ODE, IC, BC, FNS, [alpha=0.5, beta=-0.61] ):The different values of NiMlJmFscGhhRw== and NiMlJWJldGFH suggest that the shooting method has returned a different solution. Is this also a solution to the problem?odeplot(S2,[ [Z,F(Z)], [Z,G(Z)], [Z,H(Z)], [Z,-P(Z)] ], 0..10, labels=[`z*`,``], title=`Infinite Disk: alpha=0.51, beta=-0.62`, legend=["F","G","H","-1/2*H^2-2*F"] );Note that this solution does not satisfy the aforementioned sign constraints for NiMlIkZH, NiMlIkdH, NiMlIkhH, and NiMlIlBH for a physical solution to (10). Furthermore, while the boundary conditions at NiMlInpH*=10 are satisfied, it is extremely unlikely that this solution will satisfy the boundary conditions at NiMlInpH*=NiMlKWluZmluaXR5Rw==. This is a spurious solution [6, p. 167]. Automating the detection of spurious solutions is extremely difficult. In some instances, this difficulty can be avoided by the use of one of the other solution techniques for boundary value problems mentioned previously.

In this brief article we have introduced shoot, a Maple implementation of the simple shooting method for the numerical solution of two-point boundary value problems, and illustrated the application of shoot to assist with the analysis of some classic problems from chemical engineering. These examples demonstrate some of the potential that is available with Maple's combination of symbolic, numeric, and graphic computation.

[1] J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, Springer-Verlag, 1980.[2] Yonathan Bard, Nonlinear Parameter Estimation, Academic Press, 1974.[3] Mark E. Davis, Numerical Methods and Modeling for Chemical Engineers, John Wiley & Sons, Inc., 1984.[3] Uri M. Ascher, Robert M. M. Mattheij, and Robert D. Russell, Numerical Solution of Boundary Value Problems for Ordinary Differential Equations}, SIAM, 1995.[4] U. Ascher, J. Christiansen and R. D. Russell, COLSYS - a collocation code for boundary value problems, Proc. Conf. for Codes for BVPs in ODEs, Houston, Texas, 1978.[5] Frank M. White, Viscous Fluid Flow, 1st ed., McGraw--Hill, Inc., 1974.[6] P. B. Wiesz and J. S. Hicks, The behaviour of porous catalyst particles in view of internal mass and heat diffusion effects}, Chemical Engineering Science, 17, pp. 265--275, 1962.[7] Bruce A. Finlayson, Nonlinear Analysis in Chemical Engineering, McGraw-Hill, 1980.[8] T. von Karman, `Uber laminare und turbulente Reibung, Z. Angew. Math. Mech., 1, pp. 233--252, 1921.

Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.