Uri M. Ascher

Implicit-explicit methods for time-dependent partial differential equations

Implicit-explicit (IMEX) schemes have been widely used, especially in conjunction with spectral methods, for the time integration of spatially discretized partial differential equations (PDEs) of diffusion-convection type. Typically, an implicit scheme is used for the diffusion term and an explicit scheme is used for the convection term. Reaction-diffusion problems can also be approximated in this manner. In this work we systematically analyze the performance of such schemes, propose improved new schemes, and pay particular attention to their relative performance in the context of fast multigrid algorithms and of aliasing reduction for spectral methods.For the prototype linear advection-diffusion equation, a stability analysis for first-, second-, third-, and fourth-order multistep IMEX schemes is performed. Stable schemes permitting large time steps for a wide variety of problems and yielding appropriate decay of high frequency error modes are identified. Numerical experiments demonstrate that weak decay...

Collocation Software for Boundary-Value ODEs

The use of a general-purpose code, COLSYS, is described. The code is capable of solving mixed-order systems of boundary-value problems in ordinary differential equations. The method of spline collocation at Gaussian points is implemented using a B-spline basis. Approximate solutions are computed on a sequence of automatically selected meshes until a user-specified set of tolerances is satisfied. A damped Newtons method is used for the nonlinear iteration. The code has been found to be particularly effective for difficult problems. It is intended that a user be able to use COLSYS easily after reading its algorithm description. The use of the code is then illustrated by examples demonstrating its effectiveness and capabilities.

Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations

Abstract Implicit-explicit (IMEX) linear multistep time-discretization schemes for partial differential equations have proved useful in many applications. However, they tend to have undesirable time-step restrictions when applied to convection-diffusion problems, unless diffusion strongly dominates and an appropriate BDF-based scheme is selected (Ascher et al., 1995). In this paper, we develop Runge-Kutta-based IMEX schemes that have better stability regions than the best known IMEX multistep schemes over a wide parameter range.

A collocation solver for mixed order systems of boundary value problems

Implementation of a spline collocation method for solving boundary value problems for mixed order systems of ordinary differential equations is discussed. The aspects of this method considered include error estimation, adaptive mesh selection, B-spline basis function evaluation, linear system solution and nonlinear problem solution. The resulting general purpose code, COLSYS, is tested on a number of examples to demonstrate its stability, efficiency and flexibility.

On optimization techniques for solving nonlinear inverse problems

This paper considers optimization techniques for the solution of nonlinear inverse problems where the forward problems, like those encountered in electromagnetics, are modelled by differential equations. Such problems are often solved by utilizing a Gauss-Newton method in which the forward model constraints are implicitly incorporated. Variants of Newtons method which use second-derivative information are rarely employed because their perceived disadvantage in computational cost per step offsets their potential benefits of faster convergence. In this paper we show that, by formulating the inversion as a constrained or unconstrained optimization problem, and by employing sparse matrix techniques, we can carry out variants of sequential quadratic programming and the full Newton iteration with only a modest additional cost. By working with the differential equation explicitly we are able to relate the constrained and the unconstrained formulations and discuss the advantages of each. To make the comparisons meaningful we adopt the same global optimization strategy for all inversions. As an illustration, we focus upon a 1D electromagnetic (EM) example simulating a magnetotelluric survey. This problem is sufficiently rich that it illuminates most of the computational complexities that are prevalent in multi-source inverse problems and we therefore describe its solution process in detail. The numerical results illustrate that variants of Newtons method which utilize second-derivative information can produce a solution in fewer iterations and, in some cases where the data contain significant noise, requiring fewer floating point operations than Gauss-Newton techniques. Although further research is required, we believe that the variants proposed here will have a significant impact on developing practical solutions to large-scale 3D EM inverse problems.

A new basis implementation for a mixed order boundary value ODE solver

The numerical approximation of mixed order systems of multipoint value ordinary differential equations by collocation requires appropriate representation of the piecewise polynomial solutions. B-splines were originally implemented in the general purpose code COLSYS, but better alternatives exist. One promising alternative as proposed by Osborne and discussed by Ascher, Pruess and Russell. In this paper we analyze the performance of the latter solution representation for cases not previously covered, where the mesh is not necessarily dense. This analysis and other considerations have led us to implement a basis replacement in COLSYS and we discuss some implementation details. Numerical results are given which demonstrate the improvement in performance of the code.

