ASYMPTOTIC-NUMERICAL SOLVERS FOR LINEAR NEUTRAL DELAY DIFFERENTIAL EQUATIONS WITH HIGH-FREQUENCY EXTRINSIC OSCILLATIONS

. We present a method to compute efficiently and easily solutions of systems of linear neutral delay differential equations with highly oscillatory forcing terms. This method is based on asymptotic expansions in inverse powers of a perturbed oscillatory parameter. Each term of the asymptotic expansion is derived by recursion. The cost of the computation is essentially independent of the oscillatory parameter. Numerical examples are provided and show that with few terms of the asymptotic expansion, the solutions are approximated with high accuracy.


Introduction
An important and rich source of oscillatory problems with high-frequency extrinsic oscillations is computational electronic engineering.Indeed, the modeling of highly oscillating electronic circuits leads to differential equations, in particular, delay differential equations.The reason for this is that the delay occurs wherever signals are transmitted along a finite distance from one point to another.Thus, when one wishes to obtain precise communication systems, one must imperatively model the problem with delay differential equations.A wide range of applications in engineering of DDEs with highly oscillatory forcing terms can be found in [2,7,9].Among these applications, we can cite for example coupled microwave oscillators ( [3,12]), laser dynamics ( [11]) and the related secure communication techniques using chaos ( [10]).
In this paper, we are concerned with systems of linear neutral delay differential equations (NDDEs) with highly oscillatory forcing terms of the form: nature of the solution imposes a very small stepsize on standard numerical methods for solving ODEs.In effect, as has been extensively explained in [4,6,7], in any standard numerical method of order  with step ℎ, the error scales roughly like ℎ +1  (+1) ().Since the derivatives of highly oscillatory functions grow very fast, typically  (+1) () =  (︀  +1 )︀ , we are compelled to choose a very small ℎ, and therefore to require ℎ to be extremely small in order to keep the error down to an acceptable size.
The method we propose consists in solving (1.1) recursively and this on each interval [,  + 1] ,  ∈ N. Thus, on each interval [,  + 1], we get a linear ordinary differential equation of the form: Thanks to ODE solving techniques, the solution of such equation is given by: where ∆  () is the solution of the equation  ′  () =   ()+ℎ  () (which itself is the sum of the general solution of the homogeneous equation and a particular solution of the complete equation), while Ψ  (, ) is a particular solution of the equation (1.2).In other words, the term ∆  () represents the non-oscillatory part while Ψ  (, ) represents the oscillatory part of the solution of (1.2) on the interval [,  + 1].
The term ∆  () is the solution of the sequence of the following first order linear ODE's with non-oscillatory forcing term: Each equation of this sequence, can be solved numerically by standard numerical solvers of first order linear ODE's, see [8].Note also that the non-oscillatory term of the solution can also be seen as the solution of the neutral delay differential equation: This type of equations can be solved thanks to the numerical methods specific to the differential equations with delay, see [1].
The second term Ψ  (, ), representing the oscillatory part of the solution of (1.2) on the interval [,  + 1], corresponds to the solution of the sequence of the following first order linear ODE's endowed with highly oscillatory forcing term: Recall that the resolution of the equation (1.1) has been studied in [6], in the following restricted framework, The main result of this paper is to resolve (1.5).The approach we take is to divide the whole interval into the subintervals [,  + 1].On each subinterval [,  + 1], the oscillatory term of the solution is obtained by a succession of integrations by parts of (1.5).It is written by distinguishing two cases: the case where (1.1) is a scalar equation and the case where (1.1) is a system of equations.Indeed, we show that Ψ  (, ) is of the form: ) is a scalar differential equation, and of the form 1) is a system of equations of dimension  ≥ 2, and we give the exact expression of the terms Ω , (, ) and Φ , (, ) appearing in (1.7) and (1.8).
The current approach allows us to determine the coefficients Ω , (, ) and Φ , (, ) of (1.7), Ω , (, ) of (1.8) recursively and to determine the coefficient Φ , (, ) of (1.8) by solving a single non-oscillatory ODEs.This represents a great advantage in terms of computational cost compared to conventional solving methods.Indeed, unlike the classical methods, the current approach is completely independent of the size of .On the contrary, it becomes more precise when the frequency  is increased since once (1.7) (resp.(1.8)) is truncated for  ≥ 1, we obtain a numerical approximation whose error,  (︀ 1/ +1 )︀ , actually improves for growing frequency.Moreover, we show in the scalar case, when the function () appearing in (1.1) is a polynomial of degree , the exact solution is obtained in the form of a finite sum.More precisely, we obtain Ψ  (, ) where the Φ , () and the Ω , (, ) are given recursively.Finally, let us recall that the idea which consists in not using the classical methods of resolution of NDDEs, but rather to write the solution in the asymptotic form (, ) and then, determining the   (, ), has been widely and successfully used in solving some differential equations with high-frequency extrinsic oscillations, see [4,5].The paper is organized as follows: In Section 2, we justify that the Ψ  (, ) actually have the forms (1.7) and (1.8).In Section 3, we deal with the case where (1.1) is a scalar differential equation.We give all the terms of the asymptotic development (1.7).We show that these terms have a very simple expression on [0, 1] and that on the other intervals, these terms are obtained recursively without having to solve any differential equation.In Section 4, we treat in a similar way the case where (1.1) is a system of differential equations, then we study the stability of the proposed algorithm.At the end of each of Sections 3 and 4, we give several numerical examples, computed by MATLAB, to show the efficiency of the proposed algorithms.
Remark 1.If we have to resolve an equation with several oscillatory source terms, i.e., equation of the form the oscillatory term of the solution, thanks to the superposition of solutions, is in the scalar case of the form: and in the vectorial case of the form Remark 2. In order to simplify the writing, we omit to write the second parameter , in Ψ , (, ) , Ω , (, ) and Φ  (, ) .

