A spectral method for fractional nonlinear transmission lines

Nonlinear transmission lines (NLTL) have many applications for modelling wave propagation in nonlinear media (Aziz et al., 2020; El-Borai et al., 2017) and in high-frequency circuits. Examples include electrical oscillators (Yildirim, Ricketts and Ham, 2009), pulse compressors in radar systems and pulse generators for high-speed sampling systems (Nouri et al., 2017) and for cancer treatment (Abe and Minamitani, 2012). NLTLs can be formed from lumped elements or by placing a nonlinear capacitor periodically along a distributed transmission line (Chandra Paul et al., 2023). Many research papers have been written examining the behaviour of lossless and lossy NLTLs, for example Sekulic et al. (2012); Nikoo and Hashemi, (2017) and those listed in Chandra Paul et al. (2023). However, at high frequencies, fractional differential equations are considered as more realistic models for capacitors and inductors (Antonini et al., 2021) and some papers have been written on fractional transmission lines, for example R-Smith, Kartci and Brancik, (2017) and Gómez-Aguilar and Baleanu, (2014). Recently, the authors in Chandra Paul et al. (2023) examined the behaviour of solitons on the fractional nonlinear transmission line when beta derivatives are considered. They use the $(G′/G2)$ method. Experimental results showing the increased accuracy of a fractional-order transmission-line model are given in Shang et al. (2013).

Many approaches have been proposed for the numerical simulation of fractional ordinary and partial differential equations, for example Papadopoulas and Olver (Papadopoulos and Olver, 2024) examine the development of a sparse spectral method using weighted Chebyhev polynomials of the second kind and their Hilbert transforms. The authors in [Sweilam et al. (2016)] propose a sinc-Legendre collocation method. Several linear multistep methods and product integration methods are given in (Garrappa, 2018). The present paper examines the behaviour of a lossy nonlinear transmission line with fractional time derivatives. Such a study is important for examination of the effects of fractional derivatives on wave propagation in nonlinear media. It proposes using a spectral method that uses the Malmquist-Takenaka functions. Some recent works on the use of such functions for spectral methods are (Iserles and Webb, 2020; Shen and Wang, 2022). The present work proposes the use of the Malmquist-Takenaka functions in conjunction with the fractional and integer implicit trapezoidal rule. The evolution of a soliton and a wave packet shall be examined and the effect of increasing spatial frequency on the evolution of a wave packet shall be illustrated. Wave packet dynamics are of interest in several areas of physics (Karjanto, 2024). For problems on an infinite line, Hermite functions could be used as an orthogonal basis. However, while Hermite functions are suitable for solitons, the Malmquist-Takenaka functions are superior for wave packets. The Malmquist-Takenaka functions are orthonormal and complete in $L2(ℝ)$ and like the Hermite functions, result in a differentiation matrix corresponding to the first derivative that is tridiagonal. The differentiation matrix is banded for higher derivatives. However, in contrast to the Hermite functions, the expansion coefficients can be computed in an efficient manner using the fast fourier transform (FFT) and in addition, a formula for the product of expansions exists which is very useful for nonlinear partial differential equations. Furthermore, when approximating wave packets, the coefficients exhibit asymptotic exponential decay and the number of coefficients required grows linearly with frequency in contrast to the quadratic growth rate for Hermite coefficients (Iserles et al., 2023).

Consider the NLTL shown in Figure 1.

The nonlinear capacitance is governed by the following equation:

U_n is the voltage at node n, c₀ is the linear capacitance, α and β are constants. Q_n is the charge stored in the nth capacitor:

Using Kirchhoff’s Laws:

and:

r and l are the resistance and inductance of section n, respectively. I_n is the current flowing in section n. Rearranging equation (2.4) gives:

Differentiating equation (2.3) gives:

If equation (2.5) is substituted in equation (2.6), then:

Substitution of the expression for Q_n retaining only the first two terms on the assumption of a small deviation from linearity yields:

The right-hand side of equation (2.8) may be approximated as a partial derivative with respect to x. Let $Un(t)=u(x,t)$ and $Un±1(t)=u(x±δ,t)$ where δ be the length of a section. Using Taylor expansions:

So:

Terms of order $δ4$ and higher are assumed negligible and additionally, it is assumed that the time variations of the voltage are small. Hence, as shown in the Appendix, the following is the resulting partial differential equation for the voltage, u(x, t):

where $a0=1LC0,k0=−RL$ and $b0=−2αRL$⁠. $R,L,C0$ are the per-unit length parameters, i.e. $R=r/δ,C0=c0/δ,L=l/δ$⁠. This result is similar to that given in Chandra Paul et al. (2023) (but with $r2=0$⁠).

The initial and boundary conditions are:

