Large Systems Of Differential Equations Research Articles

Jacobian-free Locally Linearized Runge–Kutta method of Dormand and Prince for large systems of differential equations

This paper introduces Jacobian-free Locally Linearized Runge–Kutta formulas of Dormand and Prince for integrating large systems of initial value problems. With these new Jacobian-free formulas, an adaptive time-stepping variable-order scheme with regulated Krylov dimension is constructed with new developments in the computation of the phi-function action over vectors through Jacobian-free Krylov–Padé approximations, and in the procedures for controlling the approximation errors and for estimating the dimension of the Krylov subspaces. At each integration step, the implemented scheme performs only one Krylov subspace decomposition and a few computations of low dimensional exponential matrices, which contrasts with the implementation of other exponential-type integrators of high order. Regarding to existing schemes, the performed numerical study demonstrates a higher effectiveness of the constructed scheme in the integration of a variety of test equations.

Reproducibility of biophysical in silico neuron states and spikes from event-based partial histories.

Biophysically detailed simulations of neuronal activity often rely on solving large systems of differential equations; in some models, these systems have tens of thousands of states per cell. Numerically solving these equations is computationally intensive and requires making assumptions about the initial cell states. Additional realism from incorporating more biological detail is achieved at the cost of increasingly more states, more computational resources, and more modeling assumptions. We show that for both a point and morphologically-detailed cell model, the presence and timing of future action potentials is probabilistically well-characterized by the relative timings of a moderate number of recent events alone. Knowledge of initial conditions or full synaptic input history is not required. While model time constants, etc. impact the specifics, we demonstrate that for both individual spikes and sustained cellular activity, the uncertainty in spike response decreases as the number of known input events increases, to the point of approximate determinism. Further, we show cellular model states are reconstructable from ongoing synaptic events, despite unknown initial conditions. We propose that a strictly event-based modeling framework is capable of representing the complexity of cellular dynamics of the differential-equations models with significantly less per-cell state variables, thus offering a pathway toward utilizing modern data-driven modeling to scale up to larger network models while preserving individual cellular biophysics.

Lightning impulse modeling and simulation of dry-type and oil-immersed power- and distribution transformers

This paper presents in detail numerical methods and techniques for lightning impulse (LI) modeling and simulation of power and distribution transformers. The modeling methods are based on equivalent circuits of transformer winding entities resulting from the initial winding discretization determined by the required accuracy. The parameters of the equivalent circuit such as resistances and self- and mutual capacitances and inductances are obtained from field simulations (FEM). The circuit equations of the transformer’s equivalent circuit written in the state space form yield a large system of differential equations that is solved in time-domain by using the standard Runge-Kutta numerical integration technique. The obtained solution represents the voltage distribution over the winding in each moment of the LI-time (50μs). The results verification by comparison against measurements is presented in detail.

Analytic and numerical solutions to the seismic wave equation in continuous media.

This paper presents two approaches to mathematical modelling of a synthetic seismic pulse, and a comparison between them. First, a new analytical model is developed in two-dimensional Cartesian coordinates. Combined with an initial condition of sufficient symmetry, this provides a valuable check for the validity of the numerical method that follows. A particular initial condition is found which allows for a new closed-form solution. A numerical scheme is then presented which combines a spectral (Fourier) representation for displacement components and wave-speed parameters, a fourth-order Runge-Kutta integration method, and an absorbing boundary layer. The resulting large system of differential equations is solved in parallel on suitable enhanced performance desktop hardware in a new software implementation. This provides an alternative approach to forward modelling of waves within isotropic media which is efficient, and tailored to rapid and flexible developments in modelling seismic structure, for example, shallow depth environmental applications. Visual comparisons of the analytic solution and the numerical scheme are presented.

Model order reduction of thermo-mechanical models with parametric convective boundary conditions: focus on machine tools

