Backward Error Analysis as a Model of Computation for Numerical Methods

This thesis delineates a generally applicable perspective on numerical meth­ ods for scientific computation called residual-based a posteriori backward er­ ror analysis, based on the concepts of condition, backward error, and residual, pioneered by Turing and Wilkinson. The basic underpinning of this perspec­ tive, that a numerical method’s errors should be analyzable in the same terms as physical and modelling errors, is readily understandable across scientific fields, and it thereby provides a view of mathematical tractability readily in­ terpretable in the broader context of mathematical modelling. It is applied in this thesis mainly to numerical solution of differential equations. We examine the condition of initial-value problems for ODEs and present a residual-based error control strategy for methods such as Euler’s method, Taylor series meth­ ods, and Runge-Kutta methods. We also briefly discuss solutions of continuous chaotic problems and stiff problems

Discontinuities in numerical radiative transfer

Observations and magnetohydrodynamic simulations of solar and stellar atmospheres reveal an intermittent behavior or steep gradients in physical parameters, such as magnetic field, temperature, and bulk velocities. The numerical solution of the stationary radiative transfer equation is particularly challenging in such situations, because standard numerical methods may perform very inefficiently in the absence of local smoothness. However, a rigorous investigation of the numerical treatment of the radiative transfer equation in discontinuous media is still lacking. The aim of this work is to expose the limitations of standard convergence analyses for this problem and to identify the relevant issues. Moreover, specific numerical tests are performed. These show that discontinuities in the atmospheric physical parameters effectively induce first-order discontinuities in the radiative transfer equation, reducing the accuracy of the solution and thwarting high-order convergence. In addition, a survey of the existing numerical schemes for discontinuous ordinary differential systems and interpolation techniques for discontinuous discrete data is given, evaluating their applicability to the radiative transfer problem

Efficient goal-oriented global error estimation for BDF-type methods using discrete adjoints

This thesis develops estimation techniques for the global error that occurs during the approximation of solutions of Initial Value Problems (IVPs) on given intervals by multistep integration methods based on Backward Differentiation Formulas (BDF). To this end, discrete adjoints obtained by adjoint Internal Numerical Differentiation (IND) of the nominal integration scheme are used. For this purpose, a bridge between BDF methods and Petrov-Galerkin Finite Element (FE) methods is built by a novel functional-analytic framework. Goal-oriented global error estimators are derived in analogy to the Dual Weighted Residual methodology in Galerkin methods for Partial Differential Equations. Their asymptotic behavior, their accuracy in BDF methods with variable order and stepsize as well as their applicability for global error control are investigated. The novel results presented in this thesis include: i) a functional-analytic framework for IVPs in Ordinary Differential Equations (ODEs) in the Banach space of continuously differentiable functions. This framework is needed since the classical Hilbert space setting is not suitable to analyze the relation between the discrete values of the adjoint IND scheme and the solution of the adjoint IVP. The new framework gives rise to the definition of weak solutions of adjoint IVPs. ii) a Petrov-Galerkin FE discretization of the function spaces that allows to transform the variational formulations of the IVP and of its adjoint IVP into finite dimensional problems. The equivalence of these finite dimensional problems to BDF methods with variable but prescribed order and stepsize and their adjoint IND schemes is shown. Thus, the FE approximation of the weak adjoint is determined by the discrete values of the adjoint IND scheme and discretization and differentiation commute in the developed framework. iii) a proof that the values of the adjoint IND scheme corresponding to a BDF method with constant order and stepsize converge to the solution of the adjoint IVP on the open interval. In addition, a proof is given that demonstrates the convergence of the FE approximation to the weak solution of the adjoint IVP on the entire interval. iv) goal-oriented global error estimators for BDF methods that weight, for each integration step, a local error quantity with the corresponding value of the adjoint IND scheme and yields in sum an accurate and efficient estimate for the actual error. As local error quantity defect integrals and local truncation errors are employed, respectively. v) strategies for goal-oriented global error control in BDF methods that either adapt the locally acting relative tolerance or the given integration scheme using the stepwise error indicators. vi) an ODE model of an exothermic, self-accelerating chemical reaction with mass transfer carried out in a discontinuous Stirred Tank Reactor. With this real-world example from chemical engineering the applicability and reliability of the novel techniques for the approximation of weak adjoints and for the simulation with goal-oriented global error control are shown

SciCADE 95: International conference on scientific computation and differential equations

This report consists of abstracts from the conference. Topics include algorithms, computer codes, and numerical solutions for differential equations. Linear and nonlinear as well as boundary-value and initial-value problems are covered. Various applications of these problems are also included

Global error evaluation strategies in multistep methods applied to ordinary differential equations and index 1 differential algebraic systems

