Numerical analysis of partial differential equations using Maple and MATLAB.

*(English)*Zbl 1459.65001
Fundamentals of Algorithms 12. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM) (ISBN 978-1-61197-530-7/pbk; 978-1-61197-531-4/ebook). ix, 153 p. (2018).

The book is a part of the SIAM series on “Fundamentals and Algorithms”. The user-oriented books on state-of-the-art numerical methods are written with the aim to assist the reader to choose an appropriate method for his application, to understand the method and its limitations, and to implement the algorithm by runnable codes. It is based on one-semester courses given independently by the authors at the McGill University and at the University of Geneva between 2001 and 2012 for students from mathematics, engineering, and computer science.

The subject is focused on elliptic partial differential equations. The 153 + IX pages long monograph consists of an introduction, four chapters for the main discretization methods – the finite difference method, the finite volume method, the spectral method, and the finite element method –, bibliography, and an index. The necessary numerical analysis with complete convergence proofs for each method is treated, and codes in MATLAB for each algorithm are presented. The chapters are largely self-contained and can be studied independently, also without taking a course. The text is completed by historical remarks and quotes from the main contributors including figures. The chapters about the discretization methods are finished by concluding remarks, exercises, and projects for the reader. The authors underline that they fill a gap to provide a unified representation of the main discretization methods. Distinguished books for the several discretization techniques are well-known.

The introduction starts with notations how the scientific computing tools Maple and MATLAB handle derivatives, partial derivatives, gradients, divergence, etc. in order to prepare the user for the codes. The solution of ordinary differential equations with Maple and MATLAB is demonstrated for the pendulum problem based on Newton’s second law of motion using an adaptive fourth-order Runge-Kutta method, the Lotka-Volterra model for the predator-prey interaction, and the forecasting of crude oil production based on economic constraints. The authors continue with a classification of the partial differential equations. They shortly derive the corresponding grand challenge PDEs: heat equation, advection-reaction-diffusion equation, wave equation, Maxwell’s equations, and Navier-Stokes equations including the equations of Stokes and Euler as special cases.

The introduction is finished by the main topic of the monograph, elliptic equations. The authors formulate three classes of elliptic problems. Time-dependent problems can be solved by semidiscretization using a finite difference method, e.g., by a forward or backward Euler method, a Crank-Nicolson method etc. The spatial part of the equation can be elliptic and is solved in every time step. Elliptic problems also arise as equilibrium solutions. It is demonstrated by the steady-state solutions of the heat equation, the wave equation, the Maxwell-Faraday law in electrodynamics, the advection-diffusion equation, and the Navier-Stokes equations. In this context the Poisson and the Laplace equation are established. Time-harmonic regimes, in which the forcing term and the boundary conditions are periodic with a frequency, result in a third case of elliptic equations. Corresponding approaches for the heat equation and the second-order wave equation result in the shifted Laplace equation or the Helmholtz equation, respectively.

Chapter 2, “The finite difference method”, starts with historical remarks. The method was introduced by C. Runge [Schlömilch Z. 56, 225–232 (1908; JFM 39.0433.01)]. The first convergence proof was given by R. Amsler [Bull. Soc. Math. Fr. 56, 141–166 (1928; JFM 54.0485.01)] and the first error estimate by S. A. Gershgorin [Z. Angew. Math. Mech. 10, 373–382 (1930; JFM 56.0467.03)].

The authors start with the finite difference method for the two-dimensional Poisson equation with Dirichlet boundary conditions on the unit square, introducing the forward, the backward, and the centered approximation using a truncated Taylor series to approximate the derivatives in each variable. The corresponding linear system of equations for the second-order accurate centered approximation is derived. For the numerical solution of linear systems of equations the authors refer to other publications. A convergence analysis of the mentioned centered finite difference approximation of the Poisson equation with Dirichlet boundary conditions is established based on the truncation error estimate, the discrete maximum principle, and a Poincare-type estimate.

The authors continue with more accurate approximations and more general boundary conditions, such as Neumann and Robin boundary conditions with the use of ghost points outside the domain for the discretization. For nonrectangular domains the standard five-point difference discretization of Laplacians (centered approximation) is adapted, if a uniform grid leads to nodes which are not located on the boundary. Neumann and Robin boundary conditions for nonrectangular domains are connected with the difficulty, that the normal derivative generally does not have the same direction as the grid lines. The authors refer for this case to more appropriate algorithms given in Chapters 3 and 5.