Thermo-mechanical finite element (FE) models predict the thermal behavior of machine tools and the associated mechanical deviations. However, one disadvantage is their high computational expense, linked to the evaluation of the large systems of differential equations. Therefore, projection-based model order reduction (MOR) methods are required in order to create efficient surrogate models. This paper presents a parametric MOR method for weakly coupled thermo-mechanical FE models of machine tools and other similar mechatronic systems. This work proposes a reduction method, Krylov Modal Subspace (KMS), and a theoretical bound of the reduction error. The developed method addresses the parametric dependency of the convective boundary conditions using the concept of system bilinearization. The reduced-order model reproduces the thermal response of the original FE model in the frequency range of interest for any value of the parameters describing the convective boundary conditions. Additionally, this paper investigates the coupling between the reduced-order thermal system and the mechanical response. A numerical example shows that the reduced-order model captures the response of the original system in the frequency range of interest.

Lateral transfer in Stochastic Dollo models

Lateral transfer, a process whereby species exchange evolutionary traits through nonancestral relationships, is a frequent source of model misspecification in phylogenetic inference. Lateral transfer obscures the phylogenetic signal in the data as the histories of affected traits are mosaics of the overall phylogeny. We control for the effect of lateral transfer in a Stochastic Dollo model and a Bayesian setting. Our likelihood is highly intractable, as the parameters are the solution of a sequence of large systems of differential equations representing the expected evolution of traits along a tree. We illustrate our method on a data set of lexical traits in Eastern Polynesian languages, and obtain an improved fit over the corresponding model without lateral transfer.

Structure Preserving Model Reduction of Parametric Hamiltonian Systems

While reduced-order models (ROMs) have been popular for efficiently solving large systems of differential equations, the stability of reduced models over long-time integration presents challenges. We present a greedy approach for a ROM generation of parametric Hamiltonian systems that captures the symplectic structure of Hamiltonian systems to ensure stability of the reduced model. Through the greedy selection of basis vectors, two new vectors are added at each iteration to the linear vector space to increase the accuracy of the reduced basis. We use the error in the Hamiltonian due to model reduction as an error indicator to search the parameter space and identify the next best basis vectors. Under natural assumptions on the set of all solutions of the Hamiltonian system under variation of the parameters, we show that the greedy algorithm converges at an exponential rate. Moreover, we demonstrate that combining the greedy basis with the discrete empirical interpolation method also preserves the symplectic structure. This enables the reduction of the computational cost for nonlinear Hamiltonian systems. The efficiency, accuracy, and stability of this model reduction technique is illustrated through simulations of the parametric wave equation and the parametric Schrödinger equation.

Fourth-order two-stage explicit exponential integrators for time-dependent PDEs

Among the family of fourth-order time integration schemes, the two-stage Gauss–Legendre method, which is an implicit Runge–Kutta method based on collocation, is the only superconvergent. The computational cost of this implicit scheme for large systems, however, is very high since it requires solving a nonlinear system at every step. Surprisingly, in this work we show that one can construct and prove convergence results for exponential methods of order four which use two stages only. Specifically, we derive two new fourth-order two-stage exponential Rosenbrock schemes for solving large systems of differential equations. Moreover, since the newly schemes are not only superconvergent but also fully explicit, they turn out to be very competitive compared to the two-stage Gauss–Legendre method as well as other fourth-order time integration schemes. Numerical experiments are given to demonstrate the efficiency of the new integrators.

Coupling Elastic Bodies with an Enhanced Craig‐Bampton‐like Scheme

AbstractElastic multibody systems are widely used as description of mechanical systems that consist of separable components. These components are represented by elastic bodies that interact by the means of force and damping elements as well as boundary conditions. The description of the elasticity is usually provided by a linear finite‐element model. If the bodies are densely meshed, a large system of differential equations is obtained. The numerical effort can be reduced by model order reduction. This contribution presents reduction approaches that are especially well suited for the elastic multibody system framework. (© 2015 Wiley‐VCH Verlag GmbH & Co. KGaA, Weinheim)