Differential Equations (DEs) are among the most widely used mathematical tools in different area of sciences. Solving DEs, either analytically or numerically, has become a centre of interest for many mathematicians and a large variety of methods are nowadays available to solve DEs numerically. When solving a mathematical problem numerically, evaluating the error is of high importance in practice. Most of the methods already available for solving DEs are implemented with a mechanism to perform a local error control. However, in the real realm, it is common to require the numerical solution to approximate the exact solution with accuracy to a certain number of decimal places or significant figures. To satisfy this condition, we require the global error to be bounded by a specifically determined tolerance. In this case, a local error control is not longer efficient. On one hand, controlling the local error only cannot ensure that the required accuracy will be achieved. On the other hand, the use of such approach requires the user to do some preliminary studies on the problem, and have deep understanding of the method. Thus, we need a mechanism to control the global error in order to compute the numerical solution for a user-supplied accuracy requirement in automatic mode. The global error estimate calculated in the course of such a control can also be applied to improve the numerical solution obtained. It is straight forward since, if the error estimate is found with sufficiently high accuracy, we can just add it to the numerical solution to get a better approximation to the exact value. Thus, accurate evaluation of the the global error is crucial for the purpose mentioned above. Several techniques are already developed to compute the global error of the numerical solution. The most common algorithms include the Richardson extrapolation, Zadunaisky’s technique, Solving for the correction, and Using two different methods. These methods use two integrations to evaluate the global error, and the provided error estimate is valid if the global error admits an expansion in powers of the step size. Another approach, known as solving the linearised discrete variational equation, can also be used. This last differs from the others by the use of a truncated Taylor expansion of the defect of the method to estimate the global error; and solving the problem and estimating the error is roughly the same as one step of the underlying method. In this research, we will investigate numerically and compare the efficiency of different techniques for global error evaluation applied to multistep methods for solving ordinary differential equations (ODEs) and differential algebraic equations (DAEs). We will first study the global error evaluation techniques in multistep formulas for solving ODEs on uniform grids. In the case of nonuniform grids, both multistep methods with variable coefficients and interpolation-type multistep methods will be considered. Then, we will extend our study to multistep methods for solving DAEs. Theoretical background will accompany numerical works. The accuracy and reliability of the global error evaluation strategies will be discussed and compared for different types of multistep methods for solving ODEs and DAEs. We will analyse the efficiency in terms of accuracy obtained and CPU time spent. For that, a series of numerical experiments is conducted on a set of test problems with known solutions

Systematic construction of efficient six-stage fifth-order explicit Runge-Kutta embedded pairs without standard simplifying assumptions

This thesis examines methodologies and software to construct explicit Runge-Kutta (ERK) pairs for solving initial value problems (IVPs) by constructing efficient six-stage fifth-order ERK pairs without standard simplifying assumptions. The problem of whether efficient higher-order ERK pairs can be constructed algebraically without the standard simplifying assumptions dates back to at least the 1960s, with Cassity's complete solution of the six-stage fifth-order order conditions. Although RK methods based on the six-stage fifth-order order conditions have been widely studied and have continuing practical importance, prior to this thesis, the aforementioned complete solution to these order conditions has no published usage beyond the original series of publications by Cassity in the 1960s. The complete solution of six-stage fifth-order ERK order conditions published by Cassity in 1969 is not in a formulation that can easily be used for practical purposes, such as a software implementation. However, it is shown in this thesis that when the order conditions are solved and formulated appropriately using a computer algebra system (CAS), the generated code can be used for practical purposes and the complete solution is readily extended to ERK pairs. The condensed matrix form of the order conditions introduced by Cassity in 1969 is shown to be an ideal methodology, which probably has wider applicability, for solving order conditions using a CAS. The software package OCSage developed for this thesis, in order to solve the order conditions and study the properties of the resulting methods, is built on top of the Sage CAS. However, in order to effectively determine that the constructed ERK pairs without standard simplifying assumptions are in fact efficient by some well-defined criteria, the process of selecting the coefficients of ERK pairs is re-examined in conjunction with a sufficient amount of performance data. The pythODE software package developed for this thesis is used to generate a large amount of performance data from a large selection of candidate ERK pairs found using OCSage. In particular, it is shown that there is unlikely to be a well-defined methodology for selecting optimal pairs for general-purpose use, other than avoiding poor choices of certain properties and ensuring the error coefficients are as small as possible. However, for IVPs from celestial mechanics, there are obvious optimal pairs that have specific values of a small subset of the principal error coefficients (PECs). Statements seen in the literature that the best that can be done is treating all PECs equally do not necessarily apply to at least some broad classes of IVPs. By choosing ERK pairs based on specific values of individual PECs, not only are ERK pairs that are 20-30% more efficient than comparable published pairs found for test sets of IVPs from celestial mechanics, but the variation in performance between the best and worst ERK pairs that otherwise would seem to have similar properties is reduced from a factor of 2 down to as low as 15%. Based on observations of the small number of IVPs of other classes in common IVP test sets, there are other classes of IVPs that have different optimal values of the PECs. A more general contribution of this thesis is that it specifically demonstrates how specialized software tools and a larger amount of performance data than is typical can support novel empirical insights into numerical methods

Differential-Algebraic Equations

Differential-Algebraic Equations (DAE) are today an independent field of research, which is gaining in importance and becoming of increasing interest for applications and mathematics itself. This workshop has drawn the balance after about 25 years investigations of DAEs and the research aims of the future were intensively discussed