General setting
As it was announced in the introduction, we write the solution of equation (1.1) in the form: where ∆  () (respectively Ψ  ()) represents the non-oscillatory part (respectively the oscillatory part) of the solution on [,  + 1].The oscillatory part of the solution, verifies for  = 0, the first order linear ODEs: and for  ≥ 1, the perturbed linear ODEs: 2) The solution of equation (2.1) is given by: This result is obtained by multiplying the two members of equality (2.1), by  − and by integrating between 0 and  ∈ [0, 1].
A simple integration by parts of (2.3), gives us in the scalar case: and gives us in the vectorial case where Ω 0, () and Φ 0, () can be obtained recursively.
For  ≥ 1, by multiplying the two members of equality (2.2) by  − and then by integrating between  and  ∈ [,  + 1] , we get since Ψ  () = Ψ −1 (): (2.6) From (2.6), a reasoning by induction on , suggests to us that in the scalar case, Ψ  () is of the form and in the vectorial case, Ψ  () is of the form (2.8)

Numerical examples
Let us consider the following example: With the help of MATLAB software, the exact oscillatory part of the solution of this equation is determined by (2.3) and (2.6).Thus, we are able to compare it with the approximate oscillatory term of the solution proposed in this section.See Figures 1a and 1b.As it is seen from these figures, the asymptotic error decreases for increasing .Furthermore, the accuracy of the asymptotic method increases greatly for the same number of  levels for higher values of .Let us now consider the second example: and let us compare the results given by the current algorithm and the MATLAB routine ddensd, whether in terms of accuracy or in terms of execution time.See Figures 3a and 3b.As seen in Figure 3b, the MATLAB routine ddensd gives a worse and worse approximation as one moves away from 0. Also, the execution time becomes longer and longer when  becomes very large.On the other hand, we note on Figure 3a, an undeniable superiority of the current algorithm compared to the MATLAB routine ddensd, whether on the accuracy or on the execution time.This feature makes the current algorithm most suitable for simulation of linear neutral delay differential equations with high-frequency extrinsic oscillations.

Numerical examples
Let us consider the system modelled by the following second order linear ODE: )︂ .
The real parts of both eigenvalues of the matrix  are negative.Thus, we have asymptotic stability according to Theorem 1.
For this particular example, thanks to (2.3) and (2.6), we can compute exactly the Ψ  () with MATLAB software and compare it to the approximate solution proposed in this section.See Figures 4a-5b.As in the scalar case, we observe from these figures, the asymptotic error decreases for increasing  and the accuracy of the asymptotic method increases greatly for the same number of  levels for higher values of .5b): The minus of decimal logarithm of absolute values of the errors in the oscillatory term of the solution of (4.15) in x(t) (resp.x'(t)) on the interval [0,5], using the same equations and parameters used to obtain Figure 4a (resp.Fig. 4b), but this once for  = 1000.

Conclusion
We have shown that the solution of the equation  ′ () = ()+(−1)+ ′ (−1)+ℎ()+  ∑︀ =1   ()  , can be written as the sum of two terms.The first term represents the non-oscillatory part of the solution and is solution of a certain ODE which can be resolved by one of the classical methods of resolution of ODEs.The second term is the oscillatory part of the solution that we have expanded into asymptotic series in inverse powers of the frequencies   .We have shown that each term of the asymptotic series is obtained recursively.We have seen that with few terms of this asymptotic expansion, the oscillatory part of the solution is approximated with highly accuracy and at a lower cost.

Figure 1 .
Figure 1.(resp.Fig. 1b): The minus of decimal logarithm of absolute values of the errors in the oscillatory term of the solution of (3.21) on the interval [0,10] with  = 100, (resp. = 1000) using the equations (3.1) and (3.2) on [0,1], and the equations (3.11)-(3.14) on [1,10] with  = 1−12 (resp. = 1−8) from bottom row to top row.If we subtract the computation time of the exact solution given by (2.3) and (2.6), the execution time of the twelve curves in Figure 1a was 40.66 seconds while that of the eight curves of Figure 1b was 27.6 seconds.

Figure 2 .Remark 5 .
Figure 2. (resp.Fig. 2b): The minus of decimal logarithm of absolute values of the errors in the oscillatory term of the solution of (3.21) on the interval [0,10], using the same equations and parameters used to obtain Figures 1a and 1b, but this time taking into account Remark 4.

Figure 3 .
Figure 3. (resp.Fig. 3b): The minus of decimal logarithm of absolute values of the errors in the oscillatory term of the solution of (3.22) on the interval [0,10] with  = 200, using our algorithm, with  = 1−12 from bottom row to top row (resp.the Matlab routine ddensd).If we subtract the computation time of the exact solution given by (2.3) and (2.6), the execution time of the twelve curves in Figure 3a was 40.18 seconds while that of the curve of Figure 3b was 278.36 seconds.

Figure 5 .
Figure 5. (resp.Fig.5b):The minus of decimal logarithm of absolute values of the errors in the oscillatory term of the solution of (4.15) in x(t) (resp.x'(t)) on the interval [0,5], using the same equations and parameters used to obtain Figure4a(resp.Fig.4b), but this once for  = 1000.