Solving the wave equation with multifrequency oscillations

We explore a new asymptotic-numerical solver for the time-dependent wave equation with an interaction term that is oscillating in time with a very high frequency. The method involves representing the solution as an asymptotic series in inverse powers of the oscillation frequency. Using the new scheme, high accuracy is achieved at a low computational cost. Salient features of the new approach are highlighted by a numerical example.

1. Introduction. In this paper, we consider the time-dependent wave equation [2] where u(x, t) is the wavefunction and g ω (x, t) is a real-valued highly oscillatory term. The initial and boundary conditions are Such equations are considered when computing scattering frequencies [6,11,12]. Highly oscillatory interaction terms present a difficulty for numerical simulations as a small time step and fine space discretisation are typically required to obtain an accurate solution. The present contribution involves representing the solution as an asymptotic series in inverse powers of the temporal oscillation frequency of g ω (x, t). The coefficients of the terms in the series are independent of the oscillation frequency and hence, the computational effort in determining these coefficients and the associated asymptotic series is greatly reduced. The coefficients are determined from either recursion or partial differential equations and the solutions of the latter do not depend on the oscillation frequency of g ω (x, t). Consequently, the proposed method achieves high accuracies despite using extremely large time steps. In principle, it might appear that the error of the computation is independent of the frequency of oscillation but, actually, the situation is often even better! The more rapid the oscillation, the smaller the error of the numerical scheme. An example illustrates salient features of the underlying method. Before we describe our numerical approach, it is important to discuss briefly the well posedness of (1). This is a linear hyperbolic equation and, in general, we expect it to be well posed, but it is always a good idea to confirm this.
Letting v(x, t) = ∂ t u(x, t), we obtain a first-order hyperbolic system with the above initial and boundary conditions. Denote by E the semigroup associated with the wave equation -in other words, we consider the solution of the standard wave equation with the same initial conditions and zero Dirichlet boundary conditions by we recall that E(t) ≤ 1 in the H 1 Sobolev norm on u, which is identical to the standard L 2 norm on u v .
Applying the Duhamel principle to (2), we have therefore, by the integral form of the Grönwall Lemma and bearing in mind that E(t) ≤ 1, which, provided that max t≥0 g ω ( · , t) L∞(−L,L) is bounded, demonstrates that the solution of (2) is uniformly L 2 -bounded in every compact interval t ∈ [0, T ]. Hence well posedness. The approach of this paper is a mixture of numerical and asymptotic methodology. It continues and expands the scope of previous work in this area, both on ODEs and evolutionary PDEs [1,4,5,8,13].
A word about the stability of the underlying asymptotic-numerical approachsomething essential in the context of the Lax Equivalence Theorem. The expansion itself cannot be a source of instability. Yet, once it is solved by a time-stepping numerical method (which would be the usual state of affairs) its stability is inherited from the method in question.
2. The proposed solution. We consider equation (1) with where d1 d2 / ∈ Q and d 1 ω 1, d 2 ω 1 are the frequencies of interest. We assume that m,n |g m,n | is bounded. As a matter of principle, we could have added more frequencies, d 1 ω, d 2 ω, . . . , d q ω, say, but the current setup is sufficient to highlight all the salient points.
To expand asymptotically, we assume that there exist functions p 0,0 (x, t) and p r,m,n (x, t), r ≥ 2, n ∈ Z such that In other words, we expand u in modulated Fourier expansion [9]. Note one potential problem with our expansion: the mechanism used to produce the p r,m,n s recursively might lead to resonances because of the occurrence of small denominators. This is well known in asymptotic analysis and, at least in principle, might restrict the range of r ≥ 2 and m, n ∈ Z. Yet, given that a realistic numerical method, as described in Section 3, necessarily restricts this range, we do not anticipate this to be a substantive issue.
Differentiating term-by-term, we have We further observe that Once these equations are substituted into (1), we have We next define coefficients in two levels. Firstly we consider orders of magnitude (inverse powers of ω -signified by the values of r), and then frequencies (values of m and n) within each order of magnitude.
• The first level is when r = 0.
(a) When (m, n) = (0, 0) then The initial and boundary conditions are Practical computation of p 0,0 and other quantities of this kind is done best in a step-by-step manner. In other words, initial conditions are in general given not at the origin but at some t n ≥ 0. This presents no problem insofar as p 0,0 is concerned but a degree of care is required in computing ∂ t p 0,0 which, typically, can be computed with finite differences, consistently with the temporal order of the underlying time-stepping method. We will encounter a similar problem with other coefficients, e.g. p 3,m,n . (b) Otherwise when (m, n) = (0, 0) Note that ω plays no role in the formation of p 0,0 , which is the only differential equation that we need to compute by that stage of the algorithm.
• In the case r = 1, Let us observe, that for r = 1 and (m, n) = (0, 0): (a) when (m, n) = (0, 0) then with initial and boundary conditions 3. Constructing a numerical scheme. A common feature of eqs (3) and (4) is that they are wave equations in their own right, possessing a common form where f (t) is a source term (present only in (4)) and L(t) is a linear differential operator. Here we have deliberately suppressed the dependence on x and hidden boundary conditions for the ease of notation. Both f and L are characterised by slow variations. This slow variation in (5) allows us to use low order methods and large time steps for its solution, in contrast to (1). In particular, we will concern ourselves with order-two methods based on the application of the Strang splitting [15]. We describe how existing numerical solvers for the wave equation which work effectively for L 0 = L(τ ) for an arbitrary τ ≥ 0, can be leveraged for solving (5) up to order-two accuracy. Given initial conditions at time t = t 0 , we assume that we have access to a solver S, that provides an approximation to the solution (u(t 0 + h), u(t 0 + h) ) for (6). In the following subsection, we outline the methodS which solves (5) utilising S once. This scheme requires sampling L at the middle of the interval and f at its endpoints.
3.1.S: The extended solver. For the purpose of deriving an order-two splitting, we rewrite the equation (5) as a system of PDEs, which are order one in time.
We render the system autonomous [3,14] by adding auxiliary time variables r, s, , with the initial conditions p 0 = (u 0 , u 0 , 0, t 0 , t 0 ) . The solution of h (p 0 ). We approximate the exact flow by a Strang splitting [3,10,15], h/2 , so that the solution at time t 0 + h is approximated bỹ The flow under B is straightforward to compute exactly. For the flow under A, at t = t 0 + h, with initial conditions u(t 0 ) = u 0 and v(t 0 ) = v 0 . This system of equations is also equivalent to Assuming that the numerical solution of this wave equation is given by we may approximate [u 1 , v 1 ] by [ũ 1 ,ṽ 1 ], whereṽ 1 ≈ũ 1 − w(t a ). In other words, Composing the flows, we find that Example 2. The proposed technique is applicable to higher dimensions. To illustrate this, we consider a second example with following [7]. In practice, the operator L 0 is discretised to a matrix L 0 and the functions of this matrix are computed using Lanczos iterations. In all our examples, we utilise 25 Lanczos iterations. The proposed method A, i.e. (7), utilises the solver S twice -once directly and once in the form ofS (cf. Section 3). Table 1 compares the accuracy and costs of A with directly using the solver S for solving (1). The reference solution is computed using the sixth-order method Φ 11 [6] proposed in [2]. The number of grid points is M = 50 for all computations.
The error for A should be O(ω −3 ) for sufficiently small h. In practice, we already see this behaviour for very large time steps in Table 1 -very few calls to S suffice within the asymptotic scheme A for reaching the asymptotic accuracy (roughly 32 calls in the examples considered), while directly utilising S typically requires a very small time step to resolve high oscillations.
This highlights the low costs of the proposed scheme, which is particularly appealing in the context of large ω. For ω = 100, for instance, the method A requires two time steps (four calls to S) for reaching an accuracy of 0.1 for Example 2 (cf. Table 2). Direct utilisation of S requires 64 time steps for similar accuracy. For Example 1 under ω = 100, the accuracy achieved using two time steps of A (four calls to S) is not achieved even with 256 time steps of S (cf. Table 1).
In contrast to S, however, smaller time steps do not increase accuracy in the case of A once asymptotic accuracy is attained. Thus, once very high accuracy is required for small-to-moderate ω, the method S remains an appealing candidate. Of course, once we desire higher accuracy, we are perfectly free to use a numericalasymptotic solver along the lines of (7) that incorporates larger inverse powers of ω.
Note that the error with respect to the reference solution is well approximated as the L 1 error, where M is the number spatial grid points.
5. Comments. We have described in this paper an asymptotic-numerical approach for the solution of wave equations with external multifrequency high-oscillatory forcing. The results indicate that for large values of ω, the proposed method presents significant advantages over existing time-stepping numerical methods. A major advantage of the proposed approach is that it works equally well regardless of the frequency of oscillation (provided, paradoxically, that the oscillation is rapid enough): this represents the main advantage of methods based upon asymptotic analysis over more conventional time-stepping methods.