next up previous
Next: Optimally direction-independent numerical dispersion Up: Finite Difference Schemes for Previous: The Rectilinear Scheme

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{displaymath}\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}\end{displaymath} (A.16)

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

$\displaystyle 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}$.
\end{picture} \end{center} \vspace{0.5in}

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

$\displaystyle B_{\mbox{{\scriptsize\boldmath$\beta$}}} = -2\lambda^{2}\Big(a(\c...

and thus

$\displaystyle F_{\mbox{{\scriptsize\boldmath$\beta$}}} = a\big(\cos(\beta_{x}\D...

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]$:

$\displaystyle 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

$\displaystyle \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

$\displaystyle \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.$   (for Von Neumann stability) (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

$\displaystyle \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

$\displaystyle 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

$\displaystyle \lambda\leq \frac{1}{\sqrt{1+a}}{\rm ,}\hspace{0.3in}0\leq a \leq 1$   (for passivity)    

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$.
...ut(130,100){\mbox{\small {Not Passive}}}
\end{picture} \end{center} \end{figure}

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,

$\displaystyle \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

$\displaystyle \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.

next up previous
Next: Optimally direction-independent numerical dispersion Up: Finite Difference Schemes for Previous: The Rectilinear Scheme
Stefan Bilbao 2002-01-22