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 [4] 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}$$ (16)

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

$$a+2b = 1\hspace{1.0in}4a+4b+c = \frac{2}{\lambda^{2}}$$ (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 [4], 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.

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 (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 [1] 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 (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.$$ (18)

It is interesting to look at the interpolated scheme from a waveguide mesh point of view (see Chapter 4 of [2] 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 (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 (16), this is true only for $a\leq 1$. Requiring $Y_{c}\geq 0$ is equivalent to requiring finally, from the second of constraints (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 (17) and the passivity constraint above is striking; these bounds are graphed in Figure 3.
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 §4.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 forego 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 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 (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.