Showing papers on "Numerical analysis published in 2005"

Fourth-Order Time-Stepping for Stiff PDEs

Abstract: A modification of the exponential time-differencing fourth-order Runge--Kutta method for solving stiff nonlinear PDEs is presented that solves the problem of numerical instability in the scheme as proposed by Cox and Matthews and generalizes the method to nondiagonal operators. A comparison is made of the performance of this modified exponential time-differencing (ETD) scheme against the competing methods of implicit-explicit differencing, integrating factors, time-splitting, and Fornberg and Driscoll's "sliders" for the KdV, Kuramoto--Sivashinsky, Burgers, and Allen--Cahn equations in one space dimension. Implementation of the method is illustrated by short MATLAB programs for two of the equations. It is found that for these applications with fixed time steps, the modified ETD scheme is the best.

Applied Numerical Methods Using MATLAB

Abstract: Preface. 1. MATLAB Usage and Computational Errors. 1.1 Basic Operations of MATLAB. 1.1.1 Input/Output of Data from MATLAB Command Window. 1.1.2 Input/Output of Data Through Files. 1.1.3 Input/Output of Data Using Keyboard. 1.1.4 2-D Graphic Input/Output. 1.1.5 3-D Graphic Output. 1.1.6 Mathematical Functions. 1.1.7 Operations on Vectors and Matrices. 1.1.8 Random Number Generators. 1.1.9 Flow Control. 1.2 Computer Errors Versus Human Mistakes. 1.2.1 IEEE 64-bit Floating-Point Number Representation. 1.2.2 Various Kinds of Computing Errors. 1.2.3 Absolute/Relative Computing Errors. 1.2.4 Error Propagation. 1.2.5 Tips for Avoiding Large Errors. 1.3 Toward Good Program. 1.3.1 Nested Computing for Computational Efficiency. 1.3.2 Vector Operation Versus Loop Iteration. 1.3.3 Iterative Routine Versus Nested Routine. 1.3.4 To Avoid Runtime Error. 1.3.5 Parameter Sharing via Global Variables. 1.3.6 Parameter Passing Through Varargin. 1.3.7 Adaptive Input Argument List. Problems. 2. System of Linear Equations. 2.1 Solution for a System of Linear Equations. 2.1.1 The Nonsingular Case (M = N). 2.1.2 The Underdetermined Case (M N): Least-Squares Error Solution. 2.1.4 RLSE (Recursive Least-Squares Estimation). 2.2 Solving a System of Linear Equations. 2.2.1 Gauss Elimination. 2.2.2 Partial Pivoting. 2.2.3 Gauss-Jordan Elimination. 2.3 Inverse Matrix. 2.4 Decomposition (Factorization). 2.4.1 LU Decomposition (Factorization): Triangularization. 2.4.2 Other Decomposition (Factorization): Cholesky, QR, and SVD. 2.5 Iterative Methods to Solve Equations. 2.5.1 Jacobi Iteration. 2.5.2 Gauss-Seidel Iteration. 2.5.3 The Convergence of Jacobi and Gauss-Seidel Iterations. Problems. 3. Interpolation and Curve Fitting. 3.1 Interpolation by Lagrange Polynomial. 3.2 Interpolation by Newton Polynomial. 3.3 Approximation by Chebyshev Polynomial. 3.4 Pade Approximation by Rational Function. 3.5 Interpolation by Cubic Spline. 3.6 Hermite Interpolating Polynomial. 3.7 Two-dimensional Interpolation. 3.8 Curve Fitting. 3.8.1 Straight Line Fit: A Polynomial Function of First Degree. 3.8.2 Polynomial Curve Fit: A Polynomial Function of Higher Degree. 3.8.3 Exponential Curve Fit and Other Functions. 3.9 Fourier Transform. 3.9.1 FFT Versus DFT. 3.9.2 Physical Meaning of DFT. 3.9.3 Interpolation by Using DFS. Problems. 4. Nonlinear Equations. 4.1 Iterative Method Toward Fixed Point. 4.2 Bisection Method. 4.3 False Position or Regula Falsi Method. 4.4 Newton(-Raphson) Method. 4.5 Secant Method. 4.6 Newton Method for a System of Nonlinear Equations. 4.7 Symbolic Solution for Equations. 4.8 A Real-World Problem. Problems. 5. Numerical Differentiation/Integration. 5.1 Difference Approximation for First Derivative. 5.2 Approximation Error of First Derivative. 5.3 Difference Approximation for Second and Higher Derivative. 5.4 Interpolating Polynomial and Numerical Differential. 5.5 Numerical Integration and Quadrature. 5.6 Trapezoidal Method and Simpson Method. 5.7 Recursive Rule and Romberg Integration. 5.8 Adaptive Quadrature. 5.9 Gauss Quadrature. 5.9.1 Gauss-Legendre Integration. 5.9.2 Gauss-Hermite Integration. 5.9.3 Gauss-Laguerre Integration. 5.9.4 Gauss-Chebyshev Integration. 5.10 Double Integral. Problems. 6. Ordinary Differential Equations. 6.1 Euler's Method. 6.2 Heun's Method: Trapezoidal Method. 6.3 Runge-Kutta Method. 6.4 Predictor-Corrector Method. 6.4.1 Adams-Bashforth-Moulton Method. 6.4.2 Hamming Method. 6.4.3 Comparison of Methods. 6.5 Vector Differential Equations. 6.5.1 State Equation. 6.5.2 Discretization of LTI State Equation. 6.5.3 High-Order Differential Equation to State Equation. 6.5.4 Stiff Equation. 6.6 Boundary Value Problem (BVP). 6.6.1 Shooting Method. 6.6.2 Finite Difference Method. Problems. 7. Optimization. 7.1 Unconstrained Optimization [L-2, Chapter 7]. 7.1.1 Golden Search Method. 7.1.2 Quadratic Approximation Method. 7.1.3 Nelder-Mead Method [W-8]. 7.1.4 Steepest Descent Method. 7.1.5 Newton Method. 7.1.6 Conjugate Gradient Method. 7.1.7 Simulated Annealing Method [W-7]. 7.1.8 Genetic Algorithm [W-7]. 7.2 Constrained Optimization [L-2, Chapter 10]. 7.2.1 Lagrange Multiplier Method. 7.2.2 Penalty Function Method. 7.3 MATLAB Built-In Routines for Optimization. 7.3.1 Unconstrained Optimization. 7.3.2 Constrained Optimization. 7.3.3 Linear Programming (LP). Problems. 8. Matrices and Eigenvalues. 8.1 Eigenvalues and Eigenvectors. 8.2 Similarity Transformation and Diagonalization. 8.3 Power Method. 8.3.1 Scaled Power Method. 8.3.2 Inverse Power Method. 8.3.3 Shifted Inverse Power Method. 8.4 Jacobi Method. 8.5 Physical Meaning of Eigenvalues/Eigenvectors. 8.6 Eigenvalue Equations. Problems. 9. Partial Differential Equations. 9.1 Elliptic PDE. 9.2 Parabolic PDE. 9.2.1 The Explicit Forward Euler Method. 9.2.2 The Implicit Backward Euler Method. 9.2.3 The Crank-Nicholson Method. 9.2.4 Two-Dimensional Parabolic PDE. 9.3 Hyperbolic PDE. 9.3.1 The Explicit Central Difference Method. 9.3.2 Two-Dimensional Hyperbolic PDE. 9.4 Finite Element Method (FEM) for solving PDE. 9.5 GUI of MATLAB for Solving PDEs: PDETOOL. 9.5.1 Basic PDEs Solvable by PDETOOL. 9.5.2 The Usage of PDETOOL. 9.5.3 Examples of Using PDETOOL to Solve PDEs. Problems. Appendix A: Mean Value Theorem. Appendix B: Matrix Operations/Properties. Appendix C: Differentiation with Respect to a Vector. Appendix D: Laplace Transform. Appendix E: Fourier Transform. Appendix F: Useful Formulas. Appendix G: Symbolic Computation. Appendix H: Sparse Matrices. Appendix I: MATLAB. References. Subject Index. Index for MATLAB Routines. Index for Tables.

