    Next: 6.4.3 Shooting Method Up: 6.4 Vertical Discretization of Previous: 6.4.1 Two-Point Boundary Value

## 6.4.2 Numerical Solution Techniques for Boundary Value Problems

The numerical solution of BVPs for ODEs is a long studied and well-understood subject in numerical mathematics. Many textbooks have been published. Two books treating the topic at different depths shall vicariously be cited: In [132, chs. 16 and 17] an introduction to ``Numerical Recipes'' for initial value problems as well as for boundary value problems is given. Basic considerations are presented and simple numerical algorithms are derived. This book is especially recommendable for getting into the field, a comprehensive treatment, however, cannot be found there due to the brevity of the discussion. The classical textbook for advanced numerical methods is the monograph  by Uri M. Ascher et al.. This book provides a complete treatment of the topic including advanced recent developments. Starting with the basic theory of ODEs, the two distinct classes of numerical methods for solving two point BVPs--namely, initial value methods and finite-difference methods--are derived from the very beginning to very sophisticated numerical techniques like decoupling, automatic mesh selection and error estimation, as well as perturbation theory. We shall now give a brief survey of the two available classes of finite-difference and initial value methods. For details see .

Initial value methods. Boundary value problems were long considered some sort of a subclass of initial value problems (IVPs) whereby special care has to be taken of the initial conditions in order to get things right at the other end. This opinion has changed and IVPs are considered to be actually a special and in some sense relatively simple subclass of BVPs. The fundamental difference is that for IVPs one has the complete information about the solution at one point, i.e., the initial point, so some marching algorithm which is always local in nature may be considered. For BVPs, on the other hand, no complete information is available at any point, so the end points have to be connected by the solution algorithm in a global way. Only after stepping through the entire domain the solution at any point can be determined. In spite of the fundamental difference of the problem nature methods suited for initial values can also be applied to BVPs. The simplest initial value method is the shooting method. Values for all of the missing dependent variables at one boundary are chosen. These values must, of course, be consistent with any BCs for that boundary, but otherwise they are arranged to depend on arbitrary free parameters whose values are initially ``randomly'' chosen. Then the ODE system is integrated by any suitable initial value method, arriving at the other boundary. In general, there will be discrepancies from the desired boundary values there. However, by adjustment of the free parameters at the starting point the discrepancy can be zeroed. The BVP is thus transferred into a multidimensional root-finding problem. In case of linear ODEs a linear algebraic system has to be solved--a situation that applies to our problem. Unfortunately, this simple algorithm suffers from potential numerical instabilities. More advanced methods like reduced superposition, multiple shooting, reorthogonalization, and decoupling moderate these instabilities at the cost of increased numerical demands. The main advantage of initial value methods is the possibility of a very memory-saving implementation besides their conceptual simplicity since they somehow exploit the local nature of IVPs. In our case the number of unknowns is extremely high, the vertical variation of 4000 and more Fourier coefficients has to be determined. Hence, we implemented a very advanced IVP that follows the stabilized march algorithm derived in [200, pp. 155-164]. This algorithm has also successfully been implemented in the general purpose code MUSL  available through NETLIB.d Our implementation is tailored to the specific structure of the ODE (6.32) to enhance the performance; the core of the implementation, however, was taken from the MUSL code.

Finite-difference methods. Finite-difference or relaxation methods use a different approach. The ODEs are replaced by finite-difference equations on a mesh of points that covers the entire range of the integration interval. A trial solution of the values of the dependent variables at each mesh point is first guessed and thus generally does not satisfy the desired finite-difference equations, nor does necessarily satisfy the required BCs. Subsequent iterations also called relaxations try to adjust all the values on the mesh so as to bring them into successively closer agreement with the finite-difference equations and, simultaneously, with the BCs. Obviously the success of the relaxation depends on the quality of the guess, i.e., the closer it is to the true solution the better the performance. The number of values that have to be adjusted thus equals the number of unknowns multiplied by the number of mesh points. In our case this corresponds to roughly 4 105 values, i.e., 4000 unknowns and 100 mesh points. This number is far too big to be solved on today's workstations. Hence relaxation techniques are not suited for our problem. The dramatically high numerical demands result from the global solution strategy. In some cases--provided that the size of the problem is not too large--relaxation methods are the only means to solve the BVP since they work better than initial value methods when the BCs are especially subtle, or where they involve complicated algebraic relations that cannot easily be solved in closed form. In our case the BCs are simple linear relations and thus very easily evaluated. Furthermore, relaxation works better when the solution is smooth and not oscillatory. Such oscillations would require many grid points for accurate representation. Additionally the number and position of required points is not known a priori. Here lies another strength of initial value methods since they can employ variable step-size integrators to adjust naturally to the solution peculiarity. However, relaxation methods are advantageous in case of extraneous solutions which, while not appearing in the final solution satisfying all BCs, may wreak havoc on the initial value integrations required by shooting. A typical case is that of trying to maintain a dying exponential simultaneously with growing exponentials. The ODE system is then said to be stiff. Unfortunately, exactly this situation occurs in our problem since the waves traveling downwards the resist are damped whereas the waves traveling upwards the resist seem to be growing if one looks in the same direction both times, e.g., into the resist. Hence no simple shooting can be used in our case, and even the advanced techniques used can moderate such problems only up to a certain amount. The limitations of our algorithm thus caused are described in Section 6.5. Relaxation methods would overcome them but are prohibitively computation- and memory-intensive.

#### Footnotes

... NETLIB.d
The acronym MUSL stands for multiple shooting for linear boundary value problems. The code was written by Robert Mattheij and G. Staarink and is available at http://www.netlib.org/ode/mus*.*.    Next: 6.4.3 Shooting Method Up: 6.4 Vertical Discretization of Previous: 6.4.1 Two-Point Boundary Value
Heinrich Kirchauer, Institute for Microelectronics, TU Vienna
1998-04-17