Flexible exponential integration methods for large systems of differential equations

In this paper, we describe a flexible variant of exponential integration methods for large systems of differential equations. This version possesses the flexibility and generality which allows to further exploit the special structure of the system. By using modified B-series and bi-coloured rooted trees, we can derive the general structure of the classical order conditions for these schemes. Some numerical schemes are constructed and the order conditions are derived. Numerical experiments with reaction-diffusion type problems are included.

Application of approximate matrix factorization to high order linearly implicit Runge–Kutta methods

Linearly implicit Runge–Kutta methods with approximate matrix factorization can solve efficiently large systems of differential equations that have a stiff linear part, e.g. reaction–diffusion systems. However, the use of approximate factorization usually leads to loss of accuracy, which makes it attractive only for low order time integration schemes. This paper discusses the application of approximate matrix factorization with high order methods; an inexpensive correction procedure applied to each stage allows to retain the high order of the underlying linearly implicit Runge–Kutta scheme. The accuracy and stability of the methods are studied. Numerical experiments on reaction–diffusion type problems of different sizes and with different degrees of stiffness illustrate the efficiency of the proposed approach.

PDE approximation of large systems of differential equations

A large system of ordinary differential equations is approximated by a parabolic partial differential equation with dynamic boundary condition and a different one with Robin boundary condition. Using the theory of differential operators with Wentzell boundary conditions and similar theories, we give estimates on the order of approximation. The theory is demonstrated on a voter model where the Fourier method applied to the PDE is of great advantage.

Rosenbrock--Krylov Methods for Large Systems of Differential Equations

This paper develops a new class of Rosenbrock-type integrators based on a Krylov space solution of linear systems. The new family, called Rosenbrock--Krylov (Rosenbrock-K), is well suited for solving large-scale systems of ODEs or semidiscrete PDEs. The time discretization and the Krylov space approximation are treated as a single computational process, and the Krylov space properties are an integral part of the new Rosenbrock-K order condition theory developed herein. Consequently, Rosenbrock-K methods require a small number of basis vectors determined solely by the temporal order of accuracy. The subspace size is independent of the ODE under consideration, and there is no need to monitor the errors in linear system solutions at each stage. Numerical results show favorable properties of Rosenbrock-K methods when compared to current Rosenbrock and Rosenbrock-W schemes.

Calibration of the Stochastic Multicloud Model Using Bayesian Inference

The stochastic multicloud model (SMCM) was recently developed [B. Khouider, J. Biello, and A. J. Majda, Commun. Math. Sci., 8 (2010), pp. 187--216] to represent the missing variability in general circulation models due to unresolved features of organized tropical convection. This research aims at finding a robust calibration methodology for the SMCM to estimate key model parameters from data. We formulate the calibration problem within a Bayesian framework to derive the posterior distribution over the model parameters. The main challenge here is due to the likelihood function which requires solving a large system of differential equations (the Kolmogorov equations) as many times as there are data points, which is prohibitive in terms of both computation time and storage requirements. The most attractive numerical techniques to compute the transient solutions to large Markov chains are based on matrix exponentials, but none is unconditionally acceptable for all classes of problems. We develop a parallel ve...

On the use of exponential time integration methods in atmospheric models

Exponential integration methods offer a highly accurate approach to the time integration of large systems of differential equations. In recent years, they have attracted increased attention in a number of diverse fields due to advances in their computational efficiency. This has been as a result of the use of Krylov subspace methods for the approximation of the matrix exponentials which typically arise. In this work, we investigate the potential of exponential integration methods for use in atmospheric models. Two schemes are implemented in a shallow water model and tested against reference explicit and semi-implicit methods. In a number of experiments with standard test cases, the exponential methods are found to yield very accurate solutions with time-steps far longer than even the semi-implicit method allows. The relative efficiency of the exponential integrators, which depends mainly on the choice of the specific algorithm used for the calculation of the matrix exponent, is also discussed. The future work aimed at further improvements of the proposed methodology is outlined.

