On asymptotic expansion solvers for highly oscillatory semi-explicit DAEs

The paper is concerned with the discretization and solution of DAEs of index $1$ and subject to a highly oscillatory forcing term. Separate asymptotic expansions in inverse powers of the oscillatory parameter are constructed to approximate the differential and algebraic variables of the DAEs. The series are truncated to enable practical implementation. Numerical experiments are provided to illustrate the effectiveness of the method.


1.
Introduction. Differential Algebraic Equations (DAEs) arise in numerous applications and in particular in circuit and device simulation [14,15]. Accurate numerical solvers or discretization of such equations present many challenges [9,10]. These challenges are compounded in the presence of highly oscillatory forcing terms and it is the purpose of this paper to address these. We wish to analyse the behaviour of the system on a time scale which is much larger than the period of the forcing term and obtain an efficient method that is not restricted by its high frequency. Considering the practical application in electronic engineering, the forcing terms on the right side are typically represented as a finite set of Fourier coefficients, such as sinusoidal AC analysis [7,8,11].
We are concerned with semi-explicit time-varying highly oscillatory DAEs of the form where η ∈ Z, x(t) : R → C d1 , y(t) : R → C d2 , ω 1, while f (x, y) : C d1 × C d2 → C d1 and g(x, y) : C d1 × C d2 → C d2 are two analytic functions. We further assume that the Jacobian ∂g/∂y is nonsingular which means that the DAEs are of index 1. In Formula (1), without loss of the generality, it is convenient to set M 1 ≤ M 2 and M 1 , M 2 ∈ N. Note that the algebraic part possesses the factor ω −1 , unlike the ODE part: this is necessary for an asymptotic expansion while reflecting many problems of practical importance.
In this work, contrary to the received wisdom in the DAE theory, we convert the problem to an ODE system, since this allows us to express its solution as an asymptotic series. We introduce two asymptotic expansions of the solutions of the DAE (1) in inverse powers of ω, In these expansions, each term in this expansion can be obtained recursively using operations which, being independent of ω, do not involve high oscillation. Besides their intrinsic value as an analytic tool, our expansions can be employed as an exceedingly affordable and precise numerical method. Regardless of whether the DAE is linear or nonlinear, our asymptotic method only takes smaller computational expense. It is a crucial observation that our approach, unlike other computational methods for DAEs, requires virtually the same computational expense regardless of the frequency ω. Recall that the range of Fourier frequencies in the forcing terms of (1) is finite, restricted to η ∈ {−M j , . . . , M j }, j = 1, 2. This will imply that the number of frequencies (i.e., the range of m) in each sum in (2)(3) is also finite (and dependent upon r, as well as M 1 and M 2 ). This renders our task somewhat simpler. The error of our asymptotic expansions is denoted by z s (t, ω) and s (t, ω), where s is a suitable truncation parameter and ω is the frequency. Once ω increases, we expect the errors z s and s to decay for fixed s: asymptotic convergence! Note that the expansion (2-3) is asymptotic and its convergence can be discussed only within this context. In other words, convergence means that for every ε > 0 there exists ω 0 > 0 such that for ω ≥ ω 0 and t ∈ [0, t * ] the expansions (2-3) are within ε from the exact solution, |z s (t, ω)| < ε and | s (t, ω)| < ε for fixed s. However, a legitimate question is subject to which conditions is z s+1 (ω), say, smaller than z s (ω) for some s ≥ 1 and whether convergence (in a classical sense) takes place for fixed t and ω as s → ∞. Like in the case with other asymptotic expansions, this is a hugely nontrivial problem and we hope to return to it in future papers.
In Section 2 we formulate the asymptotic expansion for linear DAEs with highly oscillatory terms. A numerical example is provided to illustrate the theoretical results. Linear equations, of course, are important on their own merit, but the purpose of the section is also to explore the range of phenomena and the difficulties of the asymptotic expansion (2-3) in a setting in which we can do virtually everything in an explicit and accessible manner. Section 3 is concerned with nonlinear DAEs and again it concludes with a numerical example. Finally, Section 4 explores the conclusions of our asymptotic expansion method.
2.1. Theoretical analysis. Consider a set of linear DAEs where , ω 1, and the matrix functions are The condition that D(t) is invertible is imposed so that the index of the DAE set is 1. The first equation of (4) is a differential equation. Had it been independent of y, we could have expressed its solution x as an asymptotic series [3,4]. Our contention is that, although the underlying framework is considerably more complicated, this is also the case for x in (4). Although our ansatz (2-3) predicts an infinite number of terms p r,m , the linearity of (4) implies that the only nonzero terms are p r,m for |m| ≤ M 2 . When D(t) is not singular, the solution of the algebraic equation in (4) is Thus, once we can write x in a form consistent with (2-3), we can do so also for y. This motivates us to assume (and subsequently verify) the ansatz that both x and y possess an asymptotic expansion of this form. Substituting the expression for y(t) in (5) into the differential equation in (4) yields, where We now assume that the matrices A, B, C and D (hence also E) are constant. Prior to presenting a procedure for the determination of the coefficients of the asymptotic expansions, we recall from [6] that, given a smooth continuous function for η = 0 and ω 1. Using variation of constants, the solution x has the form for r ∈ N and r ≥ 2.
We substitute the terms (6-10) into (5) for y, We note again that the number of nonzero p r,m s and q r,m s is finite, this being an artefact of linearity. Indeed, Fourier modes at one level r feed exclusively to the same Fourier modes at the level r + 1 (e.g., p r+1,m , m = 0, depends only upon p s,m for s ≤ r, as well as the coefficients a m , b m and c m ). This makes the solution of linear systems considerably easier.
For ease of reference and compassion with the nonlinear setting of Section 3, we summarize the above expansion formally, Theorem 2.1. Given the DAEs (4), if the numbers of the input oscillatory forcing terms are finite, which are from −M 1 to M 1 and −M 2 to M 2 , respectively, where M 1 ≤ M 2 , then its solution can be approximated by the asymptotic series (2-3). Furthermore, the finite number of the input terms determines the finiteness of the index m in the series of (2-3). In other words, the series (2-3) have the concrete form with the coefficients p r,m and q r,m given by (6-15).

