Numerical Anisotropy in Finite Differencing

TL;DR

This paper reviews numerical anisotropy in finite difference methods for wave equations, highlighting how multi-dimensional discretization causes directional errors and discussing optimization strategies to mitigate these effects.

Contribution

It provides a comprehensive overview of numerical anisotropy in finite difference schemes and discusses recent optimization approaches to reduce this error in multi-dimensional wave simulations.

Findings

Numerical anisotropy causes directional wave speed errors.

Optimization techniques can significantly reduce anisotropy.

Review of key studies on finite difference scheme improvements.

Abstract

Numerical solutions to hyperbolic partial differential equations, involving wave propagations in one direction, are subject to several specific errors, such as numerical dispersion, dissipation or aliasing. In multi-dimensions, where the waves propagate in all directions, there is an additional specific error resulting from the discretization of spatial derivatives along grid lines. Specifically, waves or wave packets in multi-dimensions propagate at different phase or group velocities, respectively, along different directions. A commonly used term for the aforementioned multidimensional discretization error is the numerical anisotropy or isotropy error. In this review, the numerical anisotropy is briefly described in the context of the wave equation in multi-dimensions. Then, several important studies that were focused on optimizations of finite difference schemes with the objective of reducing the numerical anisotropy are discussed.

Figures

Tables

Table 1. Table 1: Weights of the selected spatial finite difference stencils

Stencil	${\alpha }_1$	${\alpha }_2$	$a_1$	$a_2$	$a_3$
$E2$	0	0	1/2	0	0
$E4$	0	0	2/3	-1/12	0
$E6$	0	0	3/4	-3/20	1/60
$DRP$	0	0	0.770882380	-0.166705904	0.020843142
$C4$	1/4	0	3/4	0	0
$Haras$	0.3534620	0	1.5669657/2	0.13995831/4	0
$Lui$	0.5381301	0.0666331	1.36757772/2	0.823428170/4	0.0185207834/6
$Lele$	0.5771439	0.0896406	1.3025166/2	0.99355/4	0.03750245/6

Equations

$$\sum _{k=1}^{N_c}{\alpha }_k(u_{j+k}^{\prime }+u_{j-k}^{\prime })+u_j^{\prime }=\frac{1}{h}\left(\sum _{k=1}^{N_e}{a}_k(u_{j+k}-u_{j-k})\right)+O(h^n),$$

$$a{u}_{j+1}^{F^{\prime }}+c{u}_{j-1}^{F^{\prime }}+(1-a-c)u_j^{F^{\prime }}=\frac{1}{h}\left[b{u}_{j+1}-(2b-1)u_j-(1-b)u_{j-1})\right],$$

$$c{u}_{j+1}^{B^{\prime }}+a{u}_{j-1}^{B^{\prime }}+(1-a-c)u_j^{B^{\prime }}=\frac{1}{h}\left[(1-b)u_{j+1}-(2b-1)u_j-b{u}_{j-1})\right],$$

$$({\hat{u}}_m^{\prime }{)}_{num}=iK(k_{m}h){\hat{u}}_m,$$

$$K(z)=\frac{{\sum }_{n=1}^{N_e}2a_n\sin (nz)}{1+{\sum }_{n=1}^{N_c}2{\alpha }_n\cos (nz)}.$$

$${\partial }_{t}u=\text{c}\nabla u,$$

$$d_{t}u=-\frac{c}{2h}\big{[}\cos\alpha(u_{{i+1},j}-u_{{i-1},j})+\sin\alpha(u_{i,{j+1}}-u_{i,{j-1}})\big{]},$$

$$\tilde{u}(\xi ,\eta ,\omega )=\frac{1}{(2\pi {)}^3}{\int }_0^{\infty }\int {\int }_{-\infty }^{\infty }u(x,y,t)e^{-i(\xi x+\eta y-\omega t)}dxdydt$$

$$\omega =cK({\cos }^2\alpha +{\sin }^2\alpha )=cK.$$

$$\omega=\frac{c}{h}\big{[}\cos\alpha\sin(Kh\cos\alpha)+\sin\alpha\sin(Kh\cos\alpha)\big{]}$$

$$c_{n}=\frac{\omega}{K}=\frac{c}{Kh}\big{[}\cos\alpha\sin(Kh\cos\alpha)+\sin\alpha\sin(Kh\cos\alpha)\big{]}.$$

$$g_{n}=\partial_{K}\omega=c\big{[}\cos^{2}\alpha\cos(Kh\cos\alpha)+\sin^{2}\alpha\cos(Kh\sin\alpha)\big{]},$$

$${\partial }_{tt}u={\partial }_{xx}u+{\partial }_{yy}u,$$

$$u_{ij}^{n+1}-u_{ij}^n+u_{ij}^{n-1}=\frac{k^2}{h^2}(u_{i+1,j}^n+u_{i-1,j}^n+u_{i,j+1}^n+u_{i,j-1}^n-4u_{i,j}^n)$$

$${\left({\partial }_{x}u\right)}_{i,j}=\frac{1}{h(1+\beta )}\sum _{\nu =-M}^{\nu =M}{a}_{\nu }\left({\text{E}}_x^{\nu }+\frac{\beta }{2}{\text{D}}_x\right)\cdot u_{i,j}$$