Algorithms for Numerical Analysis in High Dimensions

Abstract: Nearly every numerical analysis algorithm has computational complexity that scales exponentially in the underlying physical dimension. The separated representation, introduced previously, allows many operations to be performed with scaling that is formally linear in the dimension. In this paper we further develop this representation by (i) discussing the variety of mechanisms that allow it to be surprisingly efficient; (ii) addressing the issue of conditioning; (iii) presenting algorithms for solving linear systems within this framework; and (iv) demonstrating methods for dealing with antisymmetric functions, as arise in the multiparticle Schrodinger equation in quantum mechanics. Numerical examples are given.

Analysis of three‐dimensional crack initiation and propagation using the extended finite element method

Abstract: An Erratum has been published for this article in International Journal for Numerical Methods in Engineering 2005, 63(8): 1228. We present a new formulation and a numerical procedure for the quasi-static analysis of three-dimensional crack propagation in brittle and quasi-brittle solids. The extended finite element method (XFEM) is combined with linear tetrahedral elements. A viscosity-regularized continuum damage constitutive model is used and coupled with the XFEM formulation resulting in a regularized ‘crack-band’ version of XFEM. The evolving discontinuity surface is discretized through a C0 surface formed by the union of the triangles and quadrilaterals that separate each cracked element in two. The element's properties allow a closed form integration and a particularly efficient implementation allowing large-scale 3D problems to be studied. Several examples of crack propagation are shown, illustrating the good results that can be achieved. Copyright © 2005 John Wiley & Sons, Ltd.

