## Abstract

### Purpose

The paper proposes an efficient and insightful approach for solving neutral delay differential equations (NDDE) with high-frequency inputs. This paper aims to overcome the need to use a very small time step when high frequencies are present. High-frequency signals abound in communication circuits when modulated signals are involved.

### Design/methodology/approach

The method involves an asymptotic expansion of the solution and each term in the expansion can be determined either from NDDE without oscillatory inputs or recursive equations. Such an approach leads to an efficient algorithm with a performance that improves as the input frequency increases.

### Findings

An example shall indicate the salient features of the method. Its improved performance shall be shown when the input frequency increases. The example is chosen as it is similar to that in literature concerned with partial element equivalent circuit (PEEC) circuits (Bellen *et al.*, 1999). Its structure shall also be shown to enable insights into the behaviour of the system governed by the differential equation.

### Originality/value

The method is novel in its application to NDDE as arises in engineering applications such as those involving PEEC circuits. In addition, the focus of the method is on a technique suitable for high-frequency signals.

## Keywords

## Citation

Condon, M. (2024), "Numerical analysis of neutral delay differential equations with high-frequency inputs", *COMPEL - The international journal for computation and mathematics in electrical and electronic engineering*, Vol. 43 No. 1, pp. 14-23. https://doi.org/10.1108/COMPEL-12-2022-0423

## Publisher

:Emerald Publishing Limited

Copyright © 2023, Marissa Condon.

## License

Published by Emerald Publishing Limited. This article is published under the Creative Commons Attribution (CC BY 4.0) licence. Anyone may reproduce, distribute, translate and create derivative works of this article (for both commercial & non-commercial purposes), subject to full attribution to the original publication and authors. The full terms of this licence may be seen at http://creativecommons.org/licences/by/4.0/legalcode

## 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:

*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:

*A*,

*B*and

*C*are

*n*×

*n*constant matrices.

*T*is the time of simulation or time range of interest.

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

*b*

_{1,}

*(*

_{m}*t*)e

^{i}

*is differentiated, the result has two terms. In one of these,*

^{ωmt}*b*

_{1,}

*(*

_{m}*t*)i

*mω*e

^{i}

*, the*

^{mωt}*ω*cancels with the

*r =*0 layer.

Consider *r =* 0. One collects the
*ω*:

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 *r*th 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 *r*th 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*−

*C*e

^{−}

*) are less than 1. It follows then that*

^{imωτ}*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

*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*−

*C*e

^{−}

*) are less than 1.*

^{imωτ}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*is the reference model and

_{ref}*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*_{0}(*t*). If the error associated with the *p*_{0}(*t*) is removed by using the *accurate p*_{0}(*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
*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*_{0}(*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^{−7}s and the time step used for the two methods is 1 × 10^{−3}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*)^{−1}(*A* +*B*)] < 0 and if *F*(*y*, *z*) and *G*(*y*, *z*) in equation (4.5) have no common real roots
*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*(0) = 0, and so

