Heath - Scientific Computing (523150), страница 95
Текст из файла (страница 95)
But for more complicatedPDEs, domains, boundary conditions, and discretization schemes, some of the methodslisted may not be viable options.• The methods listed are not equally reliable or robust. Many of the iterative methodsdepend on judicious choices of parameters or preconditioners that may be difficult todetermine in advance, and their performance may degrade significantly with choices that354CHAPTER 11.
PARTIAL DIFFERENTIAL EQUATIONSare less than optimal.• The methods listed are not equally easy to implement—some are relatively straightforward, but others involve complicated algorithms and data structures. Because of suchdifferences, the methods may also vary significantly in the relative speeds with which theyperform an equivalent amount of computation.• In practice the work may depend on implementation details. For example, the cost ofmultigrid is sensitive to the particular strategy used for cycling through grid levels.
Thefigures given in the table assume the best possible case.• Computational complexity alone does not necessarily determine the shortest time to asolution. For example, a method with high-order accuracy usually requires more workper grid point than a low-order method but may be able to use far fewer grid points toachieve equivalent accuracy.Table 11.3: Computational cost of solving elliptic boundary value problems as a function of k, thenumber of points along each dimension in a regular two- or three-dimensional gridMethodDense CholeskyJacobiGauss-SeidelBand CholeskyOptimal SORSparse CholeskyConjugate gradientOptimal SSORPreconditioned CGOptimal ADICyclic reductionFFTMultigrid V-cycleFACRFull Multigrid2-Dk6k 4 log kk 4 log kk43k log kk3k32.5k log kk 2.52k log2 kk 2 log kk 2 log kk 2 log kk 2 log log kk23-Dk9k 5 log kk 5 log kk74k log kk6k43.5k log kk 3.53k log2 kk 3 log kk 3 log kk 3 log kk 3 log log kk3In Table 11.3 the computational work is stated in terms of the number of grid pointsper dimension in a regular finite difference grid.
In Table 11.4 the equivalent value forthe work in terms of n, the order of the matrix, is given for those methods that remainviable choices for finite element discretizations with less regular meshes. The table gives theexponent of n in the dominant term of the cost estimate. These results are applicable tomany finite element discretizations and are also consistent with the estimates in Table 11.3for finite difference discretizations (with n = k 2 or n = k 3 , depending on the dimension ofthe problem).In Table 11.3 the methods are listed in decreasing order of cost for two-dimensionalproblems. Note that the ranking is somewhat different for three-dimensional problems.The reason is that factorization methods suffer a much greater penalty in going from twodimensions to three than do the other methods.
Although factorization methods are much11.7. SOFTWARE FOR PARTIAL DIFFERENTIAL EQUATIONS355Table 11.4: Exponent of n, the order of the matrix, in the computational cost of solving ellipticboundary value problems in two or three dimensionsMethodDense CholeskyBand CholeskySparse CholeskyConjugate gradientPreconditioned CGMultigrid2-D321.51.51.2513-D32.3321.331.171less competitive in terms of the work required for three-dimensional problems, they arestill useful in some cases because of their greater robustness. This is especially true fornonsymmetric matrices, since iterative methods tend to be significantly less reliable inthat case. In addition, methods akin to factorization are often used to compute effectivepreconditioners.The tables show that multigrid methods can be optimal, in the sense that the cost ofcomputing the solution is of the same order as the cost of reading the input or writingthe output.
The FACR method is also optimal for all practical purposes, since log log k iseffectively constant for any reasonable value of k. The other fast direct methods are almostas effective in practice unless k is very large. Clearly these methods should be seriouslyconsidered whenever they are applicable, and good software is available implementing them.Unfortunately, their robustness and applicability can be quite limited, so these optimal ornearly optimal methods, though they can be quite useful in the right context, are not apanacea, and more conventional methods must often be relied upon in practice.11.7Software for Partial Differential EquationsMost of the problem categories we have studied previously are amenable to reasonablyefficient solution by general-purpose software.
Methods for the numerical solution of partialdifferential equations, on the other hand, tend to be much more problem-dependent, so thatPDEs are most often solved using custom-written software to take maximum advantage ofthe particular features of a given problem. Nevertheless, some software does exist for afew general classes of problems that occur often in practice. Table 11.5 is a list of someof the software for the numerical solution of PDEs available from major libraries.
Generalproblem-solving environments for partial differential equations are also available, includinga PDE toolbox for MATLAB.Table 11.5: Library software for partial differential equationsSource 2-D Poisson 3-D Poisson Method of linesIMSLfps2hfps3hmolchNAGd03eafd03fafd03pcfNUMALrichardsonark/arkmatSLATEC hwscrt356CHAPTER 11. PARTIAL DIFFERENTIAL EQUATIONSIn addition, we list next several individual routines and software packages, many ofwhich are available from netlib, for solving various types of PDEs and also for solvingsparse linear systems.11.7.1Software for Initial Value Problems• CLAWPACK: 1-D and 2-D hyperbolic systems• TOMS: numerous routines and packages for solving various time-dependent problems inone or two space dimensions, including-11.7.2pdeone(#494): 1-D systems using method of linespdecol(#540): 1-D systems using collocationm3rk(#553): three-step Runge-Kutta method for parabolic equationspdetwo(#565): 2-D systems using method of linesbdmg(#621): 2-D nonlinear parabolic equationsepdcol(#688): 1-D systems using collocationpdecheb(#690): Chebyshev polynomial method for parabolic equationscwres(#731): moving-grid interface for 1-D systemsSoftware for Boundary Value Problems• ELLPACK: general framework and specific algorithms for solving various elliptic boundaryvalue problems on two- and three-dimensional domains [212]• FISHPACK: fast 2-D Helmholtz solvers using various coordinate systems• MGGHAT: second-order linear elliptic PDEs using finite element method with adaptive meshrefinement and a multigrid solver• PLTMG: elliptic PDEs with grid generation, adaptive mesh refinement, and a multigridsolver [13]• TOMS: numerous routines and packages for solving various elliptic boundary value problems(e.g., Poisson, Helmholtz) in two- or three-dimensional domains, including-11.7.3gma(#527): generalized marching algorithm for elliptic problemspwscrt(#541): fast 2-D Helmholtz solvers using various coordinate systemsfft9(#543): fast 2-D Helmholtz solverhelm3d(#572): fast 3-D Helmholtz solvercmmimp(#593): fast Helmholtz solver on nonrectangular planar regiongencol(#637): collocation method for general 2-D domainintcol/hermcol(#638): collocation method for rectangular domainhfft(#651): high-order fast 3-D Helmholtz solverserrg2/b2eval(#685): separable elliptic equations on rectangular domaincapc/reccn(#732): elliptic equations on irregular 2-D domainSoftware for Sparse Linear Systems• HSL: MA chapter of Harwell subroutine library contains numerous routines for solvingsparse linear systems11.8.
HISTORICAL NOTES AND FURTHER READING357• MATLAB: as of Version 4.0, sparse matrices are supported, including reorderings and factorizations• PCGPAK: preconditioned conjugate gradient methods for linear systems• QMRPACK: quasi-minimal residual methods for nonsymmetric linear systems• SLAP: iterative methods for symmetric and nonsymmetric linear systems• SPARSKIT: iterative methods for sparse linear systems and utilities for manipulating sparsematrices• SPARSPAK: reorderings and factorizations for sparse linear systems and least squares problems• SUPERLU: Gaussian elimination with partial pivoting for sparse linear systems• SYMMLQ: iterative method for symmetric indefinite linear systems• TEMPLATES: iterative methods documented in [15]• TOMS:- gpskca(#582): reordering sparse matrices for reduced bandwidth- lsqr(#583): iterative method for linear systems and least squares problems- itpack/nspcg(#586): stationary and nonstationary iterative methods for symmetric and nonsymmetric linear systems- sblas(#692): basic linear algebra subprograms for sparse matrices- jpicc/jpicr(#740): incomplete Cholesky factorization preconditioner• UMFPACK: unsymmetric multifrontal method for sparse linear systems• YSMP: Yale Sparse Matrix Package, direct methods for linear systems• Y12M: direct method for sparse linear systems11.8Historical Notes and Further ReadingThe literature on numerical solution of partial differential equations is vast.
Two earlyclassics on finite difference methods are [84, 213]. More recent treatments focusing mainlyon finite difference methods include [7, 112, 176, 186, 234, 249, 255]. Both finite differenceand finite element methods are discussed in [35, 117, 133, 159]. For a detailed discussionof the method of lines, see [219].
For an introduction to finite element methods from anumerical point of view, see [20, 50, 140, 265]. A deeper analysis of the finite elementmethod is given in [247]. A classic engineering text on finite elements is [281]; see also theseries [17, 33, 34]. Spectral and pseudospectral methods, though not treated in this book,form another important family of methods for solving differential equations; see [81, 106].Direct methods for solving sparse linear systems are discussed in detail in [68, 93]. Fastdirect methods, which were proposed by Hockney in 1965, are surveyed in [202, 253].Of the iterative methods we discussed, the Jacobi and Gauss-Seidel methods date fromthe nineteenth century, whereas the SOR and conjugate gradient methods were both developed around 1950 by Young and by Hestenes and Stiefel, respectively.