Magnus integrators for solving linear-quadratic differential games

We consider Magnus integrators to solve linear-quadratic N$N$-player differential games. These problems require to solve, backward in time, non-autonomous matrix Riccati differential equations which are coupled with the linear differential equations for the dynamic state of the game, to be integrated forward in time. We analyze different Magnus integrators which can provide either analytical or numerical approximations to the equations. They can be considered as time-averaging methods and frequently are used as exponential integrators. We show that they preserve some of the most relevant qualitative properties of the solution for the matrix Riccati differential equations as well as for the remaining equations. The analytical approximations allow us to study the problem in terms of the parameters involved. Some numerical examples are also considered which show that exponential methods are, in general, superior to standard methods.     

Evaluating the Evans function: Order reduction in numerical methods

We consider the numerical evaluation of the Evans function, a Wronskian-like determinant that arises in the study of the stability of travelling waves. Constructing the Evans function involves matching the solutions of a linear ordinary differential equation depending on the spectral parameter. The problem becomes stiff as the spectral parameter grows. Consequently, the Gauss-Legendre method has previously been used for such problems; however more recently, methods based on the Magnus expansion have been proposed. Here we extensively examine the stiff regime for a general scalar Schrödinger operator. We show that although the fourth-order Magnus method suffers from order reduction, a fortunate cancellation when computing the Evans matching function means that fourth-order convergence in the end result is preserved. The Gauss-Legendre method does not suffer from order reduction, but it does not experience the cancellation either, and thus it has the same order of convergence in the end result. Finally we discuss the relative merits of both methods as spectral tools.

     

Magnus’ expansion for time-periodic systems: Parameter-dependent approximations

Magnus’ expansion solves the nonlinear Hausdorff equation associated with a linear time-varying system of ordinary differential equations by forming the matrix exponential of a series of integrated commutators of the matrix-valued coefficient. Instead of expanding the fundamental solution itself, that is, the logarithm is expanded. Within some finite interval in the time variable, such an expansion converges faster than direct methods like Picard iteration and it preserves symmetries of the ODE system, if present. For time-periodic systems, Magnus expansion, in some cases, allows one to symbolically approximate the logarithm of the Floquet transition matrix (monodromy matrix) in terms of parameters. Although it has been successfully used as a numerical tool, this use of the Magnus expansion is new. Here we use a version of Magnus’ expansion due to Iserles [Iserles A. Expansions that grow on trees. Not Am Math Soc 2002;49:430–40], who reordered the terms of Magnus’ expansion for more efficient computation. Though much about the convergence of the Magnus expansion is not known, we explore the convergence of the expansion and apply known convergence estimates. We discuss the possible benefits to using it for time-periodic systems, and we demonstrate the expansion on several examples of periodic systems through the use of a computer algebra system, showing how the convergence depends on parameters.     

Adiabatic Integrators for Highly Oscillatory Second-Order Linear Differential Equations with Time-Varying Eigendecomposition

Numerical integrators for second-order differential equations with time-dependent high frequencies are proposed and analysed. We derive two such methods, called the adiabatic midpoint rule and the adiabatic Magnus method. The integrators are based on a transformation of the problem to adiabatic variables and an expansion technique for the oscillatory integrals. They can be used with far larger step sizes than those required by traditional schemes, as is illustrated by numerical experiments. We prove second-order error bounds with step sizes significantly larger than the almost-period of the fastest oscillations.AMS subject classification (2000) 65L05, 65L70.Received February 2004. Accepted February 2005. Communicated by Syvert Nørsett.     

Geometric Integration Methods that Preserve Lyapunov Functions

We consider ordinary differential equations (ODEs) with a known Lyapunov function V. To ensure that a numerical integrator reflects the correct dynamical behaviour of the system, the numerical integrator should have V as a discrete Lyapunov function. Only second-order geometric integrators of this type are known for arbitrary Lyapunov functions. In this paper we describe projection-based methods of arbitrary order that preserve any given Lyapunov function. AMS subject classification (2000) 65L05, 65L06, 65L20, 65P40     

