Numerical Phase Velocity and Parasitic Modes

Because, in general, the image MDWDF of a given MDKC for a system of PDEs is a multi-step numerical integration scheme, it is reasonable to expect that parasitic modes [176] will be present in the solution. Energy in such modes often travels at speeds other than the desired wave speed in the medium, and may be highly oscillatory. If the scheme is consistent with the original system of PDEs, and stable, as is an MDWD network derived from the equivalent MDKC under the application of the trapezoid rule, then these parasitic modes must disappear in the limit as the time step is decreased (by the Lax-Richtmeyer Equivalence Theorem [176]). They have not as yet been addressed in the wave digital theory, and the subject is related to how initial conditions should be set in a MDWDF. The subject of initialization has been touched on only very briefly in [106].

Analysis of parasitic modes is easiest in the constant-coefficient case. We will examine the simplest possible non-trivial MDWDF, namely that of the constant-coefficient lossless source-free (1+1)D transmission line. Because at the stability limit, this scheme becomes equivalent to simple centered differences (see previous section), for which we do not have parasitic modes at all, we will look at the MDWDF of Figure 3.14(b) away from this limit. We have chosen $r_{0} = \sqrt{l/c}$. The MDWDF is redrawn in Figure 3.21, where we have

$$R_{1}=R_{2} = \frac{2}{\Delta}\left(v_{0}l-r_{0}\right)\hspace{1.0in}R_{0} = \frac{2r_{0}}{\Delta}$$

for some $v_{0}>1/\sqrt{lc}$.

Figure 3.21: Steady-state MDWD network for the lossless, source-free constant-coefficient (1+1)D transmission-line equations.

Note that because the system is now linear and shift-invariant, we have replaced the shifts ${\bf T}_{1}$ and ${\bf T}_{2}$ in the two directions $t_{1}$ and $t_{2}$ by their frequency domain counterparts $z_{1}^{-1}$ and $z_{2}^{-1}$. Recall also that we have, from (3.47) and (3.48), that

$$z_{1}^{-1} = z^{-1}w^{-1}\hspace{1.0in}z_{2}^{-1} = z^{-1}w$$

where $z^{-1}$ represents a unit delay in the time direction, and $w$ a unit shift in the $x$ direction. We have written the outputs of the delay registers, in an exponential state, as

$$x_{j}(k\Delta, nT) = \hat{x}_{j}z^{n}w^{k},\hspace{0.2in} {\rm for}\hspace{0.2in} j=1, \hdots, 4$$

where the $\hat{x}_{j}$ are complex amplitudes. The updating of the values in the delay registers can be written, in terms of these amplitudes, as

$$\begin{bmatrix}\hat{x}_{1}\\ \hat{x}_{2}\\ \hat{x}_{3}\\ \hat{x}_... ...n{bmatrix}\hat{x}_{1}\\ \hat{x}_{2}\\ \hat{x}_{3}\\ \hat{x}_{4}\\ \end{bmatrix}$$

which is parametrized by a reflectance

$$\alpha = \frac{R_{0}-R_{1}}{R_{0}+R_{1}}$$ (3.83)

If we introduce the variables

$$y_{1} = x_{1}+x_{4}\hspace{1.0in}y_{4} = x_{1}-x_{4}$$

the updating decouples into two subsystems, namely

${\bf A}_{21}$ and ${\bf A}_{34}$ are known as spectral amplification matrices (see Appendix A).

The symbols [176] of the two subsystems, ${\bf Q}_{24}$ and ${\bf Q}_{31}$ are defined by

$${\bf Q}_{21} = {\bf I}_{2}-z^{-1}{\bf A}_{21}\hspace{0.5in}{\bf Q}_{34} = {\bf I}_{2}-z^{-1}{\bf A}_{34}$$

where ${\bf I}_{2}$ is the $2\times 2$ identity matrix. Nontrivial solutions to the update equations (3.83) occur when the determinants of the symbols vanish. In the absence of boundary conditions, we may assume $w=e^{j\beta\Delta}$, where $\beta$ is a real wavenumber, in which case we have four solutions in terms of $z$ given by

which are simply the eigenvalues of the spectral amplification matrices. The corresponding eigenvectors of these same matrices are

$${\bf u}_{21,\pm} = \begin{bmatrix}\alpha\cos(\frac{\beta\Delta}{2... ...}(\frac{\beta\Delta}{2})}\\ (1-\alpha)e^{\frac{\beta\Delta}{2}}\\ \end{bmatrix}$$

All four eigenvalues are of unit magnitude, and thus, using $z = e^{j\omega T}$, we can rewrite solutions (3.84) as