Equation (2.10) can be rewritten as:

Spectral methods involving orthonormal bases tend to result in rapid convergence (Gao and Iserles, 2023). However, spectral methods can result in ill-conditioned dense matrices (Olver and Townsend, 2013) and care must be taken to ensure that the method is stable. An appropriate set must be chosen for the problem under consideration. For an infinite line, the basis should be orthogonal on the real line. Hermite functions and Malmquist-Takenaka functions are possible bases. Malmquist-Takenaka functions are suitable for approximating functions that have poles in the complex plane (Shen and Wang, 2022) and hence suitable for the type of functions that arise in engineering applications – functions involving decaying exponentials. Furthermore, examples in Iserles and Webb, (2020) have shown their ability to approximate certain functions is superior to Hermite functions and in particular wave-packet functions as shown in a detailed study in Iserles et al. (2023).

The Malmquist-Takenaka functions are defined as:

The system $φn(x)n∈ℤ$ forms a complete and orthonormal basis in $L2(ℝ)$ i.e.:

Expansions in Malmquist-Takenaka functions are presented as:

The expansions coefficients are determined as follows:

where $φn(x)¯$ is the complex conjugate of $φn(x)$⁠. Letting $x=tan⁡(θ/2)/2$⁠:

This is a Fourier integral and the benefit of this is that the coefficients may be computed in an efficient manner using the Fast Fourier Transform (Iserles and Webb, 2020).

The differentiation matrix for an orthonormal system is defined as:

and in general, the elements of the differentiation matrix are formed as:

However, the Malmquist-Takenaka functions satisfy:

Hence, the differentiation matrix is tridiagonal and skew hermitian and can be formed without the necessity for integration. Similarly, the second derivative may be determined from:

Therefore:

This matrix is banded with a bandwidth of 2 and is hermitian. This leads to computational efficiencies when solving equation (2.11).

Many forms of fractional derivative exist, for example, the Caputo, the Riemann-Liouville and Gr$u¨$nwald-Letnikov (de Oliveira and Machado, 2014). However, the Caputo fractional derivative is selected as only integer-order initial conditions are required (Garrappa et al., 2019):

$m−1<μ≤m$ and m is a positive integer. f(t) is a continuous function and $fm(t)$ is the mth derivative of f(t).

If fractional time derivatives are considered for the second equation in 2.11, the equation set becomes:

Once u(x, t) and v(x, t) are expanded in the Malmquist-Takenaka basis, a numerical method for determination of the time evolution of the expansion coefficients is required. For the present work, the trapezoidal rule is used for integer derivatives and the product-integration implicit trapezoidal rule shall be selected (Garrappa, 2018) for fractional derivatives. However, splitting methods such as described in Cao et al. (2015) could also be explored if the nonlinearity coefficient $b0<<1$⁠. The trapezoidal rule has an order of convergence of $O(h2)$ while the product-integration trapezoidal rule for fractional differential equations has an order of convergence of $O(hmin(1+μ,2))$ (Garrappa, 2018). The product-integration rule (Garrappa, 2018) is summarised as follows:

Consider a fractional ordinary differential equation:

$f(t,y(t))$ is assumed continuous and $y01…y0m−1$ are the $1,…,m$ derivatives of y(t) at t₀:

The product-integration rule is based on reformulating equation (4.3) as a Volterra integral equation:

The time interval $[t0,t]$ is divided into a grid $tn=t0+nh$ with a constant stepsize h. Equation (4.5) is then rewritten as:

$f(τ,y(τ))$ is approximated as:

With this approximation in equation (4.6), equation (4.4) follows.

Let:

N is the number of orthogonal functions. Let $u^=[u^−N/2,..,u^N/2−1]$ and $v^=[u^−N/2,..,v^N/2−1]$ be the vectors of the coefficients $u^n$ and $v^n$⁠.

The following property holds for the Malmquist-Takenaka functions that makes them suitable for nonlinear partial differential equations involving products of expansions:

Hence:

Equation (2.11) becomes:

where $D2$ is expressed in equation (3.9):

The banded nature of D₂ and the use of the property in equation (5.2) results in computational savings.

For comparative purposes, the method shall also be implemented using the Hermite functions. The Hermite functions (Iserles et al., 2023) are defined as:

There is however no property similar to 5.2 for products of Hermite function expansions and therefore, it is necessary to transform back and forth between the spectral space and the physical space to implement the nonlinear part of the equation (Trif and Petrila, 2009). A similar approach to that described in Trif and Petrila, (2009) is implemented for the spatial discretisation and the product integration method used for the fractional time derivatives:

where $cu(t)$ and $cv(t)$ are the columns of coefficients and T is the matrix of values of $φnh(x)$ evaluated at the roots of the Hermite polynomials. Two recurrence relations for Hermite polynomials are used for determining the spatial derivative matrix:

