The numerical solution of nonlinear two-point boundary value problems using iterated deferred correction - a survey

The use of iterated deferred correction has proved to be a very efficient approach to the numerical solution of general first order systems of nonlinear two-point boundary value problems. In particular the two high order codes TWPBVP.f, based on mono-implicit Runge-Kutta (MIRK) formulae, and TWPBVPL.f based on Lobatto Runge-Kutta formulae as well as the continuation codes ACDC.f and COLMOD.f are now widely used. In this survey we describe some of the problems involved in the derivation of efficient deferred correction schemes. In particular we consider the construction of high order methods which preserve the stability of the underlying formulae, the choice of the mesh choosing algorithm which is based both on local accuracy and conditioning, and the computation of continuous solutions

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

Approximate solution of ordinary differential equations and their systems through discrete and continuous embedded Runge-Kutta formulae and upgrading of their order

AbstractThe eight main contributions of the author to the field of approximate solutions of ordinary differential equations described herein are all application-oriented, with the purposes of simplification and the increase in efficiency and effectiveness of the Runge-Kutta processes generated. They range from the determination of an initial trial step-size to be adopted to expedite the approximation process through embedded Runge-Kutta algorithms to a more recent procedure for upgrading the order of Runge-Kutta processes. These contributions encompass all classes of differential equations of all orders, such as explicit, implicit, single or systems, and their treatment by Runge-Kutta processes of scalar or vector type (with the related equivalence conditions), of discrete or continuous kind, including the computer derivations of nonlinear algebraic equations associated with the Runge-Kutta processes. Specifically, the author developed the first fifth order Runge-Kutta formulae with fourth order embedded and the first C1 approximate solution through interpolation and Runge-Kutta formulae, which he improved by developing C1 embeddings with Runge-Kutta formulae without the use of interpolative techniques

Jet schemes for advection problems

We present a systematic methodology to develop high order accurate numerical approaches for linear advection problems. These methods are based on evolving parts of the jet of the solution in time, and are thus called jet schemes. Through the tracking of characteristics and the use of suitable Hermite interpolations, high order is achieved in an optimally local fashion, i.e. the update for the data at any grid point uses information from a single grid cell only. We show that jet schemes can be interpreted as advect-and-project processes in function spaces, where the projection step minimizes a stability functional. Furthermore, this function space framework makes it possible to systematically inherit update rules for the higher derivatives from the ODE solver for the characteristics. Jet schemes of orders up to five are applied in numerical benchmark tests, and systematically compared with classical WENO finite difference schemes. It is observed that jet schemes tend to possess a higher accuracy than WENO schemes of the same order.Comment: 26 pages, 6 figures; based on a poster presentation at the conference Fluid Dynamics, Analysis, and Numerics 201

DENSERKS: Fortran sensitivity solvers using continuous, explicit Runge-Kutta schemes

DENSERKS is a Fortran sensitivity equation solver package designed for integrating models whose evolution can be described by ordinary differential equations (ODEs). A salient feature of DENSERKS is its support for both forward and adjoint sensitivity analyses, with built-in integrators for both first and second order continuous adjoint models. The software implements explicit Runge-Kutta methods with adaptive timestepping and high-order dense output schemes for the forward and the tangent linear model trajectory interpolation. Implementations of six Runge-Kutta methods are provided, with orders of accuracy ranging from two to eight. This makes DENSERKS suitable for a wide range of practical applications. The use of dense output, a novel approach in adjoint sensitivity analysis solvers, allows for a high-order cost-effective interpolation. This is a necessary feature when solving adjoints of nonlinear systems using highly accurate Runge-Kutta methods (order five and above). To minimize memory requirements and make long-time integrations computationally efficient, DENSERKS implements a two-level checkpointing mechanism. The code is tested on a selection of problems illustrating first and second order sensitivity analysis with respect to initial model conditions. The resulting derivative information is also used in a gradient-based optimization algorithm to minimize cost functionals dependent on a given set of model parameters

Superconvergant interpolants for the collocation solution of boundary value ordinary differential equations

Publisher's version/PDFA long-standing open question associated with the use of collocation methods for boundary value ordinary differential equations is concerned with the development of a high order continuous solution approximation to augment the high order discrete solution approximation, obtained at the mesh points which subdivide the problem interval. It is well known that the use of collocation at Gauss points leads to solution approximations at the mesh points for which the global error is O(h[superscript 2k]), where k is the number of collocation points used per subinterval and h is the subinterval size. This discrete solution is said to be superconvergent. The collocation solution also yields a C[superscript 0] continuous solution approximation that has a global error of O(h[supercript k+1]). In this paper, we show how to efficiently augment the superconvergent discrete collocation solution to obtain C[superscript 1] continuous "superconvergent" interpolants whose global errors are O(h[superscript 2k]). The key ideas are to use the theoretical framework of continuous Runge-Kutta schemes and to augment the collocation solution with inexpensive monoimplicit Runge-Kutta stages. Specific schemes are derived for k = 1, 2, 3, and 4. Numerical results are provided to support the theoretical analysis

