This discussion closely follows the reference Numerical Recipes by Press et al.
Suppose that a given n-th order differential equation is written D(x)y(x)=ly(x), where D(x) is a differential operator, l is the desired eigenvalue and y(x) the eigenfunction. A set of n constraints or boundary conditions must be given in order to completely determinethe solutions. Generally these take the form of values of the function y(x) and its derivatives at the initial and final points of an interval over which solutions are required. If all n boundary conditions are given for x=x1, then the problem becomes an initial value problem. But if n1 conditions are given at the initial point x1, y1(x1), ..., yn(x1) and that n2(=n-n1) are given at the final point x2, y1(x2),... then it is a "two point boundary value problem.
An n-th order problem can always be reduced to a system of n first-order problems. This is done by introducing functions y1(x) = dy/dx, ..., yn(x) = dny/dxn which satisfy coupled first order equations , i=1,2,...,n. Here the functions f'i are known.
If it is an initial value problem, with y1(x1) and y2(x1) given, then there are several standard numerical methods of solution such as Runge-Kutta methods and the Bulirsch-Stoer method. These methods employ small but finite steps dx, dy1, and dy2 in place of differentials dx, dy1, and dy2. in the formulas for dyi/dx to compute the new values of yi from their values at the preceding point.
If it is a two-point boundary value problem, then there are two classes of numerical methods to solve the problem.
The "shooting method" resembles the problem of hitting a target by adjusting the aim of an artillery piece. We imagine doing this by guessing (perhaps randomly) values for the functions yi at the initial point (those "free" n2 values that are unknown from initial conditions). We then integrate the initial value problem to arrive at the other boundary (the target). In general the solution values of yi at the final point will deviate from the required known n2 boundary values there. Our "score" in this shooting practice is determined by the differences between target boundary values and those obtained by integrating. Now we adjust the free initial point parameters to reduce the deviations at the final point. Iterations of this procedure are continued until the desired accuracy is achieved.
Relaxation methods use a different approach. These methods are preferred over shooting method when the solutions are smooth and not highly oscillatory. The authors of "Numerical Recipes" (self confessed notorious computer gunslingers) advise us to "always shoot first, and only then relax".
As is the case for the example introduced above, the boundary value problem may contain a parameter, denoted there by l. Then solutions may not be possible for arbitrary values of l. A slight modification of the procedure is then made using the following trick. Introduce a new dependent variable, yn+1 = l, and a differential equation dyn+1/dx = 0. The solution is again found by interation starting with a guess for free initial conditions. The original n-th order eigenvalue problem with boundary conditions has been converted to n+1 coupled first order DEs with boundary values.
At the initial point there are n+1 initial values to be specified but only n1 of them are known initial conditions. That leaves n2=n-n1 free starting values that have to be guessed. Denote these n2 values by the components of a vector v. Now, before iteration begins we must guess these component values; the other n1 initial conditions are known. The shooting method requires a function that, when given v at any step of the iteration will return all n+1=n1+n2 initial conditions as components of a vector y. Thus, define y¬load(x1,v) [initial values depend on v, y, and possibly x1].
Now, the coupled DEs are numerically integrated from x1 to x2 using the known functions f'i. At the final point, the function values constitute the components of a n+1 component vector y(x2). At x2 there are n2 known final boundary conditions: bk, k=1,2,...n2. Define F to be the n2 dimensional vector whose components are the difference of y(x2) from the known final boundary conditions, say Fk = bk - yk(x2). The components of F are the "score" of the shot from initial conditions v: F¬score(x2,y). A Newton-Raphson method based on numerical derivatives can be used to estimate an improved guess of initial conditions. These improved values are substituted for v and the procedure repeated. Repetitions are continued until some convergence criterion is met. E.g., all components of F might be required to lie within TOL of zero.
For an example of the application of this method, see the mathcad example for the particle in a box with central well.
W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T. Vetterling, Numerical Recipes [Cambridge University Press, 1986]
Created or up-dated 08/03/99
by R.D. Poshusta