The Interpolated Rectilinear Scheme

This scheme, like the standard rectilinear scheme, is defined over a grid with indices $i$ and $j$, for points with $x=i\Delta$ and $y=j\Delta$. Updating, in this case, at a given point, requires access to values of the grid function at the previous time step at nearest-neighbor grid points to the north, east, west and south, as well as those to the north-east, north-west, south-east and south-west, which are more distant by a factor of $\sqrt{2}$. The scheme is referred to as ``interpolated'' in [157] because it is derived as an approximation to a hypothetical (and non-realizable) multi-directional difference scheme with minimally directionally-dependent numerical dispersion. (It is perhaps more useful to think of the scheme as interpolating between two rectilinear schemes operating on grids with a relative angle of 45 degrees.) The difference scheme will have the form

$$\begin{split}U_{i,j}(n+1)+U_{i,j}(n-1) &= \lambda^{2}a\Big(U_... ..._{i-1,j-1}(n)\Big)\\ &\quad+ \lambda^{2}cU_{i,j}(n) \end{split}$$ (A.16)

for constants $a$, $b$ and $c$ which satisfy the constraints

$$a+2b = 1\hspace{1.0in}4a+4b+c = \frac{2}{\lambda^{2}}$$ (A.17)

for consistency with the wave equation. If $b=0$, we get the standard rectilinear scheme, and if $a=0$, we get a rectilinear scheme operating on a grid of spacing $\sqrt{2}\Delta$, which is rotated by 45 degrees with respect to that of the standard scheme. This general form was put forth in [157], and the free parameter $a$ may be adjusted to give a less directionally-dependent numerical phase velocity; it may thus be used in conjunction with frequency-warping methods for reducing dispersion error. In general, the interpolated scheme cannot be decomposed into mutually exclusive subschemes.

Figure A.2: The interpolated rectilinear scheme (A.15)-- (a) numerical grid and connections for the interpolated rectilinear scheme (A.15); (b) $v_{\mbox{{\scriptsize\boldmath$\beta$}}, phase}/\gamma$ of the scheme for $a=0.62$ at the ``passivity'' bound, $\lambda = 1/\sqrt{1+a}$; (c) $v_{\mbox{{\scriptsize\boldmath$\beta$}}, phase}/\gamma$ for $a=0.62$, at the stability bound, for $\lambda = 1/\sqrt{2a}$.

It is possible to examine the stability of this method as in the previous case. We again have an amplification polynomial equation of the form of (A.5), with

$$B_{\mbox{{\scriptsize\boldmath$\beta$}}} = -2\lambda^{2}\Big(a(\c... ...os(\beta_{y}\Delta))+(1-a)\cos(\beta_{x}\Delta)\cos(\beta_{y}\Delta)-1-a\Big)-2$$

and thus

