Numerical analysis of neutral delay differential equations with high-frequency inputs

1. Introduction

The paper addresses the simulation of neutral delay differential equations (NDDE) subject to a high-frequency input, as addressed in Ait el bhira et al. (2022). NDDE are important in many engineering applications, for example, partial element equivalent circuit (PEEC) circuits in electromagnetic compatibility analysis (Ferranti et al., 2011, 2017), coupled oscillatory systems with non-instant connections (Kyrychko and Hogan, 2010), control problems (Gentile et al., 2023), neutral delayed neural networks (Chao-Jung Cheng et al., 2006), population growth models (Hui and Jibin, 2001), circuit analysis (Liao, 2016) and many more. High-frequency inputs arise, for example, in modulated signals and high-frequency noise.

Several methods exist for the numerical solution of NNDE. Bellen et al. (1999) examine the use of continuous Runge–Kutta methods with linear interpolation. They give a theorem stating that a Runge–Kutta method of uniform global error p used in conjunction with a direct-evaluation scheme has a uniform global order p. They also indicate that the Lobatto III-C method with linear interpolation is suitable and give conditions concerning its stability. Wang et al. (2009) show that Runge–Kutta methods based on the backward Euler method and used in conjunction with linear interpolation are contractive and asymptotically stable. Maleki and Davari (2021) present adaptive collocation methods and examine the convergence properties of the method. Recently, Mohammad and Trounev (2022) proposed the use of Euler wavelets for numerical solution of NNDEs. The authors of Ferranti et al. (2017) use stochastic collocation to assess parameter variability in PEEC circuits.

The focus of this contribution is on a method appropriate when high-frequency signals are present. High-frequency inputs necessitate the use of very small time steps in traditional solvers, and hence, the solution of such equations is computationally expensive. To address this, the paper shall first describe the proposed numerical method that involves an asymptotic expansion. An asymptotic expansion was used for ordinary delay differential equations (not NDDE) in the work by the authors of Condon et al. (2012). The paper shall consider NDDE and then examine the stability of the method. The value of the proposed method is that its accuracy increases with increasing input frequency. This is important in modern engineering applications where frequencies of operation are ever increasing, for example, Gad et al. (2022) and Achar (2011). The method also structures the solution such that important aspects and properties of the system considered can be readily observed.

2. Asymptotic analysis

When equations governing, for example, PEEC circuits are written taking into account retarded mutual coupling and inclusive of non-linear lumped elements, the form is as follows:

(x(t),u(ω,t)):ℝ→ℝn, ω>>1. In general, f is a non-linear function. However, in many applications, the nonlinearity is weak, and this shall be assumed in the present work. As in works such as Liu et al. (2017), small signal analysis is the primary concern, and as such a linear approximation about an equilibrium point is deemed sufficient. In addition, to make it easier to link this work to that by Ait el bhira et al. (2022) and Condon et al. (2012), the linearised equation shall be written as:

Let:

Henceforth, the summation for m shall be from −∞ to ∞ unless otherwise indicated.

The analysis proceeds by substituting the expansion for x(t) into equation (2.2):

Two kinds of scales are identified, powers of ω and frequency terms of the form e^iωmt for distinct m. Firstly, the powers of ω are separated and then the frequencies. This results in delay differential equations and recursive relations for each magnitude scale r. The r = 0 layer corresponds to the case 1ω0. Note that when b_1,m(t)e^iωmt is differentiated, the result has two terms. In one of these, b_1,m(t)imωe^imωt, the ω cancels with the 1ω resulting in a term in the r = 0 layer.

Consider r = 0. One collects the O(1) terms in ω:

Next one separates frequencies.

For m = 0:

For m ≠ 0:

Now proceed to the r = 1 layer:

For m = 0:

For m ≠ 0:

Proceeding on to the next layer:

For the rth layer:

The general expressions for the equations for b_r,m(t) are as follows:

It is important to note that there is no high-frequency input in the NDDE for each layer r.

3. Stability

Each r level involves a delay differential equation and recursive relations. Note that the form of the delay differential equation in each layer is identical. The stability of the complete numerical method is governed by the stability of each r. Following Liu et al. (2019) and Wei et al. (2008), the delay differential equation (2.2) shall be asymptotically stable if the roots of the characteristic equation have negative real parts and marginally stable if the roots have a zero real part. The characteristic equation is:

The delay differential equation for the rth layer, equation (2.13), is of the same form and so the same conditions for stability apply.

Now, consider equation (2.7). If a_m(t) is bounded, the recursive equation for b_1,m(t) shall be stable if the absolute value of the roots of det(λI − Ce^−imωτ) are less than 1. It follows then that b_1,m(t) shall be bounded. Now, b_1,m(t) ≡ 0 for t < 0. So if a_m(t) = 0 at t = 0, then b_1,m(t) shall be continuous at t = 0. However, if a_m(t) or higher-order derivatives of a_m(t) are discontinuous at t = 0, the same derivatives of b_1,m(t) shall be discontinuous at t = 0. The first derivative of b_1,m(t − τ) will have a jump at τ unless b_1,m(t) is at least C1 at 0. This impacts the recursive equation for b_2,m(t), equation (2.10). The behaviour of b_2,m(t) at t = τ depends on the behaviour of b_1,m(t). In particular, b_2,m(τ) depends on the degree of smoothness of b_1,m(t) at t = 0. Note also that the equation shall again be stable if the absolute value of the roots of det(λI − Ce^−imωτ) are less than 1.