Approximated solutions in rational form for systems of differential equations

In this paper we present a technique to study the existence of rational solutions for systems of differential equations — for an ordinary differential equation, in particular. The method is relatively straightforward; it is based on a rationality characterisation that involves matrix Padé approximants. It is important to note that, when the solution is rational, we use formal power series “without taking into account” their circle of convergence; at the end of this paper we justify this. We expound the theory for systems of linear first-order ordinary differential equations in the general case. However, the main ideas are applied in numerical resolution of partial differential equations. This revised version was published online in June 2006 with corrections to the Cover Date.     

New open modified Newton Cotes type formulae as multilayer symplectic integrators

New modified open Newton Cotes integrators are introduced in this paper. For the new proposed integrators the connection between these new algorithms, differential methods and symplectic integrators is studied. Much research has been done on one step symplectic integrators and several of them have obtained based on symplectic geometry. However, the research on multistep symplectic integrators is very poor. Zhu et al. [1] studied the well known open Newton Cotes differential methods and they presented them as multilayer symplectic integrators. Chiou and Wu [2] studied the development of multistep symplectic integrators based on the open Newton Cotes integration methods. In this paper we introduce a new open modified numerical method of Newton Cotes type and we present it as symplectic multilayer structure. The new obtained symplectic schemes are applied for the solution of Hamilton’s equations of motion which are linear in position and momentum. An important remark is that the Hamiltonian energy of the system remains almost constant as integration proceeds. We have applied also efficiently the new proposed method to a nonlinear orbital problem and an almost periodic orbital problem.     

Exponential Rosenbrock integrators for option pricing

In this paper, we are concerned with the time integration of differential equations modeling option pricing. In particular, we consider the Black-Scholes equation for American options. As an alternative to existing methods, we present exponential Rosenbrock integrators. These integrators require the evaluation of the exponential and related functions of the Jacobian matrix. The resulting methods have good stability properties. They are fully explicit and do not require the numerical solution of linear systems, in contrast to standard integrators. We have implemented some numerical experiments in Matlab showing the reliability of the new method.     

Solving ODEs arising from non‐selfadjoint Hamiltonian eigenproblems

We consider three numerical methods – one based on power series, one on the Magnus series and matrix exponentials, and one a library initial value code – for solving a linear system arising in non‐selfadjoint ODE eigenproblems. We show that in general, none of these methods has a cost or an accuracy which is uniform in the eigenparameter, but that for certain special types of problem, the Magnus method does yield eigenparameter‐uniform accuracy. This property of the Magnus method is explained by a trajectory‐shadowing result which, unfortunately, does not generalize to higher order Magnus type methods such as those in [11,12]. This revised version was published online in June 2006 with corrections to the Cover Date.     

Numerical integrators for quantum dynamics close to the adiabatic limit

Summary.  We study the numerical solution of singularly perturbed Schr?-dinger equations with time-dependent Hamiltonian. Based on a reformulation of the equations, we derive time-reversible numerical integrators which can be used with step sizes that are substantially larger than with traditional integration schemes. A complete error analysis is given for the adiabatic case. To deal with avoided crossings of energy levels, which lead to non-adiabatic behaviour, we propose an adaptive extension of the methods which resolves the sharp transients that appear in non-adiabatic state transitions. Received November 12, 2001 / Revised version received May 8, 2002 / Published online October 29, 2002 Mathematics Subject Classification (1991): 65L05, 65M15, 65M20, 65L70.     

On Galilean connections and the first jet bundle

We see how the first jet bundle of curves into affine space can be realized as a homogeneous space of the Galilean group. Cartan connections with this model are precisely the geometric structure of second-order ordinary differential equations under time-preserving transformations — sometimes called KCC-theory. With certain regularity conditions, we show that any such Cartan connection induces “laboratory” coordinate systems, and the geodesic equations in this coordinates form a system of second-order ordinary differential equations. We then show the converse — the “fundamental theorem” — that given such a coordinate system, and a system of second order ordinary differential equations, there exists regular Cartan connections yielding these, and such connections are completely determined by their torsion.     

