Elsevier

Applied Numerical Mathematics

Volume 137, March 2019, Pages 184-202
Applied Numerical Mathematics

Asymptotic-numerical solvers for highly oscillatory second-order differential equations

https://doi.org/10.1016/j.apnum.2018.11.004Get rights and content

Abstract

In this paper, we propose an approach of combination of asymptotic and numerical techniques to solve highly oscillatory second-order initial value problems. An asymptotic expansion of the solution is derived in inverse of powers of the oscillatory parameter, which develops on two time scales, a slow time t and a fast time τ=ωt. The truncation with the first few terms of the expansion results in a very effective method of discretizing the highly oscillatory differential equation and becomes more accurate when the oscillatory parameter increases. Numerical examples show that our proposed asymptotic-numerical solver is efficient and accurate for highly oscillatory problems.

Introduction

Highly oscillatory problems have appeared in many fields such as celestial mechanics, chemistry, biology, classical and quantum mechanics, and engineering. Due to the oscillatory behaviour of the solutions of such differential equation, the use of standard methods of numerical ODEs (such as Runge–Kutta) imposes an exceedingly small step size which is both too expensive for implementation and leads to an accumulation of round-off error due to the large number of steps needed to integrate the ODE in a given interval. The design and analysis of numerical methods for highly oscillatory problems has received a great deal of attention in the past decades. In this paper, we consider highly oscillatory problems of the form{d2xdt2(t)+ω2x(t)=g(t,x(t)),t(0,T],x(0)=x0,dxdt(0)=x˙0, where x(t):[0,T]R, g(t,x):(0,T]×RR is a real analytic function, and ω1 represents the spring constant. The system is highly oscillatory problem for such values, and this description is frequently used in an extended sense to characterize problems of a highly oscillatory differential equations. Here we assume that the initial values satisfy |x0|+ω1|x˙0|E, where E is a constant which is independent of ω. This paper only focuses on scalar problems, but the proposed approach can be applied to high dimensional systems. However, it will become more complicated for solving high dimensional systems of second order ODEs, like those arising from the semidiscretization with respect to the space variable of partial differential equations.

Miranker and van Veldhuizen [25] first proposed a representation of the solution to (1.1) as a modulated Fourier expansionx(t)=k=eikωtzk(t), where the coefficient functions zk(t) are smoothly varying functions (i.e., their derivatives are bounded and independent of ω). Each envelope zk(t) is obtained by solving a system of nonlinear highly oscillatory ODEs. They computed these envelopes numerically and used them to approximate the solution. As a widely analysed and essential tool developed in the past, modulated Fourier expansions have been found useful to explain various long-time phenomena in both continuous and discretized oscillatory Hamiltonian systems, ordinary differential equations as well as partial differential equations, and also regarded as the basic ingredient in the heterogeneous multiscale method. In [4], [17], this technique of modulated Fourier expansions has been developed as a tool for gaining insight into the long-time behaviour of Hamiltonian systems with highly oscillatory solutions and analyzing the long-time behaviour of numerical integrators when the time step is not small compared with the value of ω1. Cohen, Hairer and Lubich [5] used a modulated Fourier expansion in time to show long-time near conservation of the harmonic actions associated with spatial Fourier modes along the solutions of nonlinear wave equations with small initial data. Hairer and Lubich [18] employed the modulated Fourier expansions to explain some of the phenomena and theoretical results on the long-time energy behaviour of continuous and discretized oscillatory systems. Sanz-Serna [28] showed that the modulated Fourier expansion approach can be advantageously used to understand and analyze the Heterogenous multiscale methods.

For differential equations with highly oscillatory forcing terms in the case of modulated Fourier oscillator, efficient asymptotic-numerical solvers based on the asymptotic expansion in inverse powers of the oscillatory parameters and its truncation have been proposed in [6], [7], [8], [9], [10]. Condon, Iserles and Nørsett [12] made use of the variation-of-constants representation and computed the highly oscillatory integral therein by modern quadrature methods for differential equations with general highly oscillatory forcing terms. Recently, Condon, Deaño, Iserles and Kropielnicka [11] proposed efficient computation of delay differential equations with highly oscillatory terms based on the asymptotic expansion of the solution. For other types of numerical methods for highly oscillatory systems, we refer to [1], [2], [3], [13], [14], [15], [16], [19], [20], [21], [22], [23], [24], [26], [27], [29], [30] and the references therein.

In the present paper, we aim to propose an alternative different approach from the method of envelopes [25] to obtain an asymptotic expansion for the solution of (1.1) and design an efficient numerical solver to approximate the solution by truncating the first few terms of such asymptotic expansion. This approach features two fundamental advantages with respect to standard numerical ODE solvers. Firstly, the construction of the numerical solution is more efficient when the system is highly oscillatory and increasing the frequency will make our method more accurate. Secondly, the cost of the computation is essentially independent of the oscillatory parameter, in particular, one only needs to solve non-oscillatory ODEs whose right-hand side function depends on t only. This paper is organized as follows. In Section 2, we transform (1.1) into a first-order system, present the construction of the asymptotic expansion, and verify the asymptotic expansion with necessary uniform bounds of its coefficients and of the remainder term. In Section 3, we propose an asymptotic-numerical solver to obtain the approximation by truncating the first few terms of the asymptotic expansion. In Section 4, numerical experiments are carried out to demonstrate the efficiency and accuracy of our proposed method. Concluding remarks are given in Section 5.

Section snippets

Construction

We first transform (1.1) into a first-order system, instead of dealing with the original one. Let y(t)=ω1x˙(t). Equation (1.1) is equivalent to the following first-order system{(dxdt(t)dydt(t))=ωA(x(t)y(t))+(0ω1g(t,x(t))),(x(0)y(0))=(x0ω1x˙0), whereA=(0110) is a skew-symmetric matrix. Setu(t)=eiωBtU(x(t)y(t)), whereU=12(1ii1),B=(1001), “i” is the imaginary unit and U is the conjugate transpose of U. Usingx(t)=12(eiωt,ieiωt)u(t) and inserting u(t) into (2.1), we obtaindudt(t)=ω12(ieiω

Asymptotic-numerical solver

The asymptotic expansion of the solution of (1.1) has been established. The numerical implementation of the asymptotic expansion requires three different truncations. Firstly, we need to replace the right-hand side of (2.40) by a finite sum,x(t)12(eiωt,ieiωt)u0+12(eiωt,ieiωt)s=1Rωsk=eikωtps,k(t). Note that fairly small values of R are perfectly adequate for the case ω1. Secondly, we need to restrict the calculation of k=eikωtps,k(t) to the sum of a finite terms k=ρsρseikωtps,k(t)

Numerical experiments

In this section, we perform some experiments to demonstrate the efficiency and accuracy of our proposed asymptotic solver. In all cases, the exact solution may be analytically available or either computed numerically with standard Matlab routine ODE45 up to prescribed accuracy (AbsTOL=RelTOL=1013). We will provide the approximation given by the first few terms (R=1 and 2) of our proposed asymptotic-numerical solver (denoted by AF1 and AF2). For comparison, we use the Gautschi-type method with

Conclusions

We propose and verify that the solution of the highly oscillatory problem (1.1) can be written as a series in inverse powers of the oscillatory parameter ω, featuring modulated Fourier series in the expansion coefficients. The asymptotic theory yields numerical approximation by truncating the asymptotic expansion with the first few terms. Increasing ω will benefit the asymptotic-numerical solver, since the approximation with a fixed number of terms will be more accurate, while the computational

References (30)

  • M. Condon et al.

    Asymptotic solvers for ordinary differential equations with multiple frequencies

    Sci. China Math.

    (2015)
  • M. Condon et al.

    On highly oscillatory problems arising in electronic engineering

    Math. Model. Numer. Anal.

    (2009)
  • M. Condon et al.

    On second-order differential equations with highly oscillatory forcing terms

    Proc. R. Soc. Lond., Ser. A, Math. Phys. Eng. Sci.

    (2010)
  • M. Condon et al.

    Asymptotic solvers for oscillatory systems of differential equations

    SeMA J.

    (2011)
  • M. Condon et al.

    Efficient computation of delay differential equations with highly oscillatory terms

    ESAIM: Math. Model. Numer. Anal.

    (2012)
  • Cited by (0)

    1

    The work of this author is supported in part by the National Natural Science Foundation of China under Grant No. 11701378.

    2

    The work of this author is supported in part by E-Institutes of Shanghai Municipal Education Commission under Grant No. E03004, the National Natural Science Foundation of China under Grant Nos. 11671266 and 11871343, and the Natural Science Foundation of Shanghai under Grant No. 16ZR1424900.

    View full text