$$F_{\mbox{{\scriptsize\boldmath$\beta$}}} = a\big(\cos(\beta_{x}\D... ...+\cos(\beta_{y}\Delta)\big)+(1-a)\cos(\beta_{x}\Delta)\cos(\beta_{y}\Delta)-1-a$$

Note that $F_{\mbox{{\scriptsize\boldmath $\beta$}}}$ is multilinear [3] in $\cos(\beta_{x}\Delta)$ and $\cos(\beta_{y}\Delta)$, so that any extrema must occur at the corners of the region in the spatial frequency plane defined by $\vert\cos(\beta_{x}\Delta)\vert\leq 1$, and $\vert\cos(\beta_{y}\Delta)\vert\leq 1$. Thus, we need evaluate $F_{\mbox{{\scriptsize\boldmath $\beta$}}}$ only for $\beta$ $^{T}=[\beta_{x}, \beta_{y}] = [0,0]$, $[\pi/\Delta, 0]$, $[0,\pi/\Delta]$ and $[\pi/\Delta, \pi/\Delta]$:

$$F_{\mbox{{\scriptsize\boldmath$\beta$}}^{T}=[0,0]} = 0\hspace{0.5... ....5in}F_{\mbox{{\scriptsize\boldmath$\beta$}}^{T}=[\pi/\Delta,\pi/\Delta]} = -4a$$

The global maximum of $F_{\mbox{{\scriptsize\boldmath $\beta$}}}$ is non-positive (and thus condition (A.8) is satisfied) only if $a\geq 0$. The global minimum of $F_{\mbox{{\scriptsize\boldmath $\beta$}}}$, over this range of $a$ will then be

$$\min_{\mbox{{\scriptsize\boldmath$\beta$}}}F_{\mbox{{\scriptsize\... ... a\leq\frac{1}{2}\\ -4a, & \hspace{0.3in}a\geq \frac{1}{2}\\ \end{array}\right.$$

and the stability bound on $\lambda$ will be

$$\lambda \leq \left\{\begin{array}{lr} 1, & \hspace{0.3in}0\leq a\... ...}\\ \frac{1}{\sqrt{2a}}, & \hspace{0.3in}a\geq \frac{1}{2}\\ \end{array}\right.$$ (A.18)

It is interesting to look at the interpolated scheme from a waveguide mesh point of view (see Chapter 4 for details). At each grid point we will have a nine-port parallel scattering junction; four connections are made to neighboring points to the north, south, east and west, through a unit-delay bidirectional delay line of admittance $Y_{a}$, four more connections are made to the points to the north-east, south-east, north-west and south-west using waveguides of admittance $Y_{b}$, and there will be a self-loop of admittance $Y_{c}$. If the junction voltage is written as $U_{i,j}(n)$, then the difference scheme corresponding to this waveguide mesh will be exactly (A.15), with

$$\lambda^{2}a = \frac{2Y_{a}}{Y_{J}}\hspace{0.3in}\lambda^{2}b = \frac{2Y_{b}}{Y_{J}}\hspace{0.3in}\lambda^{2}c = \frac{2Y_{c}}{Y_{J}}$$

where the junction admittance $Y_{J}$ (assumed positive) will be given by

$$Y_{J} = 4Y_{a}+4Y_{b}+Y_{c}$$

The passivity condition will then be a condition on the positivity of $Y_{a}$, $Y_{b}$ and $Y_{c}$. From the previous discussion, we already require $a\geq 0$, so this ensures that $Y_{a}\geq 0$. Requiring $Y_{b}\geq 0$ is equivalent to requiring $b\geq 0$; from the first of constraints (A.16), this is true only for $a\leq 1$. Requiring $Y_{c}\geq 0$ is equivalent to requiring finally, from the second of constraints (A.16), that

$$\lambda\leq \frac{1}{\sqrt{1+a}}{\rm ,}\hspace{0.3in}0\leq a \leq 1$$

The difference between the constraints for stability from (A.17) and the passivity constraint above is striking; these bounds are graphed in Figure A.3.

Figure A.3: Stability bounds for the interpolated rectilinear scheme, as a function of the free parameter $a$. The solid line indicates the maximum value of $\lambda$ for a given value of $a$, and the dashed line the maximum value of $\lambda$ allowed in a passive waveguide mesh implementation. Note that there is a passive realization only for $0\leq a \leq 1$.

This is not the last time that we will find a discrepancy between Von Neumann stability of a scheme and passivity of the related mesh structure; it will come up again in the following section during a discussion of the triangular scheme, and in §A.3.3 when we look at the (3+1)D interpolated scheme. It is interesting to note that for a given value of $a$, with $0\leq a \leq 1$, the numerical dispersion properties can always be improved if we are willing to forgo passivity (and a mesh implementation). We have plotted the numerical phase velocities of this scheme for $a=0.62$, at both the stability limit and the passivity limit in Figure A.2.

Finally, we mention that the computational and add densities for this scheme will be, in general,

$$\rho_{interp} = \frac{v_{0}}{\Delta^{3}}\hspace{1.5in} \sigma_{interp} = \frac{10v_{0}}{\Delta^{3}}$$

over the range of $v_{0}$ allowed by the stability constraint (A.17). For the scheme at the passivity bound (for $\lambda = 1/\sqrt{1+a}$, with $0< a< 1$), we have

$$\rho_{interp}^{p} = \frac{\gamma\sqrt{1+a}}{\Delta^{3}}\hspace{1.5in} \sigma_{interp}^{p} = \frac{9\gamma\sqrt{1+a}}{\Delta^{3}}$$

We recall that for $a=0$ or $a=1$, at the stability limit, we again have the standard rectilinear scheme, for which a grid decomposition is possible; this was discussed in the previous section.