Performance assessment of exponential Rosenbrock methods for large systems of ODEs

This article studies time integration methods for stiff systems of ordinary differential equations of large dimension. For such problems, implicit methods generally outperform explicit methods because the step size is usually less restricted by stability constraints. Recently, however, a family of explicit methods, called exponential integrators, have become popular for large stiff problems due to their favourable stability properties and the rapid convergence of non-preconditioned Krylov subspace methods for computing matrix-vector products involving exponential-like functions of the Jacobian matrix. In this article, we implement the so-called exponential Rosenbrock methods using Krylov subspaces. Numerical experiments on a challenging real-world test problem reveal that these methods are a promising preconditioner-free alternative to well-established approaches based on preconditioned Newton--Krylov implementations of the backward differentiation formulas. References E. J. Carr, T. J. Moroney and I. W. Turner. Efficient simulation of unsaturated flow using exponential time integration. Appl. Math. Comput., 217(14): 6587--6596, 2011. doi:10.1016/j.amc.2011.01.041 E. J. Carr, I. W. Turner and M. Ilic. Krylov subspace approximations for the exponential Euler method: error estimates and the harmonic Ritz approximant. ANZIAM Journal, 52: C612--C627, 2011. http://journal.austms.org.au/ojs/index.php/ANZIAMJ/article/view/3938 E. J. Carr, I. W. Turner and P. Perre. A variable step size Jacobian-free exponential integrator for simulating transport in heterogeneous porous media: application to wood drying. J. Comput. Phys., 233: 66-82, 2013. doi:10.1016/j.jcp.2012.07.024 Sundials user documentation. https://computation.llnl.gov/casc/sundials/documentation/documentation.html M. Hochbruck, C. Lubich and H. Selhofer. Exponential integrators for large systems of differential equations. SIAM J. Sci. Comput., 19(5): 1552--1574, 1998. doi:10.1137/S1064827595295337 M. Hochbruck, A. Ostermann and J. Schweitzer. Exponential Rosenbrock-type methods. SIAM J. Numer. Anal., 47(1): 786--803, 2009. doi:10.1137/080717717 C. T. Kelley. Solving nonlinear equations with Newton's method. SIAM, USA, 2003. D. A. Knoll and D. E. Keyes. Jacobian-free Newton--Krylov methods: a survey of approaches and applications. J. Comput. Phys., 193(2): 357--397, 2004. doi:10.1016/j.jcp.2003.08.010 M. Tokman. Efficient integration of large stiff systems of ODEs with exponential propagation iteration (EPI) methods. J. Comput. Phys., 213: 748--776, 2006. doi:10.1016/j.jcp.2005.08.032

Epidemic Spread and Variation of Peak Times in Connected Regions Due to Travel-Related Infections---Dynamics of an Antigravity-Type Delay Differential Model

National boundaries have never prevented infectious diseases from reaching distant territories; however, the speed at which an infectious agent can spread around the world via the global airline transportation network has significantly increased during recent decades. We introduce an SEAIR-based, antigravity model to investigate the spread of an infectious disease in two regions which are connected by transportation. As a submodel, an age-structured system is constructed to incorporate the possibility of disease transmission during travel, where age is the time elapsed since the start of the travel. The model is equivalent to a large system of differential equations with dynamically defined delayed feedback. After describing fundamental but biologically relevant properties of the system, we detail the calculation of the basic reproduction number and obtain disease transmission dynamics results in terms of $\mathcal{R}_0$. We parametrize our model for influenza and use real demographic and air travel data ...

Discrete and Continuum Relaxation Dynamics of Faceted Crystal Surface in Evaporation Models