A numerical example.
To illustrate the procedure just described, we consider a circuit as shown in Fig. 1, governed by linear DAEs where the unknown functions v(t) and i(t) represent the voltage and current variables to be solved, respectively, while h(t) is an algebraic variable and e(t) is the (known) input function. This can be rewritten in the standard form The forcing term is where ω = 2πf ω and f ω is the oscillatory parameter. Thus, a m (t) The remaining values C in , R, L are the circuit capacitance, resistance and inductance, while A inj is a constant. In line with Theorem 2.1, we expect our expansions to be of the form It is instructive to compare the absolute error at different values of ω for the asymptotic method, for different values of s ≥ 1. We compute the first few terms to show that the error tends to zero for ω → ∞. In each case, we compare the pointwise error incurred by a truncated expansion with the exact solution. Before calculating the coefficients p r,m and q r,m , we let Employing (6-10) and (11)(12)(13)(14)(15), we obtain the coefficients   Fig. 4.
It is shown in Figs 2-4 that the error decreases significantly with increasing s and ω, as expected from our expansions. Note that z 0 and z 1 for v(t) and i(t) have the same behaviour since p 1,m ≡ 0 in this example.
For the linear constant coeffients DAEs, the Maple routine dsolve without any numerical parameters can get the explicit solution. Thus, we use this routine to produce reference solutions of the DAE, to compare with our asymptotic method. With the oscillatory parameter ω increasing from 200π to 2000π, the CPU time of the traditional adaptive method dsolve for the original highly oscillatory DAE (16) requires from 0.86 to 1.516 seconds and the storage from 8.8 × 10 3 to 1.5 × 10 4 , while the asymptotic method only takes from 0.203 to 0.110 seconds and from 1.1 × 10 3 to 1.2 × 10 3 kbytes to compute the solution. The traditional method, even though it computes the exact solution, needs more CPU time and storage with an increasing oscillatory parameter. On the contrary, the computational expense of the asymptotic method does not vary because we can deduce the explicit expression of the solution for the linear DAEs. This is fully in line with our theoretical analysis.
2.3. The analysis of the factor ω −1 . In this subsection, we highlight the importance of the factor ω −1 in the algebraic part: in effect, it makes the entire expansion Let us suppose that an asymptotic expansioná la Section 2 exists for x, namely Substituting this into the algebraic part for y(t), we thus have Note the presence of the M2 η=−M2 η =0 q 0,η e iηωt term at the O(1) scale. This is fine for a linear equation but a similar state of affairs is far less manageable in a nonlinear setting. Specifically, we can no longer expand f (x, y) about the point (p 0,0 , q 0,0 ) and our entire approach is no longer suitable which we will show in Section 3.
3.1. The general theory. We now proceed to the considerably more complicated case of nonlinear highly oscillatory DAEs.
Insofar as the ordinary differential system is concerned, Condon et al. have presented in [3,4] the asymptotic expansion, Matters become more complicated for the DAE (19) because we additionally need to consider the algebraic part. To this end, we differentiate the second algebraic equation formally, Since the DAE index is one, (∂g/∂y) −1 exists. Therefore we can rewrite (21) in the form This is similar to the differential equation. Therefore, the ansatz is that the variable y(t) has an asymptotic expansion of the form Before giving these explicitly, we introduce and explain notation that shall be employed in what follows. The functions f (x, y) and g(x, y) are analytic and can be expanded in Taylor series about the point (p 0,0 , q 0,0 ) (which itself depends on t). Thus, for any p ∈ C d1 , q ∈ C d2 , both sufficiently small in norm, we expand is linear in all its arguments except for p 0,0 and q 0,0 and symmetric in each of the two groups of arguments enclosed by square brackets: specifically, it is the derivative tensor f n,k (p 0,0 , q 0,0 ) = ∂ n f (x, y) ∂x k ∂y n−k (x,y)=(p 0,0 ,q 0,0 ) .
3.4. A nonlinear example. To illustrate our approach, consider the nonlinear circuit in Fig. 5. Its behaviour is governed by the equations where I nl = S inj tanh (G n v(t)/S inj ) (cf. Fig. 5).