On the Global Error of Discretization Methods for Highly-Oscillatory Ordinary Differential Equations

Commencing from a global-error formula, originally due to Henrici, we investigate the accumulation of global error in the numerical solution of linear highly-oscillating systems of the form y+ g(t)y = 0, where g(t) . Using WKB analysis we derive an explicit form of the global-error envelope for Runge-Kutta and Magnus methods. Our results are closely matched by numerical experiments.Motivated by the superior performance of Lie-group methods, we present a modification of the Magnus expansion which displays even better long-term behaviour in the presence of oscillations.     

Improved High Order Integrators Based on the Magnus Expansion

We build high order efficient numerical integration methods for solving the linear differential equation = A(t)X based on the Magnus expansion. These methods preserve qualitative geometric properties of the exact solution and involve the use of single integrals and fewer commutators than previously published schemes. Sixth- and eighth-order numerical algorithms with automatic step size control are constructed explicitly. The analysis is carried out by using the theory of free Lie algebras.     

Derivatives of the Evans function and (modified) Fredholm determinants for first order systems

The Evans function is a Wronskian type determinant used to detect point spectrum of differential operators obtained by linearizing PDEs about special solutions such as traveling waves, etc. This work is a sequel to the paper “Derivatives of (modified) Fredholm determinants and stability of standing and traveling waves”, published by F. Gesztesy, K. Zumbrun and the second author in J. Math. Pures Appl. 90 , 160–200 (2008), where the Evans and Jost functions for the Schrödinger equations have been considered. In the current work, we study the Evans function for the general case of linear ODE systems, and choose it to agree with the modified Fredholm determinant of the respective Birman‐Schwinger type integral operator. The Evans function is thus the determinant of the matrix composed of the so‐called generalized Jost solutions. These are the solutions of the homogeneous perturbed differential equation which are asymptotic to some reference solutions of the unperturbed equation. One of the main results of the current paper is a formula for the derivative of the Evans function for the first order systems. Its proof uses a matrix composed of the newly introduced modified Jost solutions. These are the solutions of an inhomogeneous perturbed differential equation with the inhomogeneous term constructed by means of the above‐mentioned generalized Jost solutions.     

Resolvent Krylov subspace approximation to operator functions

We consider the approximation of operator functions in resolvent Krylov subspaces. Besides many other applications, such approximations are currently of high interest for the approximation of φ-functions that arise in the numerical solution of evolution equations by exponential integrators. It is well known that Krylov subspace methods for matrix functions without exponential decay show superlinear convergence behaviour if the number of steps is larger than the norm of the operator. Thus, Krylov approximations may fail to converge for unbounded operators. In this paper, we analyse a rational Krylov subspace method which converges not only for finite element or finite difference approximations to differential operators but even for abstract, unbounded operators whose field of values lies in the left half plane. In contrast to standard Krylov methods, the convergence will be independent of the norm of the discretised operator and thus of the spatial discretisation. We will discuss efficient implementations for finite element discretisations and illustrate our analysis with numerical experiments.     

Yield condition for circular cylindrical shells made of a fiber-reinforced composite

A yield condition is obtained for circular cylindrical shells made of a definite class of fiber-reinforced composite material whose components possess plastic properties. It is shown that, in the plane of generalized stresses — the axial bending moment and the circumferential force (when the axial force is absent) — the yield curve consists of two linear and four curvilinear sections. By approximating the curvilinear sections, we get a piecewise linear yield condition described by a hexagon in the plane indicated. The nonlinear equations and the corresponding piecewise linear equations of the yield condition for particular cases are given in the form of tables. In solving specific boundary-value problems, we consider a circular cylindrical shell simply supported at its ends and loaded with a uniform internal pressure, for which the load-carrying capacity is determined in relation to the mechanical properties of composite components and some characteristic geometrical parameters. The results of numerical calculations are represented in the form of graphs. __________ Translated from Mekhanika Kompozitnykh Materialov, Vol. 42, No. 5, pp. 655–666, September–October, 2006.     