_{m}*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*_{0}(*t*) term. Figure 3 shows a plot of the *p*_{0}(*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*

_{0}(

*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.

## Figures

A comparison of the root mean square errors in the various methods

ω |
Finite difference method | Asymptotic model | Asymptotic model with accurate p_{0}(t) |
---|---|---|---|

20 | 2.62 × 10^{–4} |
0.0025 | 0.0025 |

100 | 2.95 × 10^{–4} |
1.54 × 10^{–4} |
2.00 × 10^{–5} |

200 | 2.98 × 10^{–4} |
1.53 × 10^{–4} |
2.72 × 10^{–6} |

400 | 3.02 × 10^{–4} |
1.53 × 10^{–4} |
3.93 × 10^{–7} |

**Source:** Author’s own work

## References

Achar, R. (2011), “Modeling of high-speed interconnects for signal integrity analysis: part I”, IEEE Microwave Magazine, Vol. 12 No. 5, pp. 61-74.

Ait el Bhira, H., Kzaz, M. and Maach, F. (2022), “Asymptotic-numerical solvers for linear neutral delay differential equations with high-frequency extrinsic oscillations”, ESAIM (to appear).

Bellen, A., Gugielmi, N. and Ruehli, A. (1999), “Methods for linear systems of circuit delay differential equations of neutral type”, IEEE Transactions on Circuits and Systems I, Vol. 46 No. 1, pp. 212-216.

Chao-Jung Cheng, J.J.Y., Liao, T.L. and Hwang, C.C. (2006), “Global asymptotic stability of a class of neutral-type neural networks with delays”, IEEE Transactions on Systems, Man and Cybernetics, Vol. 36 No. 5, pp. 1191-1195.

Condon, M., Deao, A., Iserles, A. and Kropielnicka, K. (2012), “Efficient computation of delay differential equations with highly oscillatory terms”, ESAIM. Mathematical Modelling and Numerical Analysis, Vol. 46 No. 6, pp. 1407-1420.

Ferranti, F., Nakhla, M., Antonini, G., Dhaene, T., Knockaert, L. and Ruehli, A. (2011), “Multipoint full-wave model order reduction for delayed PEEC models with large delays”, IEEE Transactions on Electromagnetic Compatibility, Vol. 53 No. 4, pp. 959-967.

Ferranti, F., Romano, D., Antonini, G. and De Camillis, L. (2017), “Stochastic collocation for uncertainty quantification of systems described by neutral delayed differential equations”, 2017 International Applied Computational Electromagnetics Society Symposium – Italy (ACES), Firenze, Italy, 2017, pp. 1-2.

Gad, E., Tao, Y. and Nakhla, M. (2022), “Fast and stable circuit simulation via interpolation-supported numerical inversion of the Laplace transform”, IEEE Transactions on Components, Packaging and Manufacturing Technology, Vol. 12 No. 1, pp. 121-130.

Gentile, F., Itovich, G. and Moiola, J. (2023), “Stability analysis of some neutral delay-differential equations with a frequency-domain approach”, Discrete and Continuous Dynamical Systems - Series B, Vol. 28 No. 3, pp. 1787-1805.

Hui, F. and Jibin, L. (2001), “On the existence of periodic solutions of a neutral delay model of single-species population growth”, Mathematical Analysis and Applications, Vol. 259 No. 1, pp. 8-17.

Kyrychko, Y. and Hogan, S. (2010), “On the use of delay differential equations in engineering applications”, Journal of Vibration and Control, Vol. 16 Nos 7/8, pp. 943-960.

Liao, X. (2016), “Dynamical behavior of Chua’s circuit with lossless transmission line”, IEEE Transactions on Circuits and Systems I: Regular Papers, Vol. 63 No. 2, pp. 245-255.

Liu, M., Dassios, I. and Milano, F. (2017), ‘“Small-signal stability analysis of neutral delay differential equations”, *IECON 2017, 43rd Annual Conference of the IEEE Industrial Electronics Society*, IEEE, pp. 5564-5649.

Liu, M., Dassios, I. and Milano, F. (2019), “On the stability analysis of systems of neutral delay differential equations”, Circuits, Systems, and Signal Processing, Vol. 38 No. 4, pp. 1639-1653.

Maleki, M. and Davari, A. (2021), “Analysis of an adaptive collocation solution for retarded and neutral delay systems”, Numerical Algorithms, Vol. 88 No. 1, pp. 67-91.

Mohammad, M. and Trounev, A. (2022), “A new technique for solving neutral delay differential equations based on Euler wavelets”, Complexity (New York, N.Y), Vol. 2022, pp. 1-8.

Wang, W., Zhang, Y. and Li, S. (2009), “Stability of continuous RungeKutta-type methods for nonlinear neutral delay-differential equations”, Applied Mathematical Modelling, Vol. 33 No. 8, pp. 3319-3329.

Wei, P., Guan, Q., Yu, W. and Wang, L. (2008), “Easily testable necessary and sufficient algebraic criteria for delay-independent stability of a class of neutral differential systems”, Systems and Control Letters, Vol. 57 No. 2, pp. 165-174.