A delay differential model of ENSO variability: Parametric instability and the distribution of extremes

We consider a delay differential equation (DDE) model for El-Nino Southern Oscillation (ENSO) variability. The model combines two key mechanisms that participate in ENSO dynamics: delayed negative feedback and seasonal forcing. We perform stability analyses of the model in the three-dimensional space of its physically relevant parameters. Our results illustrate the role of these three parameters: strength of seasonal forcing $b$, atmosphere-ocean coupling $\kappa$, and propagation period $\tau$ of oceanic waves across the Tropical Pacific. Two regimes of variability, stable and unstable, are separated by a sharp neutral curve in the $(b,\tau)$ plane at constant $\kappa$. The detailed structure of the neutral curve becomes very irregular and possibly fractal, while individual trajectories within the unstable region become highly complex and possibly chaotic, as the atmosphere-ocean coupling $\kappa$ increases. In the unstable regime, spontaneous transitions occur in the mean ``temperature'' ({\it i.e.}, thermocline depth), period, and extreme annual values, for purely periodic, seasonal forcing. The model reproduces the Devil's bleachers characterizing other ENSO models, such as nonlinear, coupled systems of partial differential equations; some of the features of this behavior have been documented in general circulation models, as well as in observations. We expect, therefore, similar behavior in much more detailed and realistic models, where it is harder to describe its causes as completely.Comment: 22 pages, 9 figure

Adaptive time-integration for goal-oriented and coupled problems

We consider efficient methods for the partitioned time-integration of multiphysics problems, which commonly exhibit a multiscale behavior, requiring independent time-grids. Examples are fluid structure interaction in e.g., the simulation of blood-flow or cooling of rocket engines, or ocean-atmosphere-vegetation interaction. The ideal method for solving these problems allows independent and adaptive time-grids, higher order time-discretizations, is fast and robust, and allows the coupling of existing subsolvers, executed in parallel. We consider Waveform relaxation (WR) methods, which can have all of these properties. WR methods iterate on continuous-in-time interface functions, obtained via suitable interpolation. The difficulty is to find suitable convergence acceleration, which is required for the iteration converge quickly. We develop a fast and highly robust, second order in time, adaptive WR method for unsteady thermal fluid structure interaction (FSI), modelled by heterogeneous coupled linear heat equations. We use a Dirichlet-Neumann coupling at the interface and an analytical optimal relaxation parameter derived for the fully-discrete scheme. While this method is sequential, it is notably faster and more robust than similar parallel methods.We further develop a novel, parallel WR method, using asynchronous communication techniques during time-integration to accelerate convergence. Instead of exchanging interpolated time-dependent functions at the end of each time-window or iteration, we exchange time-point data immediately after each timestep. The analytical description and convergence results of this method generalize existing WR theory.Since WR methods allow coupling of problems in a relative black-box manner, we developed adapters to PDE-subsolvers implemented using DUNE and FEniCS. We demonstrate this coupling in a thermal FSI test case.Lastly, we consider adaptive time-integration for goal-oriented problems, where one is interested in a quantity of interest (QoI), which is a functional of the solution. The state-of-the-art method is the dual-weighted residual (DWR) method, which is extremely costly in both computation and implementation. We develop a goal oriented adaptive method based on local error estimates, which is considerably cheaper in computation. We prove convergence of the error in the QoI for tolerance to zero under a controllability assumption. By analyzing global error propagation with respect to the QoI, we can identify possible issues and make performance predictions. Numerical results verify these results and show our method to be more efficient than the DWR method

The Physics of Financial Networks

As the total value of the global financial market outgrew the value of the real economy, financial institutions created a global web of interactions that embodies systemic risks. Understanding these networks requires new theoretical approaches and new tools for quantitative analysis. Statistical physics contributed significantly to this challenge by developing new metrics and models for the study of financial network structure, dynamics, and stability and instability. In this Review, we introduce network representations originating from different financial relationships, including direct interactions such as loans, similarities such as co-ownership and higher-order relations such as contracts involving several parties (for example, credit default swaps) or multilayer connections (possibly extending to the real economy). We then review models of financial contagion capturing the diffusion and impact of shocks across each of these systems. We also discuss different notions of ‘equilibrium’ in economics and statistical physics, and how they lead to maximum entropy ensembles of graphs, providing tools for financial network inference and the identification of early-warning signals of system-wide instabilities