Time-Varying Partial Differential Equations
2335
The State of the Art in the Numerical Solution of Time-Varying Partial Differential Equations N. L. Schryer Bell Laboratories, Murray Hill, New Jersey 07974 (Received April 26, 1977) Publication costs assisted by Bell Laboratories
Theoretically, the state of the art is in excellent shape. A two-phase scheme consisting of spatial discretization followed by time discretization has been proven to converge for virtually every classical, textbook form of time-varying partial differential equation known. The time discretization is accomplished by “stifr differential equation solvers which are robust, efficient, and easy to use. In one space variable, the spatial discretization is well understood and there are several robust and efficient schemes available. In general, for problems with more than one space variable, the spatial discretization schemes are not well understood at all. The reasons for this are discussed and possible alternative spatial discretization schemes are outlined.
1. Introduction
This paper is intended to be a general overview of the state of the art in the numerical solution of time-varying partial differential equations (PDE’s),a field of numerical analysis which is important to chemists, physicists, engineers, biologists, and the physical sciences in general. The overview is aimed at two groups. The first consists of those people implementing numerical schemes for solving various problems. The second is the users of the software produced by the first group. Where useful, references are made to the relevant technical literature, but technical details will be avoided as much as possible. Rather, the advantages and disadvantages conferred by various technical facts will be discussed and motivated. Anyone interested in these technical details, including the idly curious, is encouraged to browse through the references. A rough outline of the areas to be discussed follows. T h e Theoretical State of the Art. Virtually A n y Classical Problem can be Solved. The numerical solution of a time-varying partial differential equation is a twophase process involving spatial discretization followed by time discretization. The spatial discretization is accomplished using either the finite-element (Galerkin-Rayleigh-Ritz) method, or finite differences (possibly with deferred corrections). The time evolution then requires the solution of a “stiff” system of ordinary differential equations (ODE’s). Such “stiff’ equations are typically solved by either backwards finite differences (with perhaps extrapolation), or “stiffly stable” predictor-corrector methods (Gear). This combination of spatial and time discretization techniques has been proven to converge for most classical, textbook problems, and it is clear that it applies equally well to other, more general, problems. Thus, the theoretical state of the art is in excellent shape. Since the numerical solution of “stiff’ ODE’s is efficiently accomplished, the practicality of the scheme depends upon the existence of efficient spatial discretization techniques. Practical State of the Art. Problems in One Space Variable. For problems in one space variable, there are robust and efficient spatial discretization techniques, and there are several software packages of very wide applicability and excellent quality for solving them. With this software, real and substantial problems can be solved in a short afternoon’s work. Due to the existence of this software, we have that the cost of formulating the
equations is much greater than the cost of solving them. Thus, a person who is neither a mathematician nor a numerical analyst can solve such problems personally and easily. Impractical State of the Art. Problems in More T h a n One Space Variable. In general, for problems in more than one space variable, robust and efficient spatial discretization techniques do not yet exist. However, simple problems, such as a single linear equation on a rectangular spatial domain, may be solved cleanly and fairly efficiently using existing software. Real world problems are, unfortunately, rarely that nice. Indeed, most real problems are systems of coupled PDE’s. Such general problems cannot yet be handled efficiently. In general, these problems can only be solved by armies of numerical-analyst programmers, well supported by lots of time and money. There are software packages for solving various restricted classes of problems in two and three spatial dimensions, both static and time varying. However, none of these packages agrees with any of the others on what a problem looks like. These packages are difficult to use and understand, with documentation being measured in cubic feet of paper. Thus, for problems in more than one spatial dimension, the state of the art is murky and clumsy. Due to this state of affairs, for this class of problems, we have that the cost of formulating the equations is much less than the cost of solving them. A novice cannot solve such problems without costly expert help from mathematicians and numerical analysts. Whither the State of the Art? In two or more space variables, several techniques have proven to be orders of magnitude cheaper than the current popular spatial discretization techniques for particularly simple, linear problems like the biharmonic equation and Laplace’s equation. These techniques have been generalized to other linear problems, and even some nonlinear problems, and they show great promise. However, they have not yet matured sufficiently, either theoretically or practically, to be of substantial help in general purpose PDE solving. This is currently a very active area of research. However, the mere creation of good numerical software is not enough. The user still has to be able to use that software. To enable simple usage of the numerical software packages for solving PDE’s, various “front-end’’ preprocessors are being developed. These preprocessors allow the user to speak in terms of the physical, rather than mathematical, notation with which they are familiar. The The Journal of Physlcal Chemlstty, Vol. 81, No. 25, 1977
2338
N. L. Schryer
user’s input is then automaticaily translated by the preprocessor into input suitable for the numerical software. Since the inputs for these preprocessors are computer languages, they are certain to benefit from recent developments in the theory, design, and implementation of compilers. The preprocessors will become easier to implement and use, with their linguistic features more clearly and simply defined and presented.
2. The Theoretical State of the Art Spatial discretization techniques are used to reduce a time-varying partial differential equation to a system of ordinary differential equations in time. This reduces a very difficult problem to one which is easily and efficiently solved. The solution u of the problem is approximated in space, x = (XI,x2, . . ,), by a sum of known basis functions B i ( x ) U ( t , x) =
N
Ci(t)Bi(X)
I= 1
where the coefficients ci(t) are time dependent, and the ci(0) are known. Such representations (discretizations) typically come from some grid or mesh which has been placed on the spatial domain. When (2.1) is inserted into the time-varying PDE, it results in a system of ODE’s in time for the coefficients ci(t). In general, these equations have the form
These equations are typically “stiff”. For example, when solving the grand-daddy of all time-varyingPDEs, the heat equation
aulat = v2u
(2.3)
using either the finite-difference (FD) method or the finite-element (FE) method, the stiffness of the resulting system of ODE’s is proportional to h-2,where h is the spatial mesh size used. Thus, as the accuracy request becomes tighter (as h decreases), the stiffness of the ODE system increases. It is quite common for these problems to have stiffness ratios on the order of lo9 to 10l2. Such stiff ODE’s may be easily and efficiently solved using any of several good stiff ODE solvers. Among the methods to choose from are the “stiffly stable” predictor-corrector methods of Gear,g the standard backward-difference techniques,20 extrapolated backwards-difference schemes,26,28 and exponentially-fitted schemes.12 These schemes are discussed in detail in ref 31. This technique, spatial discretization followed by solution of the resulting ODE’s in time, has been shown to converge for a number of classical problem^.^,^^^^^ It is routinely and successfully applied22i27 to problems which are not classical, and for which proofs of convergence have not yet been obtained. Thus, the theoretical state of the art is in excellent shape. The practicality of the scheme depends upon the existence of efficient spatial discretization techniques. 3. The Practical State of the Art The method of solution outlined in section 2 relies on an unspecified spatial discretization scheme. For problems in one space dimension, there are robust, reliable, accurate, and efficient spatial discretization schemes. Such problems can be solved efficiently. For problems in more than one space variable, such schemes do not currently exist, and The Journal of Physical Chemistry, Vol. 81, No. 25, 1977
these problems cannot, in general, be solved efficiently. This section discusses problems in one space variable and some of the software which has been built to solve them. There are two basic schemes for spatial discretization. The first and oldest is finite differences,8which is a very popular and well-understood process. A comparatively recent improvement to the basic FD scheme is the deferred-correction algorithm.17 As this scheme matures, the power of the FD method in general purpose PDE solving may be greatly enhanced. The second basic scheme is the finite element meth0d,2~9~~ which is becoming very popular. In this latter technique the spatial domain is divided into a finite number of smaller regions (elements). The solution is then approximated over each element by a polynomial of degree k - 1,where the order k is any number greater than 1 the user desires. A typical value of k is 4, and piecewise cubic polynomials are used. The polynomials are required to obey certain continuity conditions across the interfaces between the elements, so that the resulting approximation is smooth. If h is the mesh spacing used, then the error in the FE solution is O(hk).When this error is contrasted to the O(h2)error in the FD method, it is clear that the FE method is the more powerful and accurate, Moreover, to obtain the same error e, the mesh spacings used by FD’s and FEs, call them hfd and hf,, respective1 will obey O(h,d2) = e = O(hf:). Thus, roughly, hfe= hfd2 This means that far fewer mesh points are used by the FE method, when k is greater than 2. Note that for k > 2, since hf, > hfd, the stiffness of the ODE system generated by FE’s will be significantly lower than that generated by FD’s. Since stiffness in a problem makes it more difficult to solve, this makes the FE method even more appealing. It also shows that the stiffness of the ODE’s is an artifact of the spatial discretization techniques, and not inherent in solving time-varying PDE’s. Indeed, as even better spatial discretization schemes are developed, it is reasonable to expect that stiffness will nearly vanish in these problems. For problems in one space variable, these discretization schemes have proven quite adequate. If the problem is nice and smooth, and not much accuracy (say 1%)is desired, the FD method is appropriate. If the solution is not very smooth (e.g., if boundary layers or fronts are present), or if high accuracy is required, the FE method is much more efficient and preferable to FD’s. Because the spatial discretization of problems with only one space variable is in such good shape, the solution of time-varying PDE’s in one space variable may be considered trivial. Indeed, several general purpose software packages are available for solving such problems. These packages are listed as follows, in their order of ability to deal with difficult problems, from easiest to most difficult. (1) PDEPAKZ7 uses finite differences in both the time and spatial discretizations. Good for solving smooth problems to low accuracy. (2) PDECOLZ7uses a finite element scheme (collocation) in space and finite differences in time. Good for solving smooth problems to high accuracy. (3) POST^^ uses a finite element scheme (Galerkin’s method) in space and an extrapolated backwards difference scheme in time. Good for solving nonsmooth problems. Capable of dealing with coupled PDE-ODE systems, and nonlocal conditions, such as periodic boundary conditions. The rather small number of such general purpose software packages, and the nice hierarchy they form, are a good indication that the numerical solution of time-varying PDE’s in one space variable is well in hand. In fact, the calling sequences for these PDE solvers are
7,‘,
2337
Time-Varying Partial Differential Equations
only slightly more complex than those for any of the current good ODE solvers. The only added complications are due to the greater complexity of a PDE as compared to an ODE. Specifically, the spatial mesh and a subroutine for specifying the boundary conditions must also be supplied to the PDE software. The solution of problems in one space variable may be considered to be a mature field of numerical analysis. The existence of general purpose and reliable software means that users who are neither numerical analysts nor applied mathematicians can solve their problems without costly expert help.
4. The Impractical State of the Art For problems with two or more spatial variables, the situation is not nearly as nice as in the one-dimensional case. Given a single linear equation over a rectangular spatial domain, the situation is rather good. However, for problems involving systems of PDE’s or irregular spatial domains, the spatial discretization schemes are very poorly understood. Great strides have been made in recent years in the formulation and solution of the systems of linear algebraic equations which are the foundation of the FE and FD methods. For example, for a single linear PDE over a rectangular spatial domain, the system of linear equations can be formulated and solved in an “optimal” manner using current techn~logy.’~J~ The result is that for an n X n grid, O(n3)multiplications are used, while O(n2log (n)) storage locations are needed. However, for systems of PDE’s, or when the spatial domain is irregular, the problem of solving the system of linear equations becomes a much more expensive and complex business. The efficient solution of these equations requires that the equations and variables be reordered to minimize the number of arithmetic operations and storage locations used. This reordering has been compared [Tarjan, 3, pp 3-22] to the solution of several problems in computer science which are known to be intractable (i.e., NP complete). While the FE method is generally a great improvement over FD’s, it still requires vast amounts of computer time and memory for such problems. This is at least partly due to the complexity of the problem itself. Time-varying PDE’s in two and three spatial dimensions are fundamentally formidable beasts. However, just as with the stiffness of the ODE’S,the massive computer run-time and memory requirements of the FE method appear to be artifacts of the numerical solution process. For example, Laplace’s equation in two dimensions has recently been solved using techniques which are orders of magnitude cheaper and faster than the FE method.2 Similar developments have also taken place for the biharmonic equation5 These techniques do not discretize the spatial domain, but only its boundary. This reduces the number of spatial variables to be discretized by one, and results in a dramatic savings in memory and time. However, these newer techniques are currently only applicable to wellunderstood, linear problems such as the heat equation. Section 5 will consider the future of these methods. Another drawback to the FE procedure is the complexity of the computer code needed to implement it. In one space dimension, the code is rather straightforward. However, for more than one space dimension, the code is truly immense and complex. The coding and documentation times, even for experts in the area, are measured in man-years, not man-months. For all the difficulties which FE’s pose in two and more
spatial dimensions, they are still the best existing general purpose spatial discretization schemes available. They are also serviceable in the sense that if the problem is important enough, money can be found to create and run a FE code to solve it. Typically, however, only national governments and oil companies have problems whose economic importance justifies the expenditure (hence the adjective “Impractical” in the title of this section). Because such software is written, typically, with the solution of a specific problem in mind, and is of, by, and for the organization which creates it, most such software is horribly narrow-minded about its purpose in life. Many large companies, and governmental agencies, can boast of literally scores of software packages for solving various specific problems in two space variables. This may seem, on the face of it, to be a PDE solvers heaven. The user’s initial response is typically, “Golly, I found just exactly the right package for my problem!”. Indeed, for the first pass through a problem the user can usually find a package designed specifically for it-because someone was there before. However, the second pass through a problem invariably takes the user through untrodden territory. What is then encountered is the fact that facet no. 1 of the problem can be handled by software package A, while facet no. 2 (the new one) can be dealt with by package B. However, the converse is not the case. Thus, no existing software package can solve the entire second pass through the problem. One of two user responses inevitably follows this discovery: Either the problem is not solved or the user writes a software package which allows all facets of his current problem. Unfortunately, the next user of that new package will encounter precisely the same type of trouble. This process has indeed been the fecund mother of nearly all the currently extant software packages in the area. Due to the genealogy of the software (alteration piled upon alteration), its documentation is typically measured in cubic feet of paper. This thorough documentation, rather than being a guide to the user, is in fact a hydralike beast with which the user cannot deal. Due to this state of affairs, a learned person in a science other than mathematics will find it virtually impossible to do any significant modeling work using PDE’s in two or more space variables without costly and time-consuming help from numerical analysts, computer scientists, and applied mathematicians.
5. Whither the State of the Art? Several joint efforts, involving specialists in computer science, numerical analysis, applied mathematics, and the physical sciences, have produced general purpose software for dealing with a specific problem area. A good example of these efforts is the DISPL package15 for solving twodimensional kinetics-diffusion problems. This trend will undoubtedly accelerate, and general purpose software will become available within each discipline of the physical sciences. In the long term, generalizations of the ideas discussed briefly in section 3, for the very efficient solution of simple problems like Laplace’s equation and the biharmonic equation, appear to hold the most promise. In these methods, the spatial domain is not discretized, only the boundary of the domain is discretized. The techniques use basis functions B,(x),see (2.1), which are solutions of the PDE, but which do not satisfy the boundary conditions. The sum (2.1) is then forced to approximately satisfy the boundary conditions at each instant in time by choosing the coefficients ci(t) appropriately. This not only reduces the number of spatial variables to be discretized by one, The Journal of Physical Chemistry, Vol. 81, No. 25, 1977
N. L. Schryer
2338
but makes it possible to give sharp, rigorous bounds on the error in the approximate computed solution4s6via the maximum prin~ip1e.l~ Thus, not only are these techniques much faster and cheaper than FE's, for these problems, but they also allow the user to know precisely how accurate the computed solution is. Also, the solutions of PDE's in more than one space variable tend, in general, to have singular derivatives near corners in the boundary of the spatial domain.16 The ability of these new methods to analytically model and remove these singularities allows very kinky problems to be solved accurately and efficiently.21 These ideas have been generalized by Bergman' and Vekua30to linear, variable coefficient PDE's in two space variables. Others have recently generalized them for linear, variable coefficient PDE's in three spatial dimensions.6J1 These techniques have been shown to be practical and efficient for these problem^.^)^^ Moreover, for nonlinear static PDE's in two spatial dimensions, these ideas have been extended and shown to be very e f f i ~ i e n t When .~~~~~ (if) these techniques mature, the solution of PDE's in two and even three spatial dimensions will be several orders of magnitude cheaper than it currently is, both in terms of computer time and memory. Once this numerical nirvana has been attained, the user of the software will still have an obstacle to overcome. The basic software will be numerically oriented rather than application oriented. That is, the basic software will speak in terms of PDE's, yet the user will probably want (or only be able) to speak in terms of the physical properties of the problem at hand. Several preprocessors which convert such physical specifications into PDE's currently exist. Since the inputs for these processors are really languages (in the same sense that Fortran and Algol are), their design, implementation, and use will be made much easier and cleaner as tools developed in computer science for dealing with such languages are used. For example, lexical analyzer generatorsla and compiler-compilers14 have been developed which enable a language to be specified in terms of its lexicographical tokens and grammar, and which will then automatically create a parser for that language. Not only does this make the design and implementation of a new language a trivial task (compared with the previous by-hand approach), but it enables the designer to concentrate on the really important part of the project, namely, the language itself. This will allow more preprocessors to be built, and their linguistic features will be both cleaner and more useful. For example, a typical such "language" today has a specification which is a long list of rules for possible constructions, with an even longer list of exceptions to those rules. When modern software tools are used to design the language, the result is likely to be a simple and concise language specification which is both more powerful and easier to learn and use. The user community should expect this from the people designing their software interfaces, and even begin demanding it from them in t h e n e a r f u t u r e . References and Notes (1) S.Bergman, "Integral Operators and Partial Differential Equations", Ergebnisse der Mathematik und Ihrer Grenzgebiete, Neue Folge, Heft 23, 1961. (2) J. L. Blue, "Boundary Integral Solution of Laplace's Equation", submitted for publication. (3) J. R. Bunch and D. J. Rose, "Sparse Matrix Computations", Academic Press, New York, N.Y., 1976. (4) J. R. Cannon, SIAM J . Appl. Math., 12, 233-237 (1964). (5) J. R. Cannon and M. M. Cecchi, SIAM J . Numer. Anal., 3, 451-466 (1968). (6) D. Colton, Bull. Am. Math. SOC.,77, 752-756 (1971). (7) S.C. Eisenstat, SIAM J . Numer. Anal., 11, 654-680 (1974). The Journal of Physical Chemlstry, Vol. 81, No. 25, 1977
G. Forsythe and W. Wasow, "Finite Difference Methods for Partial Differential Equations", Wiley, New York, N.Y., 1959. C. W. Gear, "Numerical Initial Value Problems in Ordinary Differential Equations", Prentice-Hall, Englewood Cliffs, N.J., 197 1. J. A. George, SIAM J. Numer. Anal., 10, 345-363 (1973). R. P. Gilbert and D. Kukral, Bull. AMs, 79, 96-99 (1973). W.B. Gragg, "Lecture Notes on Extrapolation Methods", presented at the SIAM National Meeting, Washington, D.C., June 1971, and at the Conference on Ordinary Differential Equations, Dundee, Scotland, August, 1971. A. J. Hoffman, M. S.Martin, and D. J. Rose, SIAM J. Numer. Anal., 10, 364-369 (1973). S. C. Johnson, "YACC-Yet Another Compiler-Complier", Bell Laboratories Computing Science Technical Report No. 32, 1975. G. K. Leaf, M. Minkoff, G. D. Bryne, D. C. Sorenson, T. H. Bleakney, "DISPL: Computer Software for Two-Dimensional Kinetics-Diffusion Problems", presented at the SIAM Fall Meeting, Atlanta, Ga., Oct, 1976. R. S. Lehman, J . Math. Mech., 8, 729-760 (1959). M. Lentini and V. Pereyra, Math. Comput., 28, 981-1004 (1975). M. E. Lesk, "LEX-A Lexical Analyzer Generator", Bell Laboratories Computing Science Technical Report No. 39, 1975. M. H. Protter and H. F. Weinberger, "Maximum Principles in Differential Equations", Prentice-Hall, Englewood Cliffs, N.J., 1967. R. D. Richtmeyer and K. W. Morton, "Difference Methods for Initial Value Problems", Interscience, New York, N.Y., 1967. N. L. Schryer, SIAM J. Numer. Anal., 9, 546-572 (1972). N. L. Schryer, "Numerical Solution of Time-Varying Partial Differential Equations in One Space Variable", Bell Laboratories Computing Science Technical Report No. 53, 1976. N. L. Schryer, Numer. Math., 18, 336-344 (1972). N. L. Schryer, "A Tutorial on Galerkin's Method, using B-splines, for Solving Differential Equations", Bell Laboratories Computing Science Technical Report No. 52, 1976. N. L. Schryer, Numer. Math., 17, 284-300 (1971). N. L. Schryer and L. R. Walker, J. Appl. phys., 45,5406-5421 (1974). R. F. Sincovec and N. K. Madsen, Scientific Computing Consulting Services, P.O. Box 335, Manhattan, Kan. 66502. H. J. Stetter, Numer. Math., 7, 18-31 (1965). G. Strang and G. Fix, "An Analysis of the Finite Element Method", Prentice-Hall, New York, N.Y., 1973. I.N. Vekua, "New Methods for Solving Elliptic Equations", NorthHolland Publishing Co., Amsterdam, 1967. D. D. Warner, J. Phys. Chem., preceding paper in this issue.
Discussion W. A. REINHARDT (NASA, Ames Research Center). Y o u are adopting a n extremely pessimistic point-of-view o n our current abilities t o numerically solve multidimensional partial differential equations. D o y o u really i n t e n d t o also discount t h e recent advances made on solving the class o f PDE's which include the Eulerian equations (inviscid fluid flow equations), the boundary layer equations, the Poisson equation, the heat equation, as well as t h e wave equation (including weak solutions)?
N. L. SCHRYER. There have been m a n y impressive advances in the technology o f f i n i t e differences and finite elements. F o r certain simple, but fundamental, problems these advances have made numerical solution achievable and affordable on large scale computers (such as the IBM 370/168). However, for "real" problems (systems o f coupled nonlinear PDE's, irregular spatial domains) this technology can presently only be implemented and run o n supercomputers (such as t h e I l l i a c IV, ASC, a n d CDC 7600). T h i s means t h a t the numerical solution o f such problems i s beyond t h e reach o f the vast m a j o r i t y o f scientists a n d even organizations. Furthermore, recent advances in the study of the computational complexity o f these discretization techniques indicate that, in general, this situation w i l l n o t improve. Depending o n t h e reader's access t o supercomputers, this state of affairs m a y or m a y n o t lead t o pessimism.
W. A. REINHARDT (NASA, Ames Research Center). Although the cost of i n d i v i d u a l computers i s increasing, t h e r a p i d i t y of computer operations is increasing even more rapidly. This results in significantly lower costs o f a unit o f computation w i t h the supercomputers (Illiac IV, T I ASC, CDC Star, CDC 7600, Cray 1). Therefore, we have the apparent paradox t h a t although computations can be done cheaper, individual laboratories are even less l i k e l y t o make the required investments. Do y o u feel t h a t m a k i n g the supercomputers a n d methods accessible ( w i t h their lower unit o f cost) via a communication network will increase
Time-Varying Partial Differential Equations the viability of numerical methods for solving “real” problems?
N. L. SCHRYER. Making supercomputers available on a network will, on a case-by-case basis, permit solution of problems which were previously thought to be impractical to solve (using the available computing power). However, it is unreasonable to expect that computation will never become cheap enough to allow the general solution of P D E s in two and three spatial dimensions using finite-difference (FD) or finite-element (FE) methods. Recent theoretical results from complexity theory have cast a pall over the future of these methods as general purpose, useful tools in the PDE modeling process. In general, for problems in more than one space variable, the efficient solution of the systems of linear algebraic equations which are the result, and basis, of the FD and FE methods, requires the solution of an intermediate reordering problem. This reordering problem is known to be N P complete (ref 3), and is hence equivalent to the traveling salesman problem, and many other problems which have so far defied practical solution. In fact, the best known algorithms for solving any member of the class of NP-complete problems have run-times which can grow exponentially with the size of the problem. For problems in two space variables, if an N X N grid is used over an irregular spatial domain, this means that the cost of the reordering scheme can grow as e@. This is a worst-case cost, and specific problems may cost considerably less. However, this horrendous cost estimate raises several points to ponder: (1) There is undoubtedly, somewhere, a problem in two space variables whose reordering will actually cost that much. So even for a 10 x 10 grid that problem would be forever unsolvable.
2339 (2) There is currently no way to detect that a given reordering problem is a worst-case example. Thus, a human being may blunder ahead in a search for a nonexistent practical solution scheme for a specific problem, without ever knowing that it is a snipe hunt. (3) In those cases for which a practical reordering scheme exists, an intelligent human being, who is well versed in the combinatoric black arts, is required to find it. This then requires that the problem formulator (chemist) seek costly and cumbersome expert help from numerical analysts and mathematicians. With a human being in the middle of the numerical solution process, there are obvious reliability, cost, and time difficulties. While there are no theorems to prove it, the same type of pessimistic comments as above (for the reordering scheme) appear to apply to other schemes for solving the systems of linear algebraic equations. For example, the question of reordering the equations does not arise when they are to be solved iteratively (rather than directly). Yet the design of a practical iterative scheme for a given linear system is incredibly complex, with very small changes in the problem requiring a completely new and novel iteration scheme. Indeed, this wee fact of life is the source of the numerous and prosperous numerical consulting firms which do oil-well modeling for the petroleum industry. In conclusion, the efficient solution of the system of linear algebraic equations arising in the FD and F E techniques is, in general, so poorly understood as to be a black art rather than a generally useful tool. Furthermore, given that some of the approaches to this problem are N P complete, it is unlikely that the situation will ever improve sufficiently to convert the solution process from a black art to a useful science.
The Journal of Physical Chemistry, Vol. 81, No. 25, 1977