We study linkages of two scales in the relaxation of an axisymmetric crystal with a facet in evaporation-condensation kinetics. The macroscale evolution is driven by the motion of concentric circular, repulsively interacting line defects (steps) which exchange atoms with the vapor. At the microscale, the step velocity is proportional to the variation of the total step free energy, leading to large systems of differential equations for the step radii. We focus on two step flow models. In one model (called M1) the discrete mobility is simply proportional to the upper-terrace width; in another model (M2) the mobility is altered by an extra geometric factor. By invoking self-similarity at long time, we numerically demonstrate that (i) in M1, discrete slopes follow closely a continuum thermodynamics approach with “natural boundary conditions” at the facet edge; (ii) in contrast, predictions of M2 deviate from results of the above continuum approach; and (iii) this discrepancy can be eliminated via a continuum boundary condition with a geometry-induced jump for top-step collapses. At the macroscale, both step models give rise to free-boundary problems for a second-order, parabolic partial differential equation, which we study via the subgradient formalism. We discuss the interpretation of the facet height as shock and prove convergence of the solution of each discrete scheme to the (weak) entropy solution of a conservation law if steps do not interact.

Multiregional SIR Model with Infection during Transportation

We present a general epidemic model to describe the spread of an infectious disease in several regions connected by transportation. We take into account that infected individuals not only carry the disease to a new place while traveling from one region to another, but transmit the disease during travel as well. We obtain that a model structured by travel time is equivalent to a large system of differential equations with multiple delays. By showing the local Lipschitz property of the dynamically defined delayed feedback function, we obtain existence and uniqueness of solutions of the system.

Krylov subspace approximations for the exponential Euler method: error estimates and the harmonic Ritz approximant

We study Krylov subspace methods for approximating the matrix-function vector product $\varphi(tA)b$ where $\varphi(z) = [\exp(z)-1]/z$. This product arises in the numerical integration of large stiff systems of differential equations by the Exponential Euler Method, where $A$~is the Jacobian matrix of the system. Recently, this method has found application in the simulation of transport phenomena in porous media within mathematical models of wood drying and groundwater flow. We develop an a posteriori upper bound on the Krylov subspace approximation error and provide a new interpretation of a previously published error estimate. This leads to an alternative Krylov approximation to $\varphi(tA)b$, the so-called Harmonic Ritz approximant, which we find does not exhibit oscillatory behaviour of the residual error. References E. J. Carr, T. J. Moroney and I. W. Turner. Efficient simulation of unsaturated flow using exponential time integration. Appl. Math. Comput., 217(14):6587--6596, 2011. doi:10.1016/j.amc.2011.01.041 E. J. Carr, I. W. Turner and P. Perre. A Jacobian-free exponential integrator for simulating transport in heterogeneous porous media: application to the drying of softwood, submitted for publication. E. J. Carr, I. W. Turner and P. Perre. A new control-volume finite-element scheme for heterogeneous porous media:application to the drying of softwood. Chem. Eng. Technol., 34(7):1143--1150, 2011. doi:10.1002/ceat.201100060 E. Celledoni and I. Moret. A Krylov projection method for systems of ODEs. Appl. Numer. Math., 24(2--3):365--378, 1997. doi:10.1016/S0168-9274(97)00033-0 N. J. Higham. Functions of matrices: theory and computation. SIAM, Philadelphia, PA, USA, 2008. M. Hochbruck, M. E. Hochstenbach. Subspace extraction for matrix functions, submitted for publication. http://www.win.tue.nl/ hochsten/pdf/funext.pdf M. Hochbruck, C. Lubich and H. Selhofer. Exponential integrators for large systems of differential equations. SIAM J. Sci. Comput., 19(5):1552--1574, 1998. doi:10.1137/S1064827595295337 B. V. Minchev and W. M. Wright. A review of exponential integrators for first order semi-linear problems. Numerics No. 2/05, Norwegian University of Science and Technology, Trondheim, Norway, 2005. http://www.math.ntnu.no/preprint/numerics/2005/N2-2005.ps P. Perre and I. Turner. A heterogeneous wood drying computational model that accounts for material property variation across growth rings. Chem. Eng. J., 86(1--2):117--131, 2002. doi:10.1016/S1385-8947(01)00270-4