On generalized successive overrelaxation methods for augmented linear systems

Abstract: For the augmented system of linear equations, Golub, Wu and Yuan recently studied an SOR-like method (BIT 41(2001)71---85). By further accelerating it with another parameter, in this paper we present a generalized SOR (GSOR) method for the augmented linear system. We prove its convergence under suitable restrictions on the iteration parameters, and determine its optimal iteration parameters and the corresponding optimal convergence factor. Theoretical analyses show that the GSOR method has faster asymptotic convergence rate than the SOR-like method. Also numerical results show that the GSOR method is more effective than the SOR-like method when they are applied to solve the augmented linear system. This GSOR method is further generalized to obtain a framework of the relaxed splitting iterative methods for solving both symmetric and nonsymmetric augmented linear systems by using the techniques of vector extrapolation, matrix relaxation and inexact iteration. Besides, we also demonstrate a complete version about the convergence theory of the SOR-like method.

Comparison of finite difference and boundary integral solutions to three-dimensional spontaneous rupture

Abstract: The spontaneously propagating shear crack on a frictional interface has proven to be a useful idealization of a natural earthquake. The corresponding boundary value problems are nonlinear and usually require computationally intensive numerical methods for their solution. Assessing the convergence and accuracy of the numerical methods is challenging, as we lack appropriate analytical solutions for comparison. As a complement to other methods of assessment, we compare solutions obtained by two independent numerical methods, a finite difference method and a boundary integral (BI) method. The finite difference implementation, called DFM, uses a traction-at-split-node formulation of the fault discontinuity. The BI implementation employs spectral representation of the stress transfer functional. The three-dimensional (3-D) test problem involves spontaneous rupture spreading on a planar interface governed by linear slip-weakening friction that essentially defines a cohesive law. To get a priori understanding of the spatial resolution that would be required in this and similar problems, we review and combine some simple estimates of the cohesive zone sizes which correspond quite well to the sizes observed in simulations. We have assessed agreement between the methods in terms of the RMS differences in rupture time, final slip, and peak slip rate and related these to median and minimum measures of the cohesive zone resolution observed in the numerical solutions. The BI and DFM methods give virtually indistinguishable solutions to the 3-D spontaneous rupture test problem when their grid spacing Δx is small enough so that the solutions adequately resolve the cohesive zone, with at least three points for BI and at least five node points for DFM. Furthermore, grid-dependent differences in the results, for each of the two methods taken separately, decay as a power law in Δx, with the same convergence rate for each method, the calculations apparently converging to a common, grid interval invariant solution. This result provides strong evidence for the accuracy of both methods. In addition, the specific solution presented here, by virtue of being demonstrably grid-independent and consistent between two very different numerical methods, may prove useful for testing new numerical methods for spontaneous rupture problems.

Convergence and stability of a discontinuous galerkin time-domain method for the 3D heterogeneous maxwell equations on unstructured meshes

Abstract: A Discontinuous Galerkin method is used for to the numerical solution of the time-domain Maxwell equations on unstructured meshes. The method relies on the choice of local basis functions, a centered mean approximation for the surface integrals and a second-order leap-frog scheme for advancing in time. The method is proved to be stable for cases with either metallic or absorbing boundary conditions, for a large class of basis functions. A discrete analog of the electromagnetic energy is conserved for metallic cavities. Convergence is proved for Discontinuous elements on tetrahedral meshes, as well as a discrete divergence preservation property. Promising numerical examples with low-order elements show the potential of the method.

Parallel Lagrange--Newton--Krylov--Schur Methods for PDE-Constrained Optimization. Part I: The Krylov--Schur Solver