Projected implicit Runge-Kutta methods for differential-algebraic equations

In this paper a new class of numerical methods, Projected Implicit Runge–Kutta methods, is introduced for the solution of index-2 Hessenberg systems of initial and boundary value differential-algebraic equations (DAEs). These types of systems arise in a variety of applications, including the modeling of singular optimal control problems and parameter estimation for differential-algebraic equations such as multibody systems. The new methods appear to be particularly promising for the solution of DAE boundary value problems, where the need to maintain stability in the differential part of the system often necessitates the use of methods based on symmetric discretizations. Previously defined symmetric methods have severe limitations when applied to these problems, including instability, oscillation, and loss of accuracy; the new methods overcome these difficulties. For linear problems we define an essential underlying boundary value ODE and prove well-conditioning of the differential (or state-space) solutio...

Stabilization of Constrained Mechanical Systems with DAEs and Invariant Manifolds

ABSTRACT Many methods have been proposed for the simulation of constrained mechanical systems. The most obvious of these have mild instabilities and drift problems. Consequently, stabilization techniques have been proposed A popular stabilization method is Baumgartes technique, but the choice of parameters to make it robust has been unclear in practice. Some of the simulation methods that have been proposed and used in computations are reviewed here, from a stability point of view. This involves concepts of differential-algebraic equation (DAE) and ordinary differential equation (ODE) invariants. An explanation of the difficulties that may be encountered using Baumgartes method is given, and a discussion of why a further quest for better parameter values for this method will always remain frustrating is presented. It is then shown how Baumgartes method can be improved. An efficient stabilization technique is proposed, which may employ explicit ODE solvers in case of nonstiff or highly oscillatory probl...

Inversion of 3D electromagnetic data in frequency and time domain using an inexact all-at-once approach

We present a general formulation for inverting frequencyor time-domain electromagnetic data using an all-at-once approach. In this methodology, the forward modeling equations are incorporated as constraints and, thus, we need to solve a constrained optimization problem where the parameters are the electromagnetic fields, the conductivity model, and a set of Lagrange multipliers. This leads to a much larger problem than the traditional unconstrained formulation where only the conductivities are sought. Nevertheless, experience shows that the constrained problem can be solved faster than the unconstrained one. The primary reasons are that the forward problem does not have to be solved exactly until the very end of the optimization process, and that permitting the fields to be away from their constrained values in the initial stages introduces flexibility so that a stationary point of the objective function is found more quickly. In this paper, we outline the all-atonce approach and apply it to electromagnetic problems in both frequency and time domains. This is facilitated by a unified representation for forward modeling for these two types of data. The optimization problem is solved by finding a stationary point of the Lagrangian. Numerically, this leads to a nonlinear system that is solved iteratively using a Gauss-Newton strategy. At each iteration, a large, indefinite matrix is inverted, and we discuss how this can be accomplished. As a test, we invert frequency-domain synthetic data from a grounded electrode system that emulates a field CSAMT survey. For the time domain, we invert borehole data obtained from a current loop on the surface.

Implicit-Explicit Methods for Time-Dependent PDE''s

Implicit-explicit (IMEX) schemes have been widely used, especially in conjunction with spectral methods, for the time integration of spatially discretized PDEs of diffusion-convection type. Typically, an implicit scheme is used for the diffusion term and an explicit scheme is used for the convection term. Reaction-diffusion problems can also be approximated in this manner. In this work we systematically analyze the performance of such schemes, propose improved new schemes and pay particular attention to their relative performance in the context of fast multigrid algorithms and of aliasing reduction for spectral methods. For the prototype linear advection-diffusion equation, a stability analysis for first, second, third and fourth order multistep IMEX schemes is performed. Stable schemes permitting large time steps for a wide variety of problems and yielding appropriate decay of high frequency error modes are identified. Numerical experiments demonstrate that weak decay of high frequency modes can lead to extra iterations on the finest grid when using multigrid computations with finite difference spatial discretization, and to aliasing when using spectral collocation for spatial discretization. When this behaviour occurs, use of weakly damping schemes such as the popular combination of Crank-Nicolson with second order Adams-Bashforth is discouraged and better alternatives are proposed. Our findings are demonstrated on several examples.