A similar analysis continues for r > 2.

4. Example

The example that is considered, similar to that in Bellen et al. (1999) where it is representative of a small-scale PEEC circuit, is as follows:

The input is:

The first step is to determine a reference model. A finite difference solution is obtained with a very small time step when there is no oscillatory input to form a reference model for the case when there is no oscillatory input. The finite difference method that is used is as follows:

The equation is then solved with a larger time step and the root mean square error is computed as follows:

N_t is the number of samples in the timespan of the simulation. y_ref is the reference model and y is the result with the larger time step.

The larger time step is accepted if the error is at an acceptable level. This is the required step length for delay differential equations at each level r that have no input. A reference solution is then determined for each value of ω. Using the same larger time step as selected earlier, the finite difference method with the oscillatory input and the asymptotic model with r = 2 stages are implemented and the error is computed. Table 1 shows the error associated with the finite difference method. Table 1 also shows the error associated with the asymptotic method. It can be observed that the error initially falls with increasing ω. However, the error then stagnates. The reason for this is that the significant error arises from the time step chosen for the determination of p₀(t). If the error associated with the p₀(t) is removed by using the accurate p₀(t) term, the increase in accuracy of the asymptotic method with frequency is clearly more evident as seen in the third column in Table 1. In addition, note that the error reduces as O(ω(−3)) as expected with r = 2 terms [see equation (2.3)] when there is no error in the numerical method for the delay differential equations in each layer r. In general, the p₀(t) could be pre-computed, as this does not change with a change in input frequency. For this example, the time of simulation is 10 s, the reference time step is 5 × 10⁻⁷s and the time step used for the two methods is 1 × 10⁻³s. Note, any method suitable for solving NDDE could be used for the differential equations in each layer in the asymptotic method, as there is no high-frequency input for these equations. The important point is that the accuracy of the proposed method increases with increasing frequency without needing a corresponding reduction in the time step size when the frequency of the input increases.

From a stability viewpoint, some of the right-most zeros of the determinant in equation (3.1) are −0.6918 ± 2.1888i,−1.9632 ± 4.8762i, –1.3853 and –2.0247. As these have negative real parts, they indicate that the equation is asymptotically stable.

Stability can also be investigated using the method presented in Wei et al. (2008). Consider:

Let f(z) and g(z) be the real and imaginary parts of d(z).

Consider:

Let F(y, z) and G(y, z) be the real and imaginary parts of D(y, z).

Then if f(z) and g(z) (equation 4.4) have no common real roots, if Reλ[(I − C)⁻¹(A +B)] < 0 and if F(y, z) and G(y, z) in equation (4.5) have no common real roots (y,z),z∈ℝ,y∈ℝ−{0}, the differential equation (2.2) is delay-independent stable [Theorem 2 (Wei et al., 2008)].

For the given system, all three of these conditions hold true confirming the stability of the equation for each level r.

The magnitudes of the eigenvalues of C are 0.0195, 0.0781 and 0.0986, again indicating that the difference equation (2.7) is asymptotically stable once a_m(t) is bounded. a_m(0) = 0, and so b_1,m(t) are continuous at t = 0. However, b_1,m’(t) are discontinuous at t = 0. This results in discontinuous b_2,m(t), see equation (2.10). Figures 1 and 2 show a plot of |b_1,1(t)| and |b_2,1(t)| to illustrate this.

A second feature of the method is the structure of the solution. In the proposed method, the solution with a constant input corresponds to the p₀(t) term. Figure 3 shows a plot of the p₀(t) term, and its oscillation frequency gives an indication of the dominant natural frequencies of the system. For example, an estimate of the dominant frequency is 2.2. This is in line with the frequency of the right-most zero, which is a complex number −0.6918 ± 2.1888i. As regards p_r,m(t), r > 0, these account for the input oscillation terms. These lie on the p₀(t) envelope. As the frequency increases, the number of r layers required for a particular level of accuracy reduces. Hence, the number of r layers can be selected based on the accuracy and efficiency requirements.

5. Conclusion

The paper addresses the simulation of NDDE subject to a high-frequency input. An asymptotic method is presented and its stability has been studied. An example is given to illustrate the salient features of the method. For a fixed level of computation, its accuracy increases with increasing frequency. The differential equations in each layer no longer have a time step governed by a high-frequency input. Furthermore, the structure of the asymptotic expansion is such that its form gives insights into the intrinsic behaviour of the system and the response of the system to high-frequency inputs. Possible future work includes extension to non-linear equations and exploration of the effect of the strength of the nonlinearity on accuracy, efficiency and its impact on stability.