Taylor series method for dynamical systems with control: Convergence and error estimates

To optimize a complicated function constructed from a solution of a system of ordinary differential equations (ODEs), it is very important to be able to approximate a solution of a system of ODEs very precisely. The precision delivered by the standard Runge-Kutta methods often is insufficient, resulting in a “noisy function” to optimize. We consider an initial-value problem for a system of ordinary differential equations having polynomial right-hand sides with respect to all dependent variables. First we show how to reduce a wide class of ODEs to such polynomial systems. Using the estimates for the Taylor series method, we construct a new “aggregative” Taylor series method and derive guaranteed a priori step-size and error estimates for Runge-Kutta methods of order r. Then we compare the 8,13-Prince-Dormand’s, Taylor series, and aggregative Taylor series methods using seven benchmark systems of equations, including van der Pol’s equations, the “brusselator,” equations of Jacobi’s elliptic functions, and linear and nonlinear stiff systems of equations. The numerical experiments show that the Taylor series method achieves the best precision, while the aggregative Taylor series method achieves the best computational time. The final section of this paper is devoted to a comparative study of the above numerical integration methods for systems of ODEs describing the optimal flight of a spacecraft from the Earth to the Moon. __________ Translated from Sovremennaya Matematika i Ee Prilozheniya (Contemporary Mathematics and Its Applications), Vol. 24, Dynamical Systems and Optimization, 2005.     

A Magnus- and Fer-Type Formula in Dendriform Algebras

We provide a refined approach to the classical Magnus (Commun. Pure Appl. Math. 7:649–673, [1954]) and Fer expansion (Bull. Classe Sci. Acad. R. Belg. 44:818–829, [1958]), unveiling a new structure by using the language of dendriform and pre-Lie algebras. The recursive formula for the logarithm of the solutions of the equations X=1+λ a ≺ X and Y=1−λ Y ≻ a in A[[λ]] is provided, where (A,≺,≻) is a dendriform algebra. Then we present the solutions to these equations as an infinite product expansion of exponentials. Both formulae involve the pre-Lie product naturally associated with the dendriform structure. Several applications are presented.      

Algebraic independence of the values ofE-functions at singular points and Siegel’s conjecture

We prove a generalization of Shidlovskii’s theorem on the algebraic independence of the values ofE-functions satisfying a system of linear differential equations that is well known in the theory of transcendental numbers. We consider the case in which the values ofE-functions are taken at singular points of these systems. Using the obtained results, we prove Siegel’s conjecture that, for the case of first-order differential equations, anyE-function satisfying a linear differential equation is representable as a polynomial in hypergeometricE-functions. Translated fromMatematicheskie Zametki, Vol. 67, No. 2, pp. 174–190, February, 2000.     

High order local linearization methods: An approach for constructing A-stable explicit schemes for stochastic differential equations with additive noise

An approach for the construction of A-stable high order explicit strong schemes for stochastic differential equations (SDEs) with additive noise is proposed. We prove that such schemes also have the dynamical property that we call Random A-stability (RA-stability), which ensures that, for linear equations with stationary solutions, the numerical scheme has a random attractor that converges to the exact one as the step size decreases. Basically, the proposed integrators are obtained by splitting, at each time step, the solution of the original equation into two parts: the solution of a linear ordinary differential equation plus the solution of an auxiliary SDE. The first one is solved by the Local Linearization scheme in such a way that A-stability is guaranteed, while the second one is approximated by any extant scheme, preferably an explicit one that yields high order of convergence with low computational cost. Numerical integrators constructed in this way are called High Order Local Linearization (HOLL) methods. Various efficient HOLL schemes are elaborated in detail, and their performance is illustrated through computer simulations. Furthermore, mean-square convergence of the introduced methods is studied.