Several generalizations of the finite difference approximations are established on the basis of general differential operators. The finite difference approximation of the three-dimensional Poisson equation is formulated. Using a stationary, nonlinear advection-reaction-diffusion equation, which contains also a first-order derivative, three choices of discretization for the first-order derivative are discussed with the result, that one get an approximate solution with non-physical solutions, if the step size is not small enough, or a solution which is always physically correct (upwind scheme), or a result, that is incorrect (downwind scheme).

Chapter 3, “The finite volume method”, starts with the remark, that the above mentioned difficulties of the finite difference approximation for nonrectangular domains have been the main reason to invent a new technique. The now so-called finite volume method works for arbitrary meshes and geometries, and was introduced by P. W. McDonald [“The computation of transonic flow through two-dimensional gas turbine cascades”, in: ASME TurboExpo: Power for Land, Sea, and Air. ASME 1971 International Gas Turbine Conference and Products Show. American Society of Mechanical Engineers. Paper 71-GT-89, 7 p. (1971; doi:10.1115/71-GT-89)]. The method consists roughly spoken of three steps. The equation is first integrated over a small control volume. Then the volume integral is transformed into a boundary integral using the divergence theorem involving fluxes. Finally the fluxes are approximated across the boundaries. The method allows one to handle more general meshes with complicated boundary conditions. For rectangular meshes the finite volume method leads to the same difference stencils as the corresponding difference approximation.

For the convergence analysis of the finite volume method new ways are required in comparison to the finite difference method. A corresponding theorem is formulated and proved. Quadratic convergence is often observed, but it has not yet been generally proven.

The subject of Chapter 4 is the spectral method. The idea is to develop the solution of the PDE as a sum of certain basis functions, such as trigonometric functions (Fourier series) or orthogonal polynomials, and then to determine the coefficients in the sum to satisfy the PDE as well as possible. The spectral method is a global one in contrast to the methods discussed above. Ritz was the first (1908) to use this form as a computational method for vibrating plates and is therefore considered as the father of the spectral method.

The authors start with the approximation of the one-dimensional Poisson equation with periodic boundary conditions inserting the Fourier series representation into the equation and determine the unknown Fourier coefficients using the orthogonality of the exponential function. Based on the decay of the Fourier coefficients the truncation error of the infinite expansion is studied and the so-called spectral convergence of the method is derived. The spectral convergence is compared with the convergence obtained for the finite difference and finite volume methods. The spectral method with discrete Fourier series including the fast Fourier transform (FFT) and the corresponding exponential convergence is described in detail.

For nonperiodic problems Chebyshev polynomials are used. The grid points of these orthogonal polynomials are not equidistant in contrast to interpolating polynomials with the disadvantage of excessive oscillations. The Chebyshev spectral method is extensively described. The numerical solutions of the Chebyshev and the finite difference method for a number of degrees of freedom are graphically illustrated for the problem \( u_{xx} = f \) and show the exponential convergence of the Chebyshev method as a function of the regularity of the right-hand side \(f\).

The approximations by the spectral methods are for smooth solutions much more accurate than the one of the finite difference and finite volume methods. But, they can only be used for simple rectangular geometries. However, the authors refer at the end of this chapter already to Chapter 5 for a combination of the spectral method with the possibility to decompose arbitrary geometries into rectangular regions resulting in the spectral finite element method.

In Chapter 5, “The finite element method”, the authors start with the assessment that the FEM is the “most flexible method for the numerical solution of PDEs” and the one “with the most solid mathematical foundation”. The main idea goes back to W. Ritz [J. Reine Angew. Math. 135, 1–61 (1908; JFM 39.0449.01)], who started directly from the principle of least action, from which the variational formulation of the equations and the boundary conditions result. It is minimized over a finite-dimensional subspace rather than over all functions of the infinite-dimensional space. The approximation of the global minimum increases with the dimension of the approximation and with the choice of appropriate functions.

Instead of using the minimization formulation of W. Ritz, B. G. Galerkin [“Rods and plates. Series occurring in various questions concerning the elastic equilibrium of rods and plates” (Russian), Engineers Bull. (Vestnik Inzhenerov) 19, 897–908 (1915); reprinted in: Selected works, Vol. I. Moskva: Izdatel’stvo “Nauka”. 168–195 (1952); see Zbl 0048.42105 for the entire volume] proposed to use the same finite-dimensional approach but to work directly with the differential equation resulting later in the Galerkin method of weighted residuals. More information about the history of the finite element method including the contribution of Courant one can find in the introduction of this chapter.

The Poisson equation in one-dimensional spatial dimension is used to introduce the finite element method. Based on the variational formulation of the equation one would get the so-called strong form. But, in order to develop the finite element discretization, the weak or variational form is used, established by multiplying the Poisson equation by any function \( v \in V \), by integrating over the domain, applying integration by parts, and by formulating and solving the minimization problem. Corresponding to the approach mentioned above, a finite-dimensional subspace \[ V_h := \text{span} \{{\phi}_1, \dots ,{\phi}_n\} \subset V \] is used for the discretization, and the Galerkin approximation is formulated including the resulting linear system of equations with the associated stiffness matrix.

Rather to compute the integrals caused by the right-hand side \(f\) of the Poisson equation one can also approximate \(f\) by a linear combination of the functions \({\phi}_i\) resulting in a system of linear equation with the associated so-called mass matrix. It is proved that stiffness and mass matrix are symmetric and positive definite. The method is demonstrated by an example with piecewise linear hat functions \({\phi}_i\). Deriving the Ritz approximation shows the equivalence with the Galerkin method. Discussed is also the inclusion of Dirichlet and Neumann boundary conditions in the FEM.

Sobolev spaces are introduced in order to analyze and generalize the FEM. The convergence analysis is based on Céa’s lemma, on an interpolation estimate, and on the Aubin-Nitsche lemma. The analysis is presented for the one-dimensional model problem and extended to a two-dimensional problem. Algorithms for mesh generation, mesh refinement, and mesh smoothing are implemented for the two-dimensional case (triangulation) using MATLAB. Finally, finite element shape functions are introduced as a restriction of the hat functions to each finite element. Finite element shape functions allow for computing the stiffness matrix and the integrals element by element. The element stiffness matrix is calculated on each element and than added at the appropriate location in the global stiffness matrix. This so-called assembly process is similarly applied for the global mass matrix.

The authors refer after all to a Wikipedia page that contains a huge list of FEM software packages. The FEM chapter is related to the node finite elements. An extension to the edge finite elements would be helpful for the readers of this well-written book.

The subject is focused on elliptic partial differential equations. The 153 + IX pages long monograph consists of an introduction, four chapters for the main discretization methods – the finite difference method, the finite volume method, the spectral method, and the finite element method –, bibliography, and an index. The necessary numerical analysis with complete convergence proofs for each method is treated, and codes in MATLAB for each algorithm are presented. The chapters are largely self-contained and can be studied independently, also without taking a course. The text is completed by historical remarks and quotes from the main contributors including figures. The chapters about the discretization methods are finished by concluding remarks, exercises, and projects for the reader. The authors underline that they fill a gap to provide a unified representation of the main discretization methods. Distinguished books for the several discretization techniques are well-known.

The introduction starts with notations how the scientific computing tools Maple and MATLAB handle derivatives, partial derivatives, gradients, divergence, etc. in order to prepare the user for the codes. The solution of ordinary differential equations with Maple and MATLAB is demonstrated for the pendulum problem based on Newton’s second law of motion using an adaptive fourth-order Runge-Kutta method, the Lotka-Volterra model for the predator-prey interaction, and the forecasting of crude oil production based on economic constraints. The authors continue with a classification of the partial differential equations. They shortly derive the corresponding grand challenge PDEs: heat equation, advection-reaction-diffusion equation, wave equation, Maxwell’s equations, and Navier-Stokes equations including the equations of Stokes and Euler as special cases.

The introduction is finished by the main topic of the monograph, elliptic equations. The authors formulate three classes of elliptic problems. Time-dependent problems can be solved by semidiscretization using a finite difference method, e.g., by a forward or backward Euler method, a Crank-Nicolson method etc. The spatial part of the equation can be elliptic and is solved in every time step. Elliptic problems also arise as equilibrium solutions. It is demonstrated by the steady-state solutions of the heat equation, the wave equation, the Maxwell-Faraday law in electrodynamics, the advection-diffusion equation, and the Navier-Stokes equations. In this context the Poisson and the Laplace equation are established. Time-harmonic regimes, in which the forcing term and the boundary conditions are periodic with a frequency, result in a third case of elliptic equations. Corresponding approaches for the heat equation and the second-order wave equation result in the shifted Laplace equation or the Helmholtz equation, respectively.

Chapter 2, “The finite difference method”, starts with historical remarks. The method was introduced by C. Runge [Schlömilch Z. 56, 225–232 (1908; JFM 39.0433.01)]. The first convergence proof was given by R. Amsler [Bull. Soc. Math. Fr. 56, 141–166 (1928; JFM 54.0485.01)] and the first error estimate by S. A. Gershgorin [Z. Angew. Math. Mech. 10, 373–382 (1930; JFM 56.0467.03)].

The authors start with the finite difference method for the two-dimensional Poisson equation with Dirichlet boundary conditions on the unit square, introducing the forward, the backward, and the centered approximation using a truncated Taylor series to approximate the derivatives in each variable. The corresponding linear system of equations for the second-order accurate centered approximation is derived. For the numerical solution of linear systems of equations the authors refer to other publications. A convergence analysis of the mentioned centered finite difference approximation of the Poisson equation with Dirichlet boundary conditions is established based on the truncation error estimate, the discrete maximum principle, and a Poincare-type estimate.

The authors continue with more accurate approximations and more general boundary conditions, such as Neumann and Robin boundary conditions with the use of ghost points outside the domain for the discretization. For nonrectangular domains the standard five-point difference discretization of Laplacians (centered approximation) is adapted, if a uniform grid leads to nodes which are not located on the boundary. Neumann and Robin boundary conditions for nonrectangular domains are connected with the difficulty, that the normal derivative generally does not have the same direction as the grid lines. The authors refer for this case to more appropriate algorithms given in Chapters 3 and 5.

Several generalizations of the finite difference approximations are established on the basis of general differential operators. The finite difference approximation of the three-dimensional Poisson equation is formulated. Using a stationary, nonlinear advection-reaction-diffusion equation, which contains also a first-order derivative, three choices of discretization for the first-order derivative are discussed with the result, that one get an approximate solution with non-physical solutions, if the step size is not small enough, or a solution which is always physically correct (upwind scheme), or a result, that is incorrect (downwind scheme).

Chapter 3, “The finite volume method”, starts with the remark, that the above mentioned difficulties of the finite difference approximation for nonrectangular domains have been the main reason to invent a new technique. The now so-called finite volume method works for arbitrary meshes and geometries, and was introduced by P. W. McDonald [“The computation of transonic flow through two-dimensional gas turbine cascades”, in: ASME TurboExpo: Power for Land, Sea, and Air. ASME 1971 International Gas Turbine Conference and Products Show. American Society of Mechanical Engineers. Paper 71-GT-89, 7 p. (1971; doi:10.1115/71-GT-89)]. The method consists roughly spoken of three steps. The equation is first integrated over a small control volume. Then the volume integral is transformed into a boundary integral using the divergence theorem involving fluxes. Finally the fluxes are approximated across the boundaries. The method allows one to handle more general meshes with complicated boundary conditions. For rectangular meshes the finite volume method leads to the same difference stencils as the corresponding difference approximation.

For the convergence analysis of the finite volume method new ways are required in comparison to the finite difference method. A corresponding theorem is formulated and proved. Quadratic convergence is often observed, but it has not yet been generally proven.

The subject of Chapter 4 is the spectral method. The idea is to develop the solution of the PDE as a sum of certain basis functions, such as trigonometric functions (Fourier series) or orthogonal polynomials, and then to determine the coefficients in the sum to satisfy the PDE as well as possible. The spectral method is a global one in contrast to the methods discussed above. Ritz was the first (1908) to use this form as a computational method for vibrating plates and is therefore considered as the father of the spectral method.

The authors start with the approximation of the one-dimensional Poisson equation with periodic boundary conditions inserting the Fourier series representation into the equation and determine the unknown Fourier coefficients using the orthogonality of the exponential function. Based on the decay of the Fourier coefficients the truncation error of the infinite expansion is studied and the so-called spectral convergence of the method is derived. The spectral convergence is compared with the convergence obtained for the finite difference and finite volume methods. The spectral method with discrete Fourier series including the fast Fourier transform (FFT) and the corresponding exponential convergence is described in detail.

For nonperiodic problems Chebyshev polynomials are used. The grid points of these orthogonal polynomials are not equidistant in contrast to interpolating polynomials with the disadvantage of excessive oscillations. The Chebyshev spectral method is extensively described. The numerical solutions of the Chebyshev and the finite difference method for a number of degrees of freedom are graphically illustrated for the problem \( u_{xx} = f \) and show the exponential convergence of the Chebyshev method as a function of the regularity of the right-hand side \(f\).

The approximations by the spectral methods are for smooth solutions much more accurate than the one of the finite difference and finite volume methods. But, they can only be used for simple rectangular geometries. However, the authors refer at the end of this chapter already to Chapter 5 for a combination of the spectral method with the possibility to decompose arbitrary geometries into rectangular regions resulting in the spectral finite element method.

In Chapter 5, “The finite element method”, the authors start with the assessment that the FEM is the “most flexible method for the numerical solution of PDEs” and the one “with the most solid mathematical foundation”. The main idea goes back to W. Ritz [J. Reine Angew. Math. 135, 1–61 (1908; JFM 39.0449.01)], who started directly from the principle of least action, from which the variational formulation of the equations and the boundary conditions result. It is minimized over a finite-dimensional subspace rather than over all functions of the infinite-dimensional space. The approximation of the global minimum increases with the dimension of the approximation and with the choice of appropriate functions.

Instead of using the minimization formulation of W. Ritz, B. G. Galerkin [“Rods and plates. Series occurring in various questions concerning the elastic equilibrium of rods and plates” (Russian), Engineers Bull. (Vestnik Inzhenerov) 19, 897–908 (1915); reprinted in: Selected works, Vol. I. Moskva: Izdatel’stvo “Nauka”. 168–195 (1952); see Zbl 0048.42105 for the entire volume] proposed to use the same finite-dimensional approach but to work directly with the differential equation resulting later in the Galerkin method of weighted residuals. More information about the history of the finite element method including the contribution of Courant one can find in the introduction of this chapter.

The Poisson equation in one-dimensional spatial dimension is used to introduce the finite element method. Based on the variational formulation of the equation one would get the so-called strong form. But, in order to develop the finite element discretization, the weak or variational form is used, established by multiplying the Poisson equation by any function \( v \in V \), by integrating over the domain, applying integration by parts, and by formulating and solving the minimization problem. Corresponding to the approach mentioned above, a finite-dimensional subspace \[ V_h := \text{span} \{{\phi}_1, \dots ,{\phi}_n\} \subset V \] is used for the discretization, and the Galerkin approximation is formulated including the resulting linear system of equations with the associated stiffness matrix.

Rather to compute the integrals caused by the right-hand side \(f\) of the Poisson equation one can also approximate \(f\) by a linear combination of the functions \({\phi}_i\) resulting in a system of linear equation with the associated so-called mass matrix. It is proved that stiffness and mass matrix are symmetric and positive definite. The method is demonstrated by an example with piecewise linear hat functions \({\phi}_i\). Deriving the Ritz approximation shows the equivalence with the Galerkin method. Discussed is also the inclusion of Dirichlet and Neumann boundary conditions in the FEM.

Sobolev spaces are introduced in order to analyze and generalize the FEM. The convergence analysis is based on Céa’s lemma, on an interpolation estimate, and on the Aubin-Nitsche lemma. The analysis is presented for the one-dimensional model problem and extended to a two-dimensional problem. Algorithms for mesh generation, mesh refinement, and mesh smoothing are implemented for the two-dimensional case (triangulation) using MATLAB. Finally, finite element shape functions are introduced as a restriction of the hat functions to each finite element. Finite element shape functions allow for computing the stiffness matrix and the integrals element by element. The element stiffness matrix is calculated on each element and than added at the appropriate location in the global stiffness matrix. This so-called assembly process is similarly applied for the global mass matrix.

The authors refer after all to a Wikipedia page that contains a huge list of FEM software packages. The FEM chapter is related to the node finite elements. An extension to the edge finite elements would be helpful for the readers of this well-written book.

Reviewer: Georg Hebermehl (Berlin)

##### MSC:

65-01 | Introductory exposition (textbooks, tutorial papers, etc.) pertaining to numerical analysis |

65N06 | Finite difference methods for boundary value problems involving PDEs |

65N08 | Finite volume methods for boundary value problems involving PDEs |

65N35 | Spectral, collocation and related methods for boundary value problems involving PDEs |

65N30 | Finite element, Rayleigh-Ritz and Galerkin methods for boundary value problems involving PDEs |

35Jxx | Elliptic equations and elliptic systems |

65N12 | Stability and convergence of numerical methods for boundary value problems involving PDEs |

65-04 | Software, source code, etc. for problems pertaining to numerical analysis |