Abstract: Large-scale optimization of systems governed by partial differential equations (PDEs) is a frontier problem in scientific computation. Reduced quasi-Newton sequential quadratic programming (SQP) methods are state-of-the-art approaches for such problems. These methods take full advantage of existing PDE solver technology and parallelize well. However, their algorithmic scalability is questionable; for certain problem classes they can be very slow to converge. In this two-part article we propose a new method for steady-state PDE-constrained optimization, based on the idea of using a full space Newton solver combined with an approximate reduced space quasi-Newton SQP preconditioner. The basic components of the method are Newton solution of the first-order optimality conditions that characterize stationarity of the Lagrangian function; Krylov solution of the Karush--Kuhn--Tucker (KKT) linear systems arising at each Newton iteration using a symmetric quasi-minimum residual method; preconditioning of the KKT system using an approximate state/decision variable decomposition that replaces the forward PDE Jacobians by their own preconditioners, and the decision space Schur complement (the reduced Hessian) by a BFGS approximation initialized by a two-step stationary method. Accordingly, we term the new method {\it Lagrange--Newton--Krylov--Schur} (LNKS). It is fully parallelizable, exploits the structure of available parallel algorithms for the PDE forward problem, and is locally quadratically convergent. In part I of this two-part article, we investigate the effectiveness of the KKT linear system solver. We test our method on two optimal control problems in which the state constraints are described by the steady-state Stokes equations. The objective is to minimize dissipation or the deviation from a given velocity field; the control variables are the boundary velocities. Numerical experiments on up to 256 Cray T3E processors and on an SGI Origin 2000 include scalability and performance assessment of the LNKS algorithm and comparisons with reduced SQP for up to $1,000,000$ state and 50,000 decision variables. In part II of the article, we address globalization and inexactness issues, and apply LNKS to the optimal control of the steady incompressible Navier--Stokes equations.

A New Finite Element Gradient Recovery Method: Superconvergence Property

Abstract: This is the first in a series of papers in which a new gradient recovery method is introduced and analyzed. It is proved that the method is superconvergent for translation invariant finite element spaces of any order. The method maintains the simplicity, efficiency, and superconvergence properties of the Zienkiewicz--Zhu patch recovery method. In addition, for uniform triangular meshes, the method is superconvergent for the linear element under the chevron pattern, and ultraconvergent at element edge centers for the quadratic element under the regular pattern. Applications of this new gradient recovery technique will be discussed in forthcoming papers.

A finite volume method for the laplace equation on almost arbitrary two-dimensional grids

Abstract: We present a finite volume method based on the integration of the Laplace equation on both the cells of a primal almost arbitrary two-dimensional mesh and those of a dual mesh obtained by joining the centers of the cells of the primal mesh. The key ingredient is the definition of discrete gradient and divergence operators verifying a discrete Green formula. This method generalizes an existing finite volume method that requires "Voronoi-type" meshes. We show the equivalence of this finite volume method with a non-conforming finite element method with basis functions being P 1 on the cells, generally called "diamond-cells", of a third mesh. Under geometrical conditions on these diamond- cells, we prove a first-order convergence both in the H 1 norm and in the L 2 norm. Superconvergence results are obtained on certain types of homothetically refined grids. Finally, numerical experiments confirm these results and also show second-order convergence in the L 2 norm on general grids. They also indicate that this method performs particularly well for the approximation of the gradient of the solution, and may be used on degenerating triangular grids. An example of application on non- conforming locally refined grids is given.

The Adomian decomposition method for solving delay differential equation

Abstract: A numerical method based on the Adomian decomposition method which has been developed by Adomian [Adomian, G. (1994). Solving Frontier Problems of Physics: The Decomposition Method. Kluwer Academic Publishers, Boston, MA] is introduced in this paper for the approximate solution of delay differential equations. The algorithm is illustrated by studying an initial value problem. The results obtained are presented and show that only a few terms are required to obtain an approximate solution which is found to be accurate and efficient.

Numerical methods for nonlinear stochastic differential equations with jumps

Abstract: We present and analyse two implicit methods for Ito stochastic differential equations (SDEs) with Poisson-driven jumps. The first method, SSBE, is a split-step extension of the backward Euler method. The second method, CSSBE, arises from the introduction of a compensated, martingale, form of the Poisson process. We show that both methods are amenable to rigorous analysis when a one-sided Lipschitz condition, rather than a more restrictive global Lipschitz condition, holds for the drift. Our analysis covers strong convergence and nonlinear stability. We prove that both methods give strong convergence when the drift coefficient is one-sided Lipschitz and the diffusion and jump coefficients are globally Lipschitz. On the way to proving these results, we show that a compensated form of the Euler–Maruyama method converges strongly when the SDE coefficients satisfy a local Lipschitz condition and the pth moment of the exact and numerical solution are bounded for some p>2. Under our assumptions, both SSBE and CSSBE give well-defined, unique solutions for sufficiently small stepsizes, and SSBE has the advantage that the restriction is independent of the jump intensity. We also study the ability of the methods to reproduce exponential mean-square stability in the case where the drift has a negative one-sided Lipschitz constant. This work extends the deterministic nonlinear stability theory in numerical analysis. We find that SSBE preserves stability under a stepsize constraint that is independent of the initial data. CSSBE satisfies an even stronger condition, and gives a generalization of B-stability. Finally, we specialize to a linear test problem and show that CSSBE has a natural extension of deterministic A-stability. The difference in stability properties of the SSBE and CSSBE methods emphasizes that the addition of a jump term has a significant effect that cannot be deduced directly from the non-jump literature.