$$(\xi h)_{opt}^{*}=\frac{2}{(1+\beta)}\sum_{n=-N}^{M}a_{n}\Big{\{}e^{nI\xi h}+\frac{\beta}{2}\big{[}e^{nI(\xi+\eta)h}+e^{nI(\xi-\eta)h}\big{]}\Big{\}},$$

$${\nabla }^{2}p+{\xi }^{2}p=f$$

$$\frac{p_{i+1,j}-2p_{i,j}+p_{i-1,j}}{h^2}+\frac{p_{i,j+1}-2p_{i,j}+p_{i,j-1}}{h^2}+{\xi }^2{p}_{i,j}=f_{i,j}$$

$$p_a(r,\theta {)}_{r_{ij}\to \infty }=\left(\frac{2\pi }{\xi }\right)\frac{\pi }{i{r}^{1/2}}{e}^{i(\xi r-\pi /4)}\bar{F}({\bar{\alpha }}_s,{\bar{\beta }}_{+}({\bar{\alpha }}_s))+O(r^{-3/2})$$

$$p_n(r_{ij},{\theta }_{ij}{)}_{r_{ij}\to \infty }=\frac{e^{i{K}_{ij}{r}_{ij}}}{r_{ij}^{1/2}}\left[G_0({\theta }_{ij}+\frac{G_1({\theta }_{ij}}{r_{ij}})+O(r_{ij}^{-5/2})\right]$$

$$D(\theta ,\xi h)=\frac{|p_a|}{|p_n|}$$

$${\partial }_{t}u+{\partial }_{x}f=0,f=cu,c>0$$

$$d_t{\bar{u}}_i+{\bar{f}}_{i+1/2}-{\bar{f}}_{i-1/2}=0$$

$$\alpha v_{i-1/2}+v_{i+1/2}+\alpha v_{i+3/2}=b\frac{v_{i+5/2}-v_{i-3/2}}{4h}+a\frac{v_{i+3/2}-v_{i-1/2}}{2h}$$

$$a=\frac{2(2+\alpha )}{3},b=\frac{-1+4\alpha }{3}$$

$$\hat{A}(w)=\frac{i\left(a\sin (w)+b\sin (2w)/2\right)}{1+2\alpha \cos w}$$

$$w_i(w,\theta )=\cos (\theta )w_d(w\cos (\theta ))+\sin (\theta )w_d(w\sin (\theta ))$$

$$E_i(\alpha ,w_{max})={\int }_0^{w_{max}}{\int }_0^{\pi /2}|w_i-w|d\theta dw$$

$$u_{i,j}^{F^{\prime }}$$

$$u_{i,j}^{B^{\prime }}$$

$$u_{i,j}^{F^{\prime }}$$

Full text

1 Introduction

Numerical anisotropy is a discretization error that is specific to numerical approximations of multidimensional hyperbolic partial differential equations (PDE). This error is often neglected, and the focus is directed toward the reduction of other types of discretization errors, such as numerical dissipation, dispersion or aliasing (e.g., Lele [23], Tam and Webb [43], Kim and Lee [18], Zingg et al. [48], Mahesh [27], Hixon [10], Ashcroft and Zhang [2], Fauconnier et al. [7] or Laizet and Lamballais [22]), or toward improving the accuracy of various time marching schemes (e.g., Hu et al. [13], Stanescu and Habashi [38], Mead and Renaut [28], Bogey and Bailly [5] or Berland et al. [4]). There are several areas, however, where the numerical anisotropy can significantly affect the numerical solution based on finite difference or finite volume schemes (example include computational acoustics, computational electromagnetics, elasticity or seismology). The numerical anisotropy can be reduced by using, for example, one-dimensional high-resolution discretization schemes, multidimensional optimized difference schemes, or sufficiently fine grids. However, by increasing the number of grid points the computational time may increase considerably, while one-dimensional high-resolution difference schemes may generate spurious waves at the boundaries of the domain. Oftentimes, optimizations of multidimensional difference schemes are more effective.

High-order finite difference schemes that are optimized in one-dimension may not preserve their wavenumber resolution in multi-dimensional problems. These schemes may experience numerical anisotropy, because the dispersion characteristics along grid lines may not be the same as the dispersion characteristics associated with the diagonal directions. Over the years, several attempts to reduce the numerical anisotropy by various techniques were reported. A comprehensive analysis of the numerical anisotropy was performed in the book of Vichnevetsky [45] where, among others, the two-dimensional wave equation was solved using two different finite difference schemes for the Laplacian operator. A considerable reduction of the numerical anisotropy was attained by weight averaging the two schemes. A slightly similar approach was previously used by Trefethen [44] who used the leap frog scheme to solve the wave equation in two dimensions. Zingg and Lomax [47] performed optimizations of finite difference schemes applied to regular triangular grids, that give six neighbor points for a given node. They conducted comparisons between the newly derived schemes and conventional schemes that were discretized on square grids, and found that the numerical anisotropy can be significantly reduced by using triangular grids. Tam and Webb  [42] proposed an anisotropy correction to the finite difference representation of the Helmholtz equation. They derived an anisotropy correction factor using asymptotic solutions to the continuous equation and its finite difference approximation.

Jo et al. [15], in the context of solving the acoustic wave equation, proposed a finite difference scheme over a stencil consisting of grid points from more than one direction, by linearly combining two discretizations of the second derivative operator. A notable reduction of the numerical anisotropy was obtined, but the numerical dispersion error was increased. Hustesdt et al. [14] proposed a two-staggered-grid finite difference schemes for the acoustic wave propagation in two dimensions, where the first derivative operator was discretized along the grid line and along the diagonal direction. Lin and Sheu [24] explored the dispersion-relation-preserving concept of Tam and Webb [43] in two dimensions to optimize the first-order spatial derivative terms of a model equation that resembles the incompressible Navier-Stokes momentum equation. They approximated the derivative using a nine-point grid stencil, resulting in nine unknown coefficients. Eight of them were determined by employing Taylor series expansions, while the ninth one was determined by requiring that the two-dimensional numerical dispersion relation is the same as the exact dispersion relation.

Kumar [21] derived isotropic finite difference schemes for the first and second derivatives in the context of symmetric dendritic solidification, and obtained a notable reduction of the numerical anisotropy. Patra and Kartunnen [29] introduced several finite difference stencils for the Laplacian, Bilaplacian, and gradient of Laplacian, with the objective of improving the isotropic characteristics. Their stencils consisted of more grid points than the conventional schemes, but it was shown that the computational cost may decrease with more than 20$\%$ due to some gain in terms of stability. Stegeman et al. [39] applied spectral analysis to evaluate the error in numerical group velocity (both the magnitude and the direction) of vorticity, entropy and acoustic waves, using the numerical solution to the linearized Euler equations in two-dimensions. They showed that a different measure of the group velocity error must be used to account for the error in the propagation direction of the waves. They also stressed that the numerical group velocity is more important than the numerical phase velocity in analyzing the errors associated with wave propagation. In a series of papers [30, 31, 32, 33], Sescu et al. proposed a technique to derive finite difference schemes in multi-dimensions with improved isotropy. The optimization performed in [30, 31, 32, 33] improved the isotropy of the wave propagation and, moreover, the stability restrictions of the multidimensional schemes in combination with either Runge-Kutta or linear multistep time marching methods were found to be more effective. They found that the stability restrictions are more favorable when using multidimensional schemes, even if they involve more grid points in the stencils. However, this was advantageous for low order schemes, such as those of second or fourth order of accuracy, but it was also shown that favorable stability restrictions can be obtained for higher order of accuracy schemes (sixth or eight) by increasing the isotropy corrector factor. The approach was extended to prefactored compact schemes by Sescu and Hixon [35, 36]. Beside reducing the numerical anisotropy, the new multidimensional compact schemes are computationally cheaper than the corresponding explicit multidimensional scheme defined on a same stencil.

In computational electromagnetics, there were many attempts to reduce the numerical anisotropy, by applying various techniques. Berini and Wu [3] conducted a comprehensive analysis of the numerical dispersion and numerical anisotropy of finite difference schemes applied to transmission-line modeling (TLM) meshes. They found that, under certain circumstances, the time domain nodes introduce anisotropy into the dispersion characteristics of isotropic media, stressing the importance of developing schemes with improved isotropy. Gaitonde and Shang [8] proposed a class of high-order compact difference-based finite-volume schemes that minimizes the dispersion and isotropy error functions for the range of wavenumbers of interest. Sun and Trueman [40] proposed an optimization of two-dimensional finite difference schemes, by considering additional nodes surrounding the point of differencing. They obtained a significant reduction in the numerical anisotropy, dispersion error and the accumulated phase errors over a broad bandwidth. Further optimizations of this scheme were performed in another paper of Sun and Trueman [41]. Koh et al. [19] derived a two-dimensional finite-difference time-domain method, discretizing the Maxwell equations, to eliminate the numerical dispersion and anisotropy. They showed that the new algorithm has isotropic dispersion and resemble the exact phase velocity, whose isotropic property is superior to that of other existing schemes. Shen and Cangellaris [37] introduced a new stencil for the spatial discretization of Maxwell’s equations. Compared to conventional second-order accurate FDTD scheme, their scheme experienced superior isotropy characteristics of the numerical phase velocity. They also showed that the Courant number cab be increased by using the newly derived schemes. Kim et al. [17] derived new three-dimensional isotropic dispersion-finite-difference time-domain schemes (ID-FDTD) based on a linear combination of the traditional central difference equation and a new difference equation using extra sampling points. Among all versions of the proposed finite-difference schemes, three of them showed improved isotropy of the wave propagation compared to the original scheme of the Yee [46]. Kong and Chu [20] introduced a new unconditionally-stable finite-difference time-domain method with low numerical anisotropy in three-dimensions. Compared with other finite-difference time-domain methods, the normalized numerical phase velocity of their proposed scheme was significantly improved, while the dispersion error and numerical anisotropy have been reduced.

This review will describe and discuss the numerical anisotropy in the framework of wave equation, and will present some of the most important optimizations of finite difference schemes in the context of reducing the numerical anisotropy. In section II, the dispersion error and the numerical anisotropy existing in finite difference discretizations of the wave equation are introduced and discussed. In section III, several approaches to reduce the numerical anisotropy, that were developed over the years by various research groups, are reviewed and discussed. Concluding remarks are included in section IV.

2 Dispersion Error and Numerical Anisotropy

Let us consider the centered finite difference approximation of the spatial derivative, which contains both the explicit and the implicit (or compact) parts:

[TABLE]

where the gridfunctions are $u_{j}=u(x_{j})$ for $1\leq j\leq N$, the derivatives are denoted by a prime, $u_{j}^{\prime}$, $h$ is the space step, and $\alpha_{k}$ and $a_{k}$ are given coefficients. If $N_{c}=0$ the scheme is termed explicit, while compact schemes (also known as implicit or Pade schemes), by contrast, have $N_{c}\neq 0$ and require the solution of a matrix equation to determine the derivatives along a grid line. Conventionally, the coefficients $\alpha_{k}$ and $a_{k}$ are chosen to provide the largest possible exponent, $n$, in the truncation error, for a given stencil width, but in some instances some of these coefficients are determined to provide improved dispersion characteristics of the scheme. Table 1 includes some of these weights for various explicit and compact finite difference schemes: explicit classical second order scheme (E2), explicit classical fourth order scheme (E4), explicit classical sixth order scheme (E6), dispersion-relation-preserving scheme of Tam and Webb [43], compact classical fourth order scheme (C4), optimized tridiagonal compact scheme of Haras and Ta’asan [9] (Haras), optimized pentadiagonal scheme of Lui and Lele [25] (Lui) and spectral-like pentadiagonal compact scheme of Lele [23] (Lele). The prefactored compact scheme of Hixon [10, 11] is also included here in the form:

[TABLE]

where $F$ and $B$ stand for ’forward’ and ’backward’, respectively (in a predictor-corrector time marching framework). For sixth order accuracy, $a=1/2-1/(2\sqrt{5})$, $b=1-1/(30a)$ and $c=0$. The leading order term in the truncation error of a finite difference scheme depends on the choice of the coefficients and the $(n+1)$st derivative of the function $u$.

To study the wavenumber characteristics of finite difference schemes, consider a periodic domain in real space, $x\in[0,L]$, with $N$ uniformly spaced points(the spatial step size is $h=L/N$). The discrete Fourier transform of $u$ is given as $\hat{u}_{m}=\frac{1}{N}\sum_{j=1}^{N}u_{j}e^{-ik_{m}x_{j}}$ with $m=-N/2,...,N/2-1,$ where the wavenumber is $k_{m}=2\pi m/L$. The $m$th component of the discrete Fourier transform of $u^{\prime}$ denoted $\hat{u}_{m}^{\prime}$ is simply $ik_{m}\hat{u}_{m}$. Taking the discrete Fourier transform of equation (1) implies that

[TABLE]

where the numerical wavenumber is given as

[TABLE]

Figure 1 shows the numerical wavenumber for various explicit and compact schemes, corresponding to those given in table 1. The numerical wavenumber is compared to the analytical wavenumber which is represented by the straight line in figure 1. As one can notice, the compact schemes are superior to the explicit schemes; however, compact schemes are computationally more demanding because large matrices have to be inverted.

In muldimensions, the numerical wavenumber and the numerical phase and group velocity are also dependent on the direction of propagation. Figure 2 shows the numerical wavenumber surface for the wave equation in two dimension, corresponding to schemes E2, E6 and Hixon as given in table 1 and equation (2), respectively. The cone represents the exact wavenumber surface, obtained by revolving the straight line from figure 1 around the vertical axis. One can clearly notice the anisotropy in the numerical wavenumber surfaces associated with the finite differencing.

A simple way to reveal the numerical anisotropy is by considering the advection equation in two dimensions,

[TABLE]

with the initial condition $u(\textbf{r},0)=u_{0}(\textbf{r})$, where $\textbf{r}=(x,y)$ is the vector of spatial coordinates, $\textbf{c}=c(\cos\alpha\hskip 5.69054pt\sin\alpha)$ is the velocity vector ($c$ is a scalar and $\alpha$ the propagation direction angle), $\nabla=(\partial_{x}\hskip 5.69054pt\partial_{y})^{T}$ and $u(\textbf{r},t)$ and $u_{0}(\textbf{r})$ are scalar functions. A simple semi-discretization of equation (5) on a square grid is obtained as

[TABLE]

where $h$ is the grid step. Consider the Fourier-Laplace transform:

[TABLE]

where $\xi=K\cos\alpha$ and $\eta=K\sin\alpha$ are the components of the wavenumber and $\omega$ is the frequency ($K$ is the wavenumber magnitude). The application of Fourier-Laplace transform to equation (5) gives the exact dispersion relation:

[TABLE]

The exact phase velocity is given by $c_{e}=\omega/K=c$. By substituting $\omega$ in equation (7) with (8), $u(\textbf{r},t)$ is obtained as a superposition of sinusoidal solutions in the plane with constant phase lines given by $x\cos\alpha+y\sin\alpha-c_{e}t=const$. As one can notice, the exact phase velocity $c_{e}$ does not depend on the propagation direction $\alpha$, which means that the wave propagates with the same phase velocity in all directions (it is isotropic). Moreover, the exact group velocity defined as $g_{e}=\partial\omega/\partial K=c$ is the same as the exact phase velocity because the dispersion relation is a linear function of $K$.

We now apply the same Fourier-Laplace transform to the numerical approximation (6) and obtain the numerical dispersion relation in the form

[TABLE]

The numerical phase velocity will be given as

[TABLE]

The constant phase lines are expressed by the equation $x\cos\alpha+y\sin\alpha-c_{n}t=const$ and move with the phase velocity $c_{n}$. The numerical anisotropy is revealed in equation (10) by the dependence of the numerical phase velocity on the propagation direction angle $\alpha$. In addition, the numerical group velocity is different from the numerical phase velocity (while previously, in the continuous case, they were the same),

[TABLE]

which is also dependent on the propagation direction. This directional dependence of both phase and group velocities defines the numerical anisotropy. As an illustration, figure 3 shows polar diagrams for two typical schemes, fourth order explicit E4 and sixth order compact C6 schemes, revealing the numerical anisotropy (the circle of radius $1$ in figure 3 represents the exact solution).

3 Reduction of the Numerical Anisotropy

In this section, several attempts to reduce the numerical anisotropy, performed by various research groups over the years, are briefly reviewed. The optimizations of the schemes are grouped according to the mathematical model: wave equation, Helmholtz equations, advection equation, Maxwell equation, and dendritic solidification equations.

3.1 Wave Equation

Although the behavior of the numerical anisotropy was often reported in various one-dimensional optimizations of finite difference schemes, one of the first systematic attempts to specifically reduce the numerical anisotropy in finite difference schemes was introduced by Trefethen [44] in the framework of wave equation. To illustrate Trefethen’s approach, let us consider the two dimensional wave equation in the form

[TABLE]

defined in $R^{2}\times[0,\infty)$, with appropriate initial and boundary conditions. Using the Fourier-Laplace transform, it is ease to find the exact dispersion relation in the form $\omega^{2}=\xi^{2}+\eta^{2}$, where $\omega$ is the frequency and $(\xi,\eta)$ is the wavenumber vector. Equation (12) was discretized by Trefethen [44] on a Cartesian grid, using second order accurate schemes for both temporal and spatial derivatives as

[TABLE]

which was labeled $LF^{2}$. Then the same scheme was used to discretize equation (12), except the spatial derivatives were approximated along the diagonal directions with the space step $\sqrt{2}h$; this latter discretization was termed $LF^{2}$. It was found that the weighted averaging $2/3LF^{2}+1/3LF_{2}$ provided a low numerical anisotropy in the order of $(\sqrt{\xi^{2}+\eta^{2}}h)^{4}$. Slightly the same approach was used by Vichnevetsky [45] who corrected the numerical isotropy of the wave propagation in two dimensions using either the linear advection equation or the wave equation.

In a series of papers, Sescu et al. [30, 31, 32] proposed a technique to derive explicit multidimensional finite difference schemes for wave equation and Euler equations. By using the transformation matrix between two orthogonal reference frames, one aligned with the grid line and the other along the diagonal direction, the multidimensional finite difference scheme was obtained as

[TABLE]

where the multidimensional space shift operator $\textbf{E}_{x}^{\nu}\cdot u_{i,j}=u_{i+\nu,j}$ (see Vichnevetsky and Bowles [45] for one dimension) is used. The coefficients $a_{n}$ are those from the classical centered explicit schemes. The operator $\textbf{D}_{x}^{\nu}\cdot$ was defined as $\textbf{D}_{x}^{\nu}\cdot=\left(\textbf{E}_{x}^{\nu}\textbf{E}_{y}^{\nu}+\textbf{E}_{x}^{-\nu}\textbf{E}_{y}^{\nu}\right)\cdot$ The parameter $\beta$ is called isotropy corrector factor (ICF). The application of the Fourier transform to the multidimensional schemes gives the numerical wavenumber

[TABLE]

Then the numerical dispersion relation corresponding to two-dimensional wave equation was considered in the form \omega^{2}-\big{[}(\xi h)_{opt}^{*\hskip 5.69054pt2}+(\eta h)_{opt}^{*\hskip 5.69054pt2}\big{]}=0, and the ICF was determined by minimizing the integrated error between the phase or group velocities defined along $x$ and $x=y$ directions. Two curves in wavenumber-frequency space were considered: one was the intersection between the numerical dispersion relation surface and $\eta=0$ plane, and the other was the intersection between the numerical dispersion relation surface and the $\xi=\eta$ plane. These two curves were superposed in the $(Kh,\omega)$ plane, where Kh=\big{[}(\xi h)^{2}+(\eta h)^{2}\big{]}^{\frac{1}{2}}. Assuming that the equations of the two curves in $(Kh,\omega)$ plane are $\omega_{1}=\omega_{1}(Kh,\beta)$ and $\omega_{2}=\omega_{2}(Kh,\beta)$, the integrated error between the phase velocities was then calculated on a specified interval as C(\beta)=\int_{0}^{\eta}\big{|}c_{1}(Kh,\beta)-c_{2}(Kh,\beta)\big{|}^{2}d(Kh), where $c_{1}(Kh,\beta)$ and $c_{2}(Kh,\beta)$ are the numerical phase velocities. The minimization was done by equating the first derivative of $C(\beta)$ or $G(\beta)$ with zero, which provided the value of ICF, $\beta$.

Sescu et al. [33, 34] conducted a comprehensive stability analysis of the multidimensional schemes combined with either linear-multistep or multi-stage time marching schemes, and obtained several noteworthy results. For the Leap-Frog scheme applied to the advection equations, it was shown that the stability restriction corresponding to multidimensional schemes differs from the corresponding stability restriction via conventional schemes by the factor $(2\beta+2)/(\beta+2)$, where $\beta$ is the isotropy corrector factor. The conclusion was that the stability restrictions corresponding to multidimensional schemes are more convenient compared to the conventional schemes. For an arbitrary direction of the convection velocity with $|c_{x}|\geq|c_{y}|$, the stability restriction for conventional stencils was given by $\sigma_{x}+\sigma_{y}\leq CFL$, where $\sigma_{x}=k|c_{x}|/h$ and $\sigma_{y}=k|c_{y}|/h$. For multidimensional stencils the stability restriction was given by $(1+\beta)\sigma_{x}+\sigma_{y}\leq CFL(1+\beta)$ (where, for example, $CFL$ is $1$, $0.72874$ or $0.63052$ corresponding to E2, E4 or E6 scheme, respectively). Adams-Bashforth and Runge-Kutta time marching schemes in combination with conventional and multidimensional schemes were also analyzed, and it was found that the multidimensional schemes provide less restrictive stability limits.

3.2 Helmholtz Equation

Tam and Webb [42] performed an anisotropy correction of the finite difference representation of the Helmholtz equation,

[TABLE]

where $p$ is the pressure perturbation, $\nabla^{2}$ is the Laplacian operator, $f$ is the source distribution (e.g., a monopole), $\xi=2\pi/\lambda$ is the wavenumber, and $\lambda$ is the acoustic wavelength. Tam and Webb [42] showed that the finite difference discretization of the Helmholtz equation,

[TABLE]

with five grid points per wavelength introduces significant numerical anisotropy (equally-spaced grid is assumed in both x- and y-directions, and the spatial step is denoted as before by $h$). They constructed an anisotropy correction factor using asymptotic solutions to the continuous equation (16) and its finite difference approximation (17) as

[TABLE]

and

[TABLE]

respectively, where $(r_{ij},\theta_{ij})$ are polar coordinates, $K_{ij}=\alpha_{s}(\theta_{ij})\cos{\theta_{ij}}+\beta_{s}(\theta_{ij})\sin{\theta_{ij}}$ (with $\alpha_{s}$ and $\beta_{s}$ being the wavenumber components from the Fourier transform), and $G_{0}(\theta_{ij})$ and $G_{1}(\theta_{ij})$ are functions depending on $\alpha_{s}$, $\beta_{s}$, $\theta$ and the Fourier transform $\bar{F}$ of the source term (for more details see equations (19) and (21) in Tam and Webb [42]). The anisotropy corrector factor was then defined by the ratio between the absolute values of the two,

[TABLE]

The correction factor is independent of the distribution of sources, meaning that it can be computed once and for all types of sources. Significant reduction of the anisotropy error was obtained.

3.3 Advection Equation

Gaitonde and Shang [8] proposed a class of high-order compact difference-based finite-volume schemes which minimized the dispersion and isotropy error functions for the range of wavenumbers of interest. The starting point was the one dimensional advection equation,

[TABLE]

which was discretized using a finite volume approach as

[TABLE]

where $\bar{u}$ is the average value of $u$ inside a cell, $\bar{u}=1/h\int_{x_{i-1/2}}^{x_{i+1/2}}udx$, and $\bar{f}$ is the flux function approximating $f$, which is dependent on the values of $\bar{u}$ from neighbor cells. The reconstruction can be done by considering a primitive function $v=\int_{0}^{x}$ which must be discretized at the cell interface. Gaitonde and Shang [8] considered a five-point compact stencil in the form

[TABLE]

where $\alpha$, $a$, and $b$ are constants which determine the order of accuracy of the scheme. Using Taylor series expansions, they sacrificed the order of accuracy of the schemes by writing $a$ and $b$ as functions of $\alpha$,

[TABLE]

The spectral function associated with the scheme (23) is given as

[TABLE]

where $w=2\pi\xi h/L$ is the scaled wave number. The dispersion error is associated with the imaginary part of the spectral function, $w_{d}(w)=Im(\hat{A}(w))$. A scaled isotropy wavenumber was defined as

[TABLE]

where $\theta$ is the angle that the direction of propagation makes with the x-axis. An isotropy error function was defined by Gaitonde and Shang [8] in the form

[TABLE]

which was minimized to find the value of $\alpha_{opt}$ that gives the lowest numerical anisotropy. Numerical examples confirmed a considerable reduction of the isotropy error.

Sescu and Hixon [35, 36] extended the previous optimization performed in [31] to prefactored compact finite difference schemes [10, 11] applied to the advection equation. The prefactored compact schemes are defined on a three-point stencil and can return up to eight order of accuracy (see equations (2)). They can be used within a predictor-corrector type time marching scheme framework (MacCormack [26]), because the numerical derivatives are determined by sweeping from one boundary to the other, in both directions. Following the same analysis as in the case of explicit schemes, the multidimensional prefactored compact schemes were obtained as

[TABLE]

for fourth order of accuracy, and

[TABLE]

for sixth order of accuracy. $\beta$ is the isotropy corrector factor (ICF) and its magnitude can be determined by minimizing the dispersion error corresponding to the wave-front propagating along a grid line and the wave-front propagating along a diagonal direction.

Using Fourier analysis, the numerical wavenumbers and the numerical dispersion relation corresponding to the two dimensional wave equation were found. The individual (forward or backward) numerical wavenumber has both real and imaginary parts: the real part of the forward operator is equal to the real part of the backward operator, and the imaginary parts are opposite. As a result, in a MacCormack predictor-corrector scheme the overall imaginary part is zero. The real parts of the numerical wavenumbers corresponding to multidimensional schemes, for derivatives along $x$-direction, were given by:

[TABLE]

where $m=4$ for fourth and $m=6$ for sixth order of accuracy, $f_{4}(\eta_{x})=3\sin{\eta_{x}}/(2+\cos{\eta_{x}})$, $f_{6}(\eta_{x})=(28\sin{\eta_{x}}+\sin{2\eta_{x}})/(18+12\cos{\eta_{x}})$, $\eta_{x}=\xi h$, $\eta_{y}=\eta h$ and $\xi$ and $\eta$ are the components of the wavenumber.

In terms of numerical stability, more efficient stability restrictions were obtained as in the case of multidimensional explicit schemes. For example, multidimensional MacCormack schemes were found to provide a stability restriction in the form

[TABLE]

if $|c_{x}|\geq|c_{y}|$, and

[TABLE]

if $|c_{y}|\geq|c_{x}|$. For diagonal directions, with respect to the grid, ($|c_{x}|=|c_{y}|=|c|$) the stability restriction becomes

[TABLE]

It is obvious that the right hand side of equation (35) is greater than $1/(2\xi_{max})^{3/2}$ when $\beta>0$, and goes to $1/(\xi_{max})^{3/2}$ when $\beta\rightarrow\infty$. This generated more efficient stability restrictions by using multidimensional compact schemes. Test cases showed that the multidimensional compact schemes were more efficient for both fourth and sixth order accurate schemes.

3.4 Maxwell Equations

Sun and Trueman [40] performed an optimization of finite difference schemes applied to Maxwell equations, in terms of reducing the dispersion and isotropy errors. For brevity, we show here the numerical dispersion relations (for finite differencing representations of the Maxwell equations, see equations (1), (2) and (4) in Sun and Trueman [40]):

[TABLE]

corresponding to a grid line, and

[TABLE]

corresponding to the diagonal direction, where $w$ is a weighting factor, $\beta_{a}$ is the numerical phase constant along the grid line, $\beta_{d}$ is the numerical phase constant along the diagonal direction, $\omega$ is the frequency, and $k$ is the time step (an equally-spaced grid is considered again). The optimization in terms of reducing the numerical anisotropy was done by eliminating the time step terms in equations (36) and (37) to obtain

[TABLE]

This optimal weight $w_{i}$ is a function of mesh density only, and is not dependent of the time step size or the frequency of the signal. This method theoretically provides a uniform phase velocity in all directions. Further optimizations of this scheme were performed in another paper of Sun and Trueman [41].

Koh et al. [19] derived a two dimensional finite-difference time-domain method, discretizing the Maxwell equations, to eliminate the numerical dispersion and anisotropy. The proposed scheme is given as

[TABLE]

where $d_{t}^{2}$ is the central difference operator with respect to time,

[TABLE]

with $p$ or $q$ being either $x$ or $y$, and

[TABLE]

where $f$ is a generic function. In equation (3.4), $E$ is the electric field, $H$ is the magnetic field strength, $\sigma$, $\mu$ and $\epsilon$ are the conductivity, permeability and the permittivity, respectively, of the domain, $k$ is the time step, and $h$ is the spatial step in all directions. For a nonconductive media $\sigma=0$, the numerical dispersion relation of can be obtained as

[TABLE]

where $C_{+}=\sin^{2}(\xi h/2)+\sin^{2}(\eta h/2)$, $C_{\times}=\sin^{2}(\xi h/2)\sin^{2}(\eta h/2)$, and $\xi$ and $\eta$ are the components of the wavenumber. Equation (43) is a quadratic equation in $\alpha$, and the solution is given as

[TABLE]

An optimal value for $\alpha$, achieving an isotropic numerical phase velocity, can be simply estimated as the mean value of $\alpha$ over the azimuthal angles, and it was found that it remains constant (approximately, $0.167$) for a wide range of grid sizes, and it is insensitive to the value of the Courant number.

Kim et al. [17] derived new three-dimensional isotropic dispersion-finite-difference time-domain schemes (ID-FDTD) based on a a linear combination of the traditional central difference equation and a new difference equation based on the extra sampling points. They used the same scaling factors as for the two-dimensional case to attain isotropic dispersion and exact phase velocity. Based on the weighting factors, seven different FDTD schemes were formulated, including the Yee scheme [46]. Among the seven proposed FDTD schemes, three showed improved isotropy of the dispersion compared to the dispersion of the Yee scheme. For the sake of brevity, the complete expressions of the schemes are not included here (see Kim et al. [17] for more details), and only the numerical dispersion relation is briefly presented. Plane wave solutions were introduced in discretized forms as

[TABLE]

[TABLE]

where $I=\sqrt{-1}$, $\omega$ is the frequency, $(\xi,\eta,\zeta)$ is the numerical wavenumber vector, and $\textbf{E}_{0}$ and $\textbf{H}_{0}$ are constant vectors. After inserting (44) and (45) into the discretized form of the Maxwell equations (see equation (10) in Kim et al. [17]), matrix equations are obtained as $C\textbf{H}_{0}=S_{t}\epsilon_{0}\textbf{E}_{0},\hskip 2.84526ptC\textbf{E}_{0}=S_{i}\mu_{0}\textbf{H}_{0}$ where

[TABLE]

and $K_{p}=S_{p}/h[\alpha(P_{p}-Q_{p})-\beta Q_{p}/2+1]$ ($p$ being either $x$, $y$ or $z$), $S_{x}=\sin(\xi h/2)$, $S_{y}=\sin(\eta h/2)$, $S_{z}=\sin(\zeta h/2)$, $P_{x}=SySz$, $P_{y}=SxSz$, $P_{z}=SxSy$, $Q_{x}=S_{y}^{2}+S_{z}^{2}$, $Q_{y}=S_{x}^{2}+S_{z}^{2}$, $Q_{z}=S_{x}^{2}+S_{y}^{2}$, and $S_{t}=\sin{\omega k/2}/k$. An eigenvalue equation was obtain as

[TABLE]

and the numerical dispersion relation was obtained by vanishing the associated determinant,

[TABLE]

where $c_{0}=1/\sqrt{\epsilon_{0}\mu_{0}}$. The isotropy correction was performed by defining the values of the weighting factors $\alpha$ and $\beta$, which unlike the two-dimensional case are not unique. Kim et al. [17] used the scaling factor from the two-dimensional case, and modified the numerical dispersion relation to estimate the weighting factors.

3.5 Dendritic Solidification

Kumar [21] derived isotropic finite difference schemes for the first and second derivatives in the context of symmetric dendritic solidification. The first derivative was discretized as

[TABLE]

which involves grid points not only along $x$-direction, but also along $y$-direction. The Taylor expansion of the scheme (3.5) can be written as $(\partial_{x}u)_{I,i,j}=(1+h^{2}/6\nabla^{2})(\partial_{x}u)_{i,j}$, where the leading order term involves the Laplacian only, implying no directional dependence. The second derivative was discretized as

[TABLE]

where the Taylor expansion is given by $(\partial_{xx}u)_{I,i,j}=(1+h^{2}/12\nabla^{2})(\partial_{xx}u)_{i,j}$, being again a function of the Laplacian only. The conventional cross derivative $(\partial_{xy}u)_{I,i,j}$ was found to be intrinsically isotropic according to the criterion developed by Kumar [21]. The Laplacian can be obtained by combining the isotropic derivatives along x- and y-directions, $(\nabla^{2}u)_{i,j}=(\partial_{xx}u)_{I,i,j}+(\partial_{yy}u)_{I,i,j}$. Significant reduction of the numerical anisotropy was obtained by using these schemes. Shen and Cangellaris [37] exploited further this approach to develop new isotropic finite-difference time-domain schemes modeling electromagnetic wave propagation.

4 Concluding Remarks

Numerical anisotropy in finite difference discretizations of partial differential equations was discussed and reviewed. In some instances, the numerical anisotropy can be neglected, and the focus is directed toward other types of one-dimensional errors, such as numerical dispersion, dissipation or aliasing. These errors can be analyzed in the context of one dimensional differencec equations, while the extension to multidimensional discretizations is straightforward. By increasing the accuracy of one dimensional schemes or by increasing the number of grid points in the grid, the isotropic characteristics of the waves in multi-dimensions can be improved. These two practices, however, are not always effective since an increase in accuracy may require larger stencils which may introduce spurious waves at the boundaries of the domain, while by increasing of the resolution of the grid may increase the computational time. It is necessary then to analyzed the schemes in multi-dimensions and design specific optimizations with the specific objective of reducing the numerical anisotropy, and at the same time of conserving the dispersion characteristics of the corresponding one dimensional schemes. Various attempts to reduce the numerical anisotropy in finite differencing applied to various model equations were presented and discussed.

Future directions should focus on optimizations of existing compact finite difference schemes in terms of reducing the numerical anisotropy, or derivations of novel compact schemes with low numerical anisotropy. Optimizations and derivations of finite volume schemes (in terms of reducing the numerical anisotropy) applied to either structured or unstructured grids should be also taken into account, especially in the framework of wave propagation problems. Filtering schemes, as applied, for example, in large eddy simulations to separate the small scales from the large scales, may experience numerical anisotropy since they are effective at high wavenumber ranges. Optimizations of such filters in terms of reducing the numerical anisotropy is also another future area of research.