$$e^{j\omega_{21\pm}T} = \pm e^{\frac{j\Delta}{2}\left(\beta \pm \n... ...7in}e^{j\omega_{34\pm}T} = \pm e^{-\frac{j\Delta}{2}\left(\beta \pm \nu\right)}$$ (3.86)

for some real $\nu$ defined by $\sin(\nu\Delta/2) = \alpha\sin(\beta\Delta/2)$. ($\nu$ always exists because we have $\vert\alpha\vert\leq 1$, from (3.82).) For small wavenumbers, we have

$$\nu \approx \alpha\beta$$

and we thus have in this limit, for the roots subscripted with $+$ in (3.85),

$$e^{j\omega_{21+}T} \approx e^{\frac{j\Delta (1+\alpha)}{2}\beta} \Rightarrow \frac{\omega_{21+}}{\beta}\approx \frac{1}{\sqrt{lc}}$$

$$e^{j\omega_{34+}T} \approx e^{-\frac{j\Delta (1+\alpha)}{2}\beta} \Rightarrow \frac{\omega_{34+}}{\beta}\approx -\frac{1}{\sqrt{lc}}$$

where we have used the fact that $(1+\alpha)/2 = 1/(v_{0}\sqrt{lc})$, which follows from (3.82) and the definitions of the port resistances in (3.64) as well as $v_{0} = \Delta/T$. The quantities $\omega_{21+}/\beta$ and $\omega_{34+}/\beta$ are called numerical phase velocities [176]; they approach the propagation speed in the medium, from (3.58), and these two solutions are to be interpreted as approximations to the traveling wave solution to the transmission line equations. The other modes, however, are parasitic, in that they do not propagate near the physical velocity. They are not problematic, provided initial conditions are set properly; indeed, in the limit as $\Delta$ becomes small, any reasonable initial conditions tend to align the system with the dominant traveling modes of the system.

Clearly, if we are at the passivity limit, where $v_{0}=1/\sqrt{lc}$, then $R_{1}=0$, and thus $\alpha = 1$, which implies, finally that $\nu = \beta$, so that we have, from (3.85), that $\omega_{21+}/\beta = 1/\sqrt{lc}$ and $\omega_{34+}/\beta = -1/\sqrt{lc}$; wave propagation is thus dispersionless. As mentioned in the previous section, at this limit, the MDWD network reduces to an exact digital traveling wave solution (this was also noted in §3.7.5). It is also interesting to note that when $R_{1} = R_{2} = R_{0}$, so that $\alpha$ and $\nu$ are zero, then (3.85) implies that wave propagation is also dispersionless in this case as well. It is easy to see here, from Figure 3.21, that because $R_{1} = R_{2} = R_{0}$, there will be no scattering through the adaptors; the pure time delays may thus be shifted directly into the lattice two port, and we can perform a manipulation similar to that of §3.7.5 to give a simplified digital ``traveling wave'' network, with doubled time delays. Here, we are in effect implementing a traveling wave solution on a different grid, but the implication is that for the corresponding problem with material variation, the MDWD network gives a good approximation to the numerical phase velocity even for certain values of $v_{0}$ which are far from the local physical wave speed. This is not true for digital waveguide networks, where the numerical phase velocities degrade considerably away from the passivity limit. We will return to these expressions (which provide complete information regarding the numerical dispersion properties of the scheme) in §4.3.8 in a comparison with the digital waveguide network for the same system. In anticipation of the discussion in §3.10, we mention that for constant $v_{0}$, we have for the eigenvectors corresponding to the dominant modes, that

$$\lim_{\Delta\rightarrow 0} {\bf u}_{21 +} = \lim_{\Delta\rightarrow 0}{\bf u}_{34,+} = \begin{bmatrix}\alpha +1\\ 1-\alpha\\ \end{bmatrix}$$

Because we also have, from Figure 3.21, that $\hat{a}_{1} = -\hat{x}_{1} = -(\hat{y}_{1}+\hat{y}_{4})/2$, and $\hat{a}_{2} = -(\hat{x}_{2}+\hat{x}_{3})/2$, we can also write, for the dominant mode,

$$\lim_{\Delta\rightarrow 0} \begin{bmatrix}\hat{a}_{1}\\ \hat{a}_{... ...d{bmatrix} = -\frac{2}{R_{1}+R_{0}}\begin{bmatrix}R_{1}\\ R_{0}\\ \end{bmatrix}$$

Thus in this limit, the wave variables incident on the left adaptor occur in the same ratio as the port resistances, and are in fact aligned with an eigenvector of the scattering matrix corresponding to the adaptor. A similar statement holds for the quantities incident on the right adaptor. We will return to this observation in the next section.