4.
Conclusions. Differential equations with highly oscillatory forcing terms are ubiquitous in many applications, since the forcing term corresponds to a highfrequency input into a dynamical system. Their solution by standard numerical methods is greatly restricted by the requirement that the step size should scale as the reciprocal of the largest frequency. This explains recent interest in computation which combines numerical and asymptotic insight, e.g. [1,3,4,12]. This has been focussed on ordinary and partial differential equations [13], inclusive of the multifrequency case [2], with a single paper on delay-differential equations [5]. However, many realistic problems occurring in circuit simulation require the solution of DAEs: this is among the first papers addressing the solution of DAE systems with highfrequency input using asymptotic-numerical techniques. We have demonstrated that the solution of such equations can be approximated to an exceedingly high precision using solely non-oscillatory computations -whether the solution of DAEs with no forcing terms or straightforward recursion. The oscillation is introduced into the computed solution only once the non-oscillatory ingredients are synthesised in a simple manner. Conversely (unless one is familiar with an asymptotic 'frame of mind'), the precision for fixed computational cost increases with frequency. This analysis emphasises the important point that this paper presents an initial foray into the novel and difficult area of DAEs with highly oscillatory input terms. In this paper we have presented a comprehensive analysis of one such DAE model, leading to a very effective numerical method. Yet, open questions abound, both in understanding further our approach (e.g., providing rigorous error bounds and convergence results for asymptotic expansions in a numerical setting) and extending it to other DAE models.