From these, a banded $D2$ matrix can be formed. Equation (2.11) becomes:

* is the multiplication of the individual elements of the two vectors (dot multiplication).

To show the efficacy of the proposed method, the partial differential equation with μ = 1 shall be solved using the leapfrog method (Iserles, 2009) implemented with a very small timestep:

$r=δt/δx.$ Subscript l denotes the spatial step index. Superscript n refers to the time index. The error is defined as:

where u_ex is the result obtained using the leapfrog method.

A similar finite-difference method shall be implemented for the case $μ<1$⁠. The spatial derivative is implemented as a second-order finite difference and the time derivatives shall be implemented with the product-integration trapezoidal rule.

Consider equation (2.11) when $a0=1,b0=−0.1$ and $k0=−1$⁠. The initial condition is $u(x,0)=sech2(x12)$⁠. The simulation is performed for a duration of $T=1s$⁠. The timestep is $h=1×10−3s$⁠. The choice of timestep and the number of orthogonal functions is selected so as to meet a specified tolerance. Table 1 shows the fall in the error as the number of orthogonal functions is increased for both the Malmquist-Takenaka functions and the Hermite functions. As is evident from the results in Table 1, Hermite functions are more efficient for representing solitons. However, this is not the case for wave packets.

Consider the case when the initial condition is a wavepacket $e−x2cos⁡(ωx)$ (Iserles and Webb, 2020). The spatial frequency ω shall be varied. As the frequency increases, the number of Hermite basis functions required to represent a wave packet increases quadratically. Issues with ill-conditioned matrices arise when determing a numerical solution to the differential equation. Table 2 shows the increase in the error as ω increases when n = 64 basis functions are included. This is expected as the roots of the Hermite functions are insufficiently close together to capture the increasing spatial frequency.

Noting the results in Table 2, the Malmquist-Takenaka functions are selected and the number of expansion coefficients is chosen so that the error is less than a specified tolerance when compared to an implementation of the finite-difference method with a very small time step. The number of orthogonal functions included in the simulation is $N=26$⁠. Figure 2 shows the results at $T=1s$ when the initial condition is the soliton. The results for the integer time derivative and the fractional time derivative with $μ=0.75$ are shown. As evident from the results, the fractional derivative has little impact on the evolution of the soliton. Figure 3 shows the same results obtained with the finite difference method to confirm the validity of the results.

Consider the case when the initial condition is a wavepacket $e−x2cos⁡(ωx)$⁠. $N=28$ basis functions are used for the following simulations. Figure 4 shows the magnitude of the coefficients of the Malmquist-Takenaka functions $u^(T)$ with $T=0.1s$ and when $ω=10rad/m$⁠. Figure 5 shows the magnitude of the coefficients when $μ=0.9$⁠. The fall-off or rate of decay of the magnitude of the coefficients is similar in both cases. In Iserles et al. (2023), it is shown that the expansion coefficients for a wave packet exhibit asymptotic exponential decay.

Figure 6 shows u(x, T) when the initial condition is the wavepacket. In this case, there is a more visual change when the fractional derivative is included with increased damping. However, when the spatial frequency is increased to $ω=20rad/m$⁠, the situation is reversed and there is less damping as shown in Figure 7. When the frequency is further increased to $ω=40rad/m$⁠, the damping is again increased as shown in Figure 8. Varying a₀ in equation (2.10) will change the spatial frequencies at which this reduced-damping effect occurs.

The paper has shown that an implementation of a spectral method using the Malmquist-Takenaka functions is effective in simulating the behaviour of wave-packets of increasing spatial frequency on a nonlinear transmission line. The form of the Malmquist-Takenaka functions enables computational efficiencies to be achieved. The differentiation matrix is a narrow-band matrix. Furthermore, expansion coefficients can be computed using the FFT. In addition, it is not necessary to transform between the spectral and physical space. The number of functions required to represent a wavepacket increases linearly with frequency unlike Hermite functions where the number increases quadratically (Iserles et al., 2023).

Fractional derivatives are seen to have a significant effect on the evolution of wavepackets on the line with the spatial frequency of the wavepacket significantly influencing the nature of the behaviour. This feature is important as fractional derivatives are increasingly seen as more accurate models for THz applications (Antonini et al., 2021; Shang et al., 2013).

From equations 2.8 and 2.9:

Expanding this:

If α and $∂u∂t$ are of order $ϵ≪1$⁠, then the second last term can be neglected. If $∂2u∂t2$ is also of order $ϵ2$⁠, then the last term can be neglected:

Let $R=r/δ,C0=c0/δ,L=l/δ$ be the per-unit length parameters: