Next  |  Prev  |  Up  |  Top  |  Index  |  JOS Index  |  JOS Pubs  |  JOS Home  |  Search


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{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} (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}}$ (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. \begin{figure}[h]
\begin{center}
\begin{picture}(560,200)
\par
\put(30,10){\ep...
... $a=0.62$, at the stability bound, for $\lambda = 1/\sqrt{2a}$.}}
\end{figure}

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

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

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

$\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 (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) (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

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

$\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 (17) and the passivity constraint above is striking; these bounds are graphed in Figure 3. \begin{figure}[h]
\begin{center}
\begin{picture}(250,130)
\par
\put(0,0){\epsf...
... that there is a passive realization only for $0\leq a \leq 1$}.}
\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 §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,

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



Subsections
Next  |  Prev  |  Up  |  Top  |  Index  |  JOS Index  |  JOS Pubs  |  JOS Home  |  Search

Download vonn.pdf

``Spectral Analysis of Finite Difference Meshes'', by .
Copyright © 2005-12-28 by Julius O. Smith III<jos_email.html>
Center for Computer Research in Music and Acoustics (CCRMA),   Stanford University
CCRMA  [Automatic-links disclaimer]