\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}
\usepackage{subfigure}

\AtBeginDocument{{\noindent\small Seventh Mississippi State - UAB
Conference on Differential Equations and Computational
Simulations, {\em Electronic Journal of Differential Equations},
Conf. 17 (2009),  pp. 213--225.\newline ISSN: 1072-6691. URL:
http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2009 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document} \setcounter{page}{213}
\title[\hfilneg EJDE-2009/Conf/17 \hfil Optimized difference schemes]
{Optimized difference schemes for multidimensional hyperbolic
partial differential equations}

\author[A. Sescu, A. A. Afjeh, R. Hixon\hfil EJDE/Conf/17 \hfilneg]
{Adrian Sescu, Abdollah A. Afjeh, Ray Hixon} % not in alphabetical order

\address{Adrian Sescu \newline
University of Toledo, Toledo, OH 43606, USA}
\email{asescu@utnet.utoledo.edu}

\address{Abdollah A. Afjeh \newline
University of Toledo, Toledo, OH 43606, USA}
\email{aafjeh@utnet.utoledo.edu}

\address{Ray Hixon \newline
University of Toledo, Toledo, OH 43606, USA}
\email{dhixon@utnet.utoledo.edu}


\thanks{Published April 15, 2009.}
\subjclass[2000]{76Q05, 65M06, 65M20, 65Q05}
\keywords{Computational aeroacoustics; isotropy error; difference methods}

\begin{abstract}
 In numerical solutions to hyperbolic partial differential equations 
 in multidimensions, in addition to  dispersion and dissipation errors, 
 there is a grid-related error (referred to as isotropy error or 
 numerical anisotropy) that  affects the directional dependence of 
 the wave propagation.
 Difference schemes are mostly analyzed and optimized in one
 dimension, wherein the anisotropy correction may not be effective enough.
 In this work, optimized multidimensional difference schemes with
 arbitrary order of accuracy are designed to have improved isotropy
 compared to conventional schemes. The derivation is performed based
 on Taylor series expansion and Fourier analysis. The schemes are
 restricted to equally-spaced Cartesian grids, so the generalized
 curvilinear transformation method and Cartesian grid methods are
 good candidates.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\allowdisplaybreaks

\section{Introduction}

The numerical anisotropy is a type of error occurring in the
numerical approximation of hyperbolic partial differential equations (PDE), using structured grids.  This error can be reduced by using, for example, high-resolution difference schemes or sufficiently dense grids. While the former possibility requires special care at the
boundaries, the latter might be computationally expensive. This
work proposes a way to derive difference schemes in
multidimensions that use points from more than one direction: the
result is a improved isotropy of the wave propagation.

The optimization of the centered spatial differencing schemes in
terms of lowering the dispersion error especially for
Computational Aeroacoustics, Large Eddy Simulations and Direct
Numerical Simulations is an actual field of research. Among
others,  the works of Lele \cite{Lel} and Tam and
Webb \cite{Tam93} are considered good starting points; the former
conducted optimizations of Pad\'{e} schemes using Fourier
analysis, and the latter used the so-called
dispersion-relation-preserving (DRP) concept to derive explicit
high-resolution finite difference stencils. Kim and Lee \cite{Kim}
performed an analytic optimization of the compact finite
difference schemes. They showed that an analytic optimization
produces the maximum spatial resolution characteristics of the
compact finite difference approximation in the evaluation of the
spatial finite derivatives. Li \cite{Li} has proposed new
wavenumber-extended high-order upwind-biased schemes up to
11th-order by means of Fourier analysis. He showed that both the
upwind-biased scheme of order 2N - 1 and the corresponding
centered differencing scheme of order 2N have the same dispersion
characteristics. Mahesh \cite{Mah} derived a family of compact
finite difference schemes for the spatial derivatives in the
Navier-Stokes equations based on Hermite interpolation. He
simultaneously solved for the first and second derivatives,
getting higher-order of accuracy and better spectral resolution.
Hixon \cite{Hix} derived prefactored high-order compact schemes
which use three-point stencils and return up to eighth-order
accuracy. His schemes combine the tridiagonal compact formulation
with the optimized split derivative operators of an explicit
MacCormack type scheme. The tridiagonal matrix inversion was
avoided by using bidiagonal matrices for the forward and backward
operators. The optimization of Hixon's \cite{Hix} schemes in terms
of dispersion error was performed by Ashcroft and Zhang \cite{Ash}
who used Fourier analysis to select the coefficients of the biased
operators, such that the dispersion characteristics match those of
the original centered compact scheme and their numerical
wavenumbers have equal and opposite imaginary components.

All of the above optimizations were performed in one-dimensional
space and they may suffer from the isotropy error in multidimensions.
An extended analysis of the isotropy error was performed
by Vichnevetsky \cite{Vic} who also solved the two-dimensional wave
equation using two different schemes for the Laplacian operator, and
averaged the two solutions. Considerable improvement of the isotropy
of wave propagation was obtained based on variation of the weighted
average. The same idea was considered by Trefethen \cite{Tre}
 who used the leap frog scheme to solve the wave equation in
two dimensions. Zingg and Lomax \cite{Zin} performed optimizations
of finite difference schemes applied to regular triangular grids,
that give six neighbor points for a given node. Tam and Webb
 \cite{Tam94} proposed an anisotropy correction for Helmholtz
equation; they found the anisotropy correction factor applicable to
all noise radiation problems, irrespective of the complexity of the
noise sources. Lin and Sheu \cite{Lin} used the idea of
DRP of Tam and Webb \cite{Tam93} 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 the
nine-point grid system, resulting in nine unknown coefficients. Eight
of them were determined by employing Taylor series expansions, and
the remaining one was determined by requiring that the
two-dimensional numerical dispersion relation is the same as the
exact dispersion relation.

In this paper, multidimensional optimized schemes are derived
using the weighted averaging technique and the transformation
matrix between two orthogonal bases. Because the optimized schemes
are linear combinations of classical finite difference schemes,
they have the same order of accuracy as the corresponding
classical schemes. Using Fourier analysis, their advantage is
revealed in terms of isotropy error: compared to classical
schemes, they have improved isotropy.

The organization of the paper is as follow. In Sec. II, the
definition of the isotropy error is presented. In Sec. III, the
dispersion relation for hyperbolic systems is determined.  In Sec.
IV, the procedure of deriving the optimized schemes is presented.
In Sec. V the isotropy corrected factor is determined, and in Sec.
VI, results for a problem from the First Computational
Aeroacoustics Workshop \cite{Har} are reported. Concluding remarks
are given in Sec. VII and the Taylor series expansions for the
second and fourth order optimized schemes are written in the
appendix.

\section{Isotropy error}

Wave propagation is an inherent feature of the solutions of
hyperbolic partial differential equations. In multidimensional
space most of the waves or wave packets are propagating in all
directions with the same phase or group velocity, respectively. A
type of error occurring in the numerical approximation of
hyperbolic equations in multidimensions is the numerical
anisotropy. In addition to being dependent upon frequency, the
numerical  phase or group velocity is also dependent upon
direction. The easiest way to illustrate this is by considering
the advection equation in two dimensions:
\begin{gather}\label{1}
\frac {\partial{u}}{\partial{t}} =\textbf{c} \nabla u; \\
\label{2}
u(\textbf{r},0)=u_0(\textbf{r})
\end{gather}
where $\textbf{r}=(x,y)$ is the vector of spatial coordinates,
$\textbf{c}=c(\cos\alpha   \sin\alpha)$ is the velocity vector
($c$ is scalar), $\nabla=(\partial/\partial{x}
\partial/\partial{y})^T$ and $u(x,y,t)$ and $u_0(\textbf{r})$ are
scalar functions.

A simple semi-discretization of the equation (\ref{1}) is obtained
with Cartesian coordinates on a square grid,
\begin{equation}\label{3}
    \frac{du}{dt}=-\frac{c}{2 h} \big[\cos \alpha
(u_{{i+1},j}-u_{{i-1},j})+ \sin \alpha
(u_{i,{j+1}}-u_{i,{j-1}})\big]
\end{equation}
where $h$ is the grid step (the same in both directions). Consider
the Fourier-Laplace transform:
\begin{equation}\label{4}
    \tilde{u}(k_1,k_2, \omega)=\frac{1}{(2\pi)^{3}}\int_{0}^{\infty} \int
\int_{-\infty}^{\infty}u(x,y,t) e^{-i (kx \cos\alpha +ky
\sin\alpha - \omega t)} dx dy dt
\end{equation}
and its inverse
\begin{equation}\label{5}
    u(x,y,t)=\int_\Gamma \int \int_{-\infty}^{\infty}
\tilde{u}(k_1, k_2,
    \omega) e^{i(kx\cos\alpha +ky \sin\alpha
    - \omega t)}dk_1dk_2 d\omega,
\end{equation}
where $k \cos \alpha$ and $k \sin \alpha$ are the components of the
wave number and $\omega$ is the frequency. $\Gamma$ is a line
parallel to the real axis in the complex $\omega$-plane above all
poles and singularities of the integrand (Tam and Webb,
 \cite{Tam93}). The application of Fourier-Laplace transform to Eq.
(\ref{1}) gives the exact dispersion relation:
\begin{equation}\label{6}
    \omega=ck(\cos^2 \alpha+\sin^2 \alpha)=ck
\end{equation}
such that the phase velocity is obtained as
\begin{equation}\label{7}
    c_{e}=\frac{\omega}{k}=c
\end{equation}
Plugging (\ref{6}) back into (\ref{4}), $u(x,y,t)$ is obtained as a
superposition of sinusoidal solutions in the plane with constant
phase lines given by
\begin{equation}\label{8}
    x\cos\alpha+y\sin\alpha-c_{e}t=\textrm{const.}
\end{equation}
and, according to Eq. (\ref{7}), the phase velocity, $c_{e}$,
does not depend on direction (it is isotropic).

If the same is done for the numerical approximation,  (\ref{2}),
the numerical dispersion relation takes the form:
\begin{equation}\label{9}
    \omega=\frac{c}{h}\big[\cos\alpha \sin(kh\cos\alpha)
    +\sin\alpha \sin(kh\cos\alpha)\big]
\end{equation}
and the numerical phase velocity will be given by:
\begin{equation}\label{10}
    c_{n}=\frac{\omega}{k}=\frac{c}{kh}\big[\cos\alpha \sin(kh\cos\alpha)
+\sin\alpha \sin(kh\cos\alpha)\big]
\end{equation}
The constant phase lines are expressed by the equation
\begin{equation}\label{11}
    x\cos\alpha+y\sin\alpha-c_{n}t=const.
\end{equation}
and moves with the phase velocity $c_{n}$. The numerical
anisotropy is the effect of the numerical phase velocity which is
dependent on direction.

The same considerations are valid for the group velocity defined as
\begin{equation}\label{12}
    g_{e}=\frac{\partial\omega}{\partial k}=c
\end{equation}
showing that the exact group velocity is the same as the phase
velocity because the dispersion relation is a linear function of
$k$. But the numerical group velocity is different from the
numerical phase velocity,
\begin{equation}\label{13}
    g_{n}=\frac{\partial\omega}{\partial k}=c\big[\cos^2
\alpha\cos(kh\cos\alpha)+\sin^2 \alpha\cos(kh\sin\alpha )\big],
\end{equation}
which is also dependent on direction. This directional dependence
(in both phase and group velocities) produces the numerical anisotropy.

\section{Dispersion Relations for First Order Hyperbolic PDEs}

Consider a hyperbolic set of first order partial differential equations
defined in multidimensional space. Let $\Omega$  be an open subset in
$\mathbf{R}^p$, and let $\mathbf{f}_j$, $1\leq j \leq d$ be $d$
smooth functions from $\Omega$ to $\mathbf{R}^p$. The general form
is given by
\begin{equation}\label{14}
    \frac{\partial \textbf{u}}{\partial t}
 +\sum_{j=1}^{d}\frac{\partial}{\partial
    x_j}\textbf{f}_j(\textbf{u})=0
\end{equation}
where
\begin{equation}\label{15}
    \textbf{u}=(u_1 \dots u_p)^T
\end{equation}
is a vector valued function and
\begin{equation}\label{16}
    \mathbf{f}_j=(f_{1j}\dots f_{pj})^T
\end{equation}
is called flux vector. To simplify the analysis, suppose that the flux vector is a linear function of $\mathbf{u}$.
Equation (\ref{14}) written in primitive form is:
\begin{equation}\label{17}
    \frac{\partial u_i}{\partial t}+\sum_{j=1}^{d}\Big[\frac{\partial f_{ij}(u_i)}{\partial
    u_k}\Big]\frac{\partial u_i}{\partial x_j}=0, \quad 1\leq i,k
    \leq p
\end{equation}
where the Einstein summation convention is used. The set of
equations (\ref{14}) is said to be hyperbolic if the eigenvalues
associated with the matrix
\begin{equation}\label{18}
    \mathbf{B}(\mathbf{u},\mathbf{\alpha})=\sum_{j=1}^d \alpha _j \Big[\frac{\partial f_{ij}(u_i)}{\partial
    u_k}\Big]_{1\leq i,k\leq p}
\end{equation}
are all real and the associated eigenvectors are linearly independent. The hyperbolicity implies that the initial set of
equations admits wave-form solutions. The Fourier-Laplace transform to Eq. (\ref{16}) yields a matrix equation:
\begin{equation}\label{19}
    \mathbf{A} \hat{\mathbf{u}}=\tilde{\mathbf{G}}
\end{equation}
where $\tilde{\mathbf{G}}$ may result from the initial conditions.
The dispersion relations associated with
different waves are determined by making the determinant of the
matrix $\mathbf{A}$ zero. This occurs when any of its eigenvalues is
zero. If $(\lambda_i)_{q<i\leq m}$ are the $m-q+1$ eigenvalues
corresponding to a type of wave, then the dispersion relation
associated with this wave is defined as
\begin{equation}\label{20}
    \prod_{i=q}^m \lambda_i =0
\end{equation}
All properties associated with that wave are encoded in the dispersion
relation \cite{Whi}, and the specific errors (numerical dispersion,
dissipation or anisotropy) that result from a particular numerical
approximation are most commonly analyzed by comparing the numerical
dispersion relation with the exact dispersion relation.


\subsection*{Example}

Consider waves in a two-dimensional uniform,
isentropic, subsonic compressible fluid flow $(\bar{u},0)$ along the $x$ direction. For small perturbations in the density and the
velocity components we may linearize and nondimensionalize the Euler
equations of gas-dynamics, so that we the Linearized Euler
Equations (LEE) are obtained:
\begin{equation}\label{21}
    \frac{\partial\mathbf{Q}}{\partial t}
+\frac{\partial\mathbf{E}}{\partial x}+\frac{\partial\mathbf{F}}{\partial
    y}=0
\end{equation}
where
\begin{equation}\label{22}
    \mathbf{Q}= \begin{pmatrix}
                   \rho' \\
                   u' \\
                   v' \\
                   p'
                 \end{pmatrix},
\quad \mathbf{E}=\begin{pmatrix}
                   M \rho'+u' \\
                   M u'+p' \\
                   M v' \\
                   M p'+u'
                 \end{pmatrix},
\quad \mathbf{F}=\begin{pmatrix}
                   v' \\
                   0 \\
                   p' \\
                   v'
                 \end{pmatrix}
\end{equation}
$\rho'$, $u'$, $v'$ and $p'$ are the perturbations of
density, $x$-component velocity, $y$-component velocity and pressure,
respectively. $M$  is the Mach number corresponding to the mean velocity.

The application of the Fourier-Laplace transform to
(\ref{21}) will generate the matrix
\begin{equation}\label{23}
    \mathbf{A}=\begin{pmatrix}
                         \omega - k_1 M & -k_1 & -k_2 & 0 \\
                         0 & \omega - k_1 M & 0 & -k_1 \\
                         0 & 0 & \omega - k_1 M & -k_2 \\
                         0 & -k_1 & -k_2 & \omega - k_1 M
                       \end{pmatrix}
\end{equation}
where $k_1$  and $k_2$  are the components of the wave number and
$\omega$ is the frequency. It is easy to show that the eigenvalues
of matrix $\mathbf{A}$ are
\begin{equation}\label{24}
   \begin{gathered}
  \lambda_1=\lambda_2=\omega -k_1M \\
  \lambda_3=(\omega -k_1M)+\sqrt{k_1^2+k_1} \\
  \lambda_4=(\omega -k_1M)-\sqrt{k_1^2+k_1}
\end{gathered}
\end{equation}

The first and the second eigenvalues ($\lambda_1$  and $\lambda_2$)
correspond to entropy and vorticity waves, respectively. The third
and the fourth ($\lambda_3$ and $\lambda_4$) correspond to acoustic
waves, and the dispersion relation is given by
\begin{equation}\label{25}
    \lambda_3\lambda_4=(\omega -k_1M)^2-(k_1^2+k_2^2)=0
\end{equation}
For a stationary mean flow ($M=0$) the dispersion relation (\ref{25})
becomes
\begin{equation}\label{26}
    \omega^2-(k_1^2+k_2^2)=0
\end{equation}

\section{The Derivation of the Optimized Schemes}

An equally-spaced, two-dimensional Cartesian grid is considered with
$i$ index on the $x$-direction and  $j$ index on the $y$-direction;
the grid step is denoted by $h$. This kind of grid is typical
to generalized curvilinear transormation methods, wherein the
physical domain is mapped into the computational domain.
Two orthogonal bases are considered, one
$(xOy)$ related to Cartesian grid directions and the other
$(x'Oy')$ positioned at  $45 ^\circ$ with respect to the first. The
transformation matrix between these two orthogonal bases is used to
derive the optimized schemes:
\begin{equation}\label{27}
   \begin{pmatrix}
     x' \\
     y'
   \end{pmatrix}
=\begin{pmatrix}
             \cos \alpha & \sin \alpha \\
             -\sin \alpha & \cos \alpha
  \end{pmatrix}
  \begin{pmatrix}
    x \\
    y
  \end{pmatrix},
\end{equation}
where $\alpha$  is the angle between $x$ and $x'$ axes ($45 ^\circ$
in this case). The relation between the derivatives of a function
$(x,y)$ based on  (\ref{27}) is given by:
\begin{equation}\label{28}
      \begin{pmatrix}
        \frac{\partial u}{\partial x} \\
        \frac{\partial u}{\partial y}
      \end{pmatrix}
 =\begin{pmatrix}
                      \cos \alpha & -\sin \alpha \\
                      \sin \alpha & \cos \alpha
   \end{pmatrix}
   \begin{pmatrix}
     \frac{\partial u}{\partial x'} \\
     \frac{\partial u}{\partial y'}
   \end{pmatrix}.
\end{equation}
The general forms of classical difference schemes for the $x$- and
$y$- derivatives are written as:
\begin{equation}\label{29}
\Big( \frac {\partial{v}}{\partial{x}} \Big)_{i,j}
=
\sum_{\nu=-M}^{\nu=M} a_\nu
\textbf{E}_{x}^{\nu} \cdot v_{i,j}
\end{equation}
and
\begin{equation}\label{30}
\Big( \frac {\partial{v}}{\partial{y}} \Big)_{i,j}
=
\sum_{\nu=-M}^{\nu=M} a_\nu
\textbf{E}_{y}^{\nu}\cdot v_{i,j}
\end{equation}
where multidimensional space shift operators
(see Vichnevetsky and Bowles \cite{Vic} for one dimension) are used:
\begin{equation}\label{31}
\textbf{E}_{x}^{\nu}\cdot v_{i,j} = v_{i+\nu,j}; \quad
\textbf{E}_{y}^{\nu}\cdot v_{i,j} = v_{i,j+\nu}.
\end{equation}
The $a_n$  coefficients are given in Table \ref{table1} for several
frequently used centered schemes. By using the transformation matrix
and a weighted averaging we get the optimized schemes in the form:
\begin{equation}\label{32}
\Big( \frac {\partial{v}}{\partial{x}} \Big)_{i,j}
=
\frac{1}{h(1+\beta)}
\sum_{\nu=-M}^{\nu=M} a_\nu
\Big(
\textbf{E}_{x}^{\nu} +
\frac{\beta}{2}
\textbf{D}_{x}
\Big)\cdot v_{i,j}
\end{equation}
and
\begin{equation}\label{33}
\Big( \frac {\partial{v}}{\partial{y}} \Big)_{i,j}
=
\frac{1}{h(1+\beta)}
\sum_{\nu=-M}^{\nu=M} a_\nu
\Big( \textbf{E}_{y}^{\nu} +
\frac{\beta}{2}
\textbf{D}_{y}
\Big)\cdot v_{i,j}
\end{equation}
where the operators $\textbf{D}_{x}^{\nu}\cdot$ and
$\textbf{D}_{y}^{\nu}\cdot$ are defined as
\begin{equation}\label{34}
\textbf{D}_{x}^{\nu}\cdot =
\left(
\textbf{E}_{x}^{\nu}\textbf{E}_{y}^{\nu} +
\textbf{E}_{x}^{-\nu}\textbf{E}_{y}^{\nu}
\right)\cdot;
\quad
\textbf{D}_{y}^{\nu}\cdot =
\left(
\textbf{E}_{x}^{\nu}\textbf{E}_{y}^{\nu} +
\textbf{E}_{x}^{\nu}\textbf{E}_{y}^{-\nu}
\right)\cdot
\end{equation}
The parameter $\beta$ is called
isotropy corrector factor (ICF).

\begin{table}[ht]
\caption{Coefficients of various explicit centered
finite difference schemes.}
\begin{center}
\begin{tabular}{|lrrr|}
  \hline
  Central scheme & $a_1=-a_{-1}$ & $a_2=-a_{-2}$ & $a_3=-a_{-3}$ \\
  \hline
  $2^{nd}$ order centered & 1/2 & 0 & 0 \\
  \hline
  $4^{th}$ order centered & 2/3 & -1/12 & 0 \\
  \hline
  $6^{th}$ order centered & 3/4 & -3/20 & 1/60 \\
  \hline
  $4^{th}$ order DRP & 0.7708824 & -0.1667059 & 0.0208431 \\
  \hline
\end{tabular}
\end{center}
\label{table1}
\end{table}

The application of the Fourier transform to
the schemes gives the numerical wavenumbers
\begin{equation}\label{35}
    (k_1 h)_c^*=\sum_{n={-N}}^M a_n e^{nIk_1 h}, \quad %\label{36}
    (k_2 h)_c^*=\sum_{n={-N}}^M a_n e^{nIk_2 h},
\end{equation}
for classical schemes and
\begin{gather}\label{37}
    (k_1 h)_{opt}^*=\frac{2}{(1+\beta)} \sum _{n=-N}^M a_n \Big\{e^{nIk_1 h}+\frac{\beta}{2}\big[e^{nI(k_1+k_2) h}+e^{nI(k_1-k_2)
    h}\big]\Big\}, \\
 \label{38}
    (k_2 h)_{opt}^*=\frac{2}{(1+\beta)} \sum _{n=-N}^M a_n \Big\{e^{nIk_2 h}+\frac{\beta}{2}\big[e^{nI(k_1 \Delta x+k_2 \Delta y) h}-e^{nI(k_1-k_2)
    h}\big]\Big\},
\end{gather}
for the optimized schemes ($I=\sqrt{-1}$) .

According to (\ref{26}), and assuming that the time
integration is free of numerical dissipation and dispersion, the
numerical dispersion relation corresponding to two-dimensional wave
equation is:
\begin{equation}\label{41}
    \omega^2-\big[(k_1h)_{opt}^{* \hspace{2 mm} 2} +(k_2h)_{opt}^{* \hspace{2 mm} 2}\big]=0.
\end{equation}
which will be used in determining ICF (next section).

\section{ICF Calculation}

ICF is found by minimizing the integrated error between the phase or
group velocities on $x$ and $x=y$ directions. Let's consider two
curves in wavenumber-frequency space: one of them is the
intersection between the numerical dispersion relation surface and
$k_2=0$ plane, and the other is the intersection between the
numerical dispersion relation surface and the $k_1=k_2$ plane. These
two curves are superposed in the $(kh, \omega)$  plane, where
\begin{equation}\label{42}
    k h=\big[(k_1 h)^2+(k_2 h)^2\big]^{1/2}
\end{equation}
Suppose that the equations of the two curves in $(k \Delta x,\omega)$
plane are
\begin{equation}\label{43}
    \omega_1=\omega_1 (kh,\beta), \quad %\label{44}
    \omega_2=\omega_2 (kh,\beta).
\end{equation}
The integrated error between the phase velocities is then calculated
on a specified interval,
\begin{equation}\label{45}
    C(\beta)=\int_0^ \eta\big|c_1(kh,\beta)-c_2(kh,\beta)\big|^2d(kh),
\end{equation}
where
\begin{equation}\label{46}
    c_1(kh,\beta)=\frac{\omega_1(kh,\beta)}{kh}, \quad %\label{47}
   c_2(kh,\beta)=\frac{\omega_2(kh,\beta)}{kh}
\end{equation}
are the numerical phase velocities. The upper limit of the
integration, $\eta$, is a real number between $0$ and $\pi$, and
depends on the type of the scheme being optimized. The optimization
can be also performed using the integrated error between the numerical group
velocities,
\begin{equation}\label{48}
    G(\beta)=\int_0^ \eta \big|g_1(kh,\beta)-g_2(kh,\beta)\big|^2d(kh)
\end{equation}
where
\begin{equation}\label{49}
    g_1(kh,\beta)=\frac{\partial \omega_1(kh,\beta)}{\partial (kh)},\quad %\label{50}
    g_2(k \Delta x,\beta)=\frac{\partial \omega_2(k \Delta x,\beta)}{\partial (k \Delta x)}
\end{equation}
are the group velocities.


\begin{figure}[htp]
  \begin{center}
\includegraphics[width=0.23\textwidth]{fig1a} 
\includegraphics[width=0.23\textwidth]{fig1b} 
\includegraphics[width=0.23\textwidth]{fig1c}
\includegraphics[width=0.23\textwidth]{fig1d} 
  \end{center}
  \caption{Polar diagram of normalized phase velocities (a and b)
   and group velocities (c and d) as a function of points per wavelength (PPW)
   and the direction of propagation: (a) and (c) using second-order classical
   schemes; (b) and (d) using second-order optimized schemes.}
  \label{Fig1}
\end{figure}


\begin{figure}[htp]
  \begin{center}
\includegraphics[width=0.23\textwidth]{fig2a} 
\includegraphics[width=0.23\textwidth]{fig2b} 
\includegraphics[width=0.23\textwidth]{fig2c} 
\includegraphics[width=0.23\textwidth]{fig2d} 
  \end{center}
  \caption{Polar diagram of normalized phase velocities (a and b) and group velocities (c and d) as a
   function of points per wavelength (PPW) and the direction of propagation:
   (a) and (c) using fourth-order classical schemes; (b) and (d)
   using fourth-order optimized schemes.}
  \label{Fig2}
\end{figure}

The minimization is done by equalizing the
first derivative of $C(\beta)$ or  $G(\beta)$ with zero:
\begin{equation}\label{51}
    \frac{dC(\beta)}{d\beta}=0 \quad  or  \quad \frac{dG(\beta)}{d\beta}=0
\end{equation}
which gives the value of ICF, $\beta$ . Corresponding to dispersion
relation of the two-dimensional wave equation, for the second-order
centered schemes $\beta\cong0.53$, for the fourth-order centered
schemes $\beta\cong 0.282$, and for the sixth-order centered schemes
$\beta\cong0.152$. In Figures \ref{Fig2} and  \ref{Fig3} polar
diagrams of normalized phase or group velocities for second and fourth order
centered classical and optimized schemes are shown. The diagrams are
plotted for different numbers of points per wavelength (PPW).

\section{Numerical Tests}

A two-dimensional problem from First Computational
Aeroacoustics Workshop \cite{Har} is considered first. An acoustic wave
initially situated in $O(0,0)$ point and an entropy wave in
$P(67,67)$ point are convected with the mean flow along the $x =
y$ direction. The linearized two-dimensional
Euler equations are considered:
\begin{equation}\label{52}
    \frac{\partial \mathbf{Q}}{\partial t}
+\frac{\partial \mathbf{E}}{\partial x}
+\frac{\partial \mathbf{F}}{\partial y}=0
\end{equation}
where
\begin{equation}\label{53}
    \mathbf{Q}=\begin{pmatrix}
                   \rho' \\
                   u' \\
                   v' \\
                   p'
                 \end{pmatrix},\quad
  \mathbf{E}=\begin{pmatrix}
                   M_x \rho'+u' \\
                   M_x u'+p' \\
                   M_x v' \\
                   M_x p'+u'
                 \end{pmatrix}, \quad
  \mathbf{F}=\begin{pmatrix}
                   M_y \rho'+v' \\
                   M_y u' \\
                   M_y v'+p' \\
                   M_y p'+v'
                 \end{pmatrix}
\end{equation}
$M_x$ and $M_y$ are constant mean flow Mach number components in $x$
and $y$ directions, respectively. The computational domain embedded
in free space is $-100<x<100$,$-100<y<100$, and $x$- and
$y$-component of the Mach number are
\begin{equation}\label{53b}
    M_x=M_y=0.5 \cos \big(\frac {\pi}{4} \big)
\end{equation}
The initial conditions are Gaussian impulses:
\begin{gather}\label{54}
    p'=e^{-(\ln 2) \big(\frac{x^2+y^2}{9} \big)}, \\
\label{55}
    \rho'=e^{-(\ln 2) \big(\frac{x^2+y^2}{9} \big)}+0.1  e^{-(\ln 2)
  \big(\frac{(x-67)^2+(y-67)^2}{25} \big)}, \\
\label{56}
    u'=0.04 (y-67)  e^{-(\ln 2) \big(\frac{(x-67)^2+(y-67)^2}{25}
    \big)}, \\
\label{57}
    v'=0.04 (x-67)  e^{-(\ln 2) \big(\frac{(x-67)^2+(y-67)^2}{25}
    \big)}.
\end{gather}

The problem is solved using two types of optimized centered finite
difference schemes (second and fourth order). The higher-order schemes
(higher than fourth order) are not referred in this section, although
the anisotropy optimization is acting properly. The domain was discretized
using equally-spaced grid on $x$- and $y$-directions
(table \ref{table2} summarizes the number of points of the grids).
One and two rows of grid
points were considered outside the boundary for the application of
non-reflecting boundary conditions Tam and Webb \cite{Tam93, Tam94}.

\begin{table}[ht]
\caption{Number of grid points used for different cases.}
\begin{center}
\begin{tabular}{|lr|}
  \hline
  Schemes & Number of grid points \\
  \hline
  Second-order  & 400\\
  \hline
  Fourth-order  & 300  \\
  \hline
\end{tabular}
\end{center}
\label{table2}
\end{table}

Filtering techniques with constant coefficients  \cite{Ken} were included to
annihilate the spurious waves: sixth order filters are used.
Low-Dissipation and -Dispersion 4-6
Runge-Kutta Scheme of Hu et. al  \cite{Hu}  is used for time
integration. Results for time equal to $80$ are given in
Figures \ref{Fig3} and  \ref{Fig4}. The numerical results are
compared to analytical results.

\begin{figure}[htp]
  \begin{center}
\includegraphics[width=0.4\textwidth]{fig3a}
\includegraphics[width=0.55\textwidth]{fig3b}  
  \end{center}
  \caption{(a) Density contours and (b) density distribution along the x=y
   direction for the Computational Aeroacoustics Workshop problem using
  optimized second order schemes.}
  \label{Fig3}
\end{figure}

\begin{figure}[htp]
  \begin{center}
 \includegraphics[width=0.4\textwidth]{fig4a}
 \includegraphics[width=0.55\textwidth]{fig4b}
%    \mbox{
%      \subfigure {\label{Fig4a}\includegraphics[scale=0.8]{fig4a.eps}}
%      \subfigure {\label{Fig4b}\includegraphics[scale=0.8]{fig4b.eps}}
%      }
 \end{center}
  \caption{ (a) Density contours and (b) density distribution along the $x=y$
   direction for the Computational Aeroacoustics Workshop problem using
  optimized fourth order schemes.}
  \label{Fig4}
\end{figure}


In addition, some numerical tests on the stationary fluid
($Mx = My = 0$) were conducted: acoustic wave
is propagating from origin, and the entropy wave is neglected. The
domain was extended to $-400<x<400$, $-400<y<400$, such that the
wave can propagate for a longer distance. The front waves of the
acoustic pressure on $x$, and $x=y$ directions are compared
(Figure \ref{Fig5}). The second order optimized centered schemes and
their corresponding classical schemes are tested and compared this
time in order to reveal the anisotropy correction.
Figure \ref{Fig5a} shows that the front waves computed using
classical schemes do not coincide. In Figure \ref{Fig5b} it is seen
that the two front waves are matching by using
optimized schemes.

\begin{figure}[htp]
  \begin{center}
\subfigure {\label{Fig5a}\includegraphics[scale=0.8]{fig5a}}
\subfigure {\label{Fig5b}\includegraphics[scale=0.8]{fig5b}}
 \end{center}
  \caption{Superposed front waves on the $x$, and $x=y$ directions 
  using (a) classical and (b) optimized second order schemes.}
  \label{Fig5}
\end{figure}

\subsection*{Concluding Remarks}
Anisotropy correction of multidimensional finite difference
schemes was carried out for interior stencils. The optimized
schemes incorporate a parameter called isotropy corrector factor
(ICF) which can lower  the isotropy error to a large extent. Based
on Fourier analysis ICF was found and it was shown that, in terms
of isotropy error, the optimized schemes are more effective
compared to classical schemes. It was shown that the optimized
schemes preserve the characteristics of the corresponding
classical onedimensional schemes for all spatial directions. The
optimized schemes are restricted to generalized curvilinear
transformation and to Cartesian grid methods.

\section*{Appendix}

In this appendix, it is shown that the order of accuracy of the
second- and fourth-order optimized centered schemes is the same as
that of the corresponding classical schemes. Only the coefficients
in the Taylor series expansions and their summations are written.
The grid is equally-spaced.

The coefficients of the Taylor series expansions for the
second-order optimized scheme on the $x$ direction are:
\begin{gather}\label{a1}
    \frac {\partial ^0 u}{\partial x^0}: \quad
\frac{1}{2h(1+\beta)}\big[u_{i,j}-u_{i,j}+\frac{\beta}{2}(u_{i,j}-u_{i,j}+u_{i,j}-u_{i,j}
) \big], \\
\label{a2}
    \frac {\partial u}{\partial x}: \quad
\frac{h}{2h(1+\beta)}\big[1+1+\frac{\beta}{2}(1+1+1+1)
\big]=\frac{2h(1+\beta)}{2h(1+\beta)}=1, \\
\label{a3}
    \frac {\partial u}{\partial y}: \quad
\frac{h}{2h(1+\beta)}\big[0+\frac{\beta}{2}(1-1-1+1) \big]=0,\\
\label{a4}
    \frac {\partial^2 u}{\partial x^2}: \quad
\frac{h^2}{2h(1+\beta)}\big[\frac{1}{2!}-\frac{1}{2!}+\frac{\beta}{2}\frac{1}{2!}(1-1+1-1)
\big]=0, \\
\label{a5}
    \frac {\partial^2 u}{\partial x \partial y}: \quad
\frac{h^2}{2h(1+\beta)}\big[0+\frac{\beta}{2}\frac{1}{2!}(2-2-2+2)
\big]=0, \\
\label{a6}
    \frac {\partial^2 u}{\partial y^2}: \quad
\frac{h^2}{2h(1+\beta)}\big[0+\frac{\beta}{2}\frac{1}{2!}(1-1+1-1)
\big]=0.
\end{gather}

The coefficients of the Taylor series expansions for the
fourth-order optimized scheme on $x$ direction are:
\begin{gather}\label{a7}
    \frac {\partial ^0 u}{\partial x^0}: \quad
\frac{1}{12h(1+\beta)}\big[1-8+8-1+\frac{\beta}{2}(1-8+8-1+1-8+8-1)
 \big]=0,\\
\label{a8} \begin{aligned}
    \frac {\partial u}{\partial x}: \quad
&\frac{h}{12h(1+\beta)}\big[-2+8+8-2+\frac{\beta}{2}(-2+8+8-2-2+8+8-2)
\big]\\
&=\frac{12h(1+\beta)}{12h(1+\beta)}=1,
\end{aligned} \\
\label{a9}
    \frac {\partial u}{\partial y}: \quad
\frac{h}{12h(1+\beta)}\big[0+\frac{\beta}{2}(-2+8+8-2+2-8-8+2)
\big]=0, \\
\label{a10}
    \frac {\partial^2 u}{\partial x^2}: \quad
\frac{h^2}{12h(1+\beta)2!}\big[4-8+8-4+\frac{\beta}{2}(4-8+8-4+4-8+8-4)
\big]=0, \\
\label{a11}
    \frac {\partial^2 u}{\partial x \partial y}: \quad
\frac{h^2}{12h(1+\beta)2!}\big[0+\frac{\beta}{2}(8-16+16-8+8-16+16-8)
\big]=0, \\
\label{a12}
    \frac {\partial^2 u}{\partial y^2}: \quad
\frac{h^2}{12h(1+\beta)2!}\big[0+\frac{\beta}{2}(4-8+8-4+4-8+8-4)
\big]=0, \\
\label{a13}
    \frac {\partial^3 u}{\partial x^3}: \quad
\frac{h^3}{12h(1+\beta)3!}\big[-8+8+8-8+\frac{\beta}{2}(-8+8+8-8-8+8+8-8)
\big]=0, \\
\label{a14}
    \frac {\partial^3 u}{\partial x^2 \partial y}: \quad
\frac{h^3}{12h(1+\beta)3!}\big[0+\frac{\beta}{2}(-24+24+24-24+24-24-24+24)
\big]=0, \\
\label{a15}
    \frac {\partial^3 u}{\partial x \partial y^2}: \quad
\frac{h^3}{12h(1+\beta)3!}\big[0+\frac{\beta}{2}(-24+24+24-24-24+24+24-24)
\big]=0, \\
\label{a16}
    \frac {\partial^3 u}{\partial y^3}: \quad
\frac{h^3}{12h(1+\beta)3!}\big[0+\frac{\beta}{2}(-8+8+8-8+8-8-8+8)
\big]=0, \\
\label{a17}
    \frac {\partial^4 u}{\partial x^4}: \quad
\frac{h^4}{12h(1+\beta)4!}\big[16-8+8-16+\frac{\beta}{2}(16-8+8-16+16-8+8-16)
\big]=0, \\
\label{a18}
    \frac {\partial^4 u}{\partial x^3 \partial y}: \quad
\frac{h^4}{12h(1+\beta)4!}\big[0+\frac{\beta}{2}(64-32+32-64-64+32-32+64)
\big]=0, \\
\label{a19}
    \frac {\partial^4 u}{\partial x^2 \partial y^2}: \quad
\frac{h^4}{12h(1+\beta)4!}\big[0+\frac{\beta}{2}(96-48+48-96+96-48+48-96)
\big]=0, \\
\label{a20}
    \frac {\partial^4 u}{\partial x \partial y^3}: \quad
\frac{h^4}{12h(1+\beta)4!}\big[0+\frac{\beta}{2}(64-32+32-64-64+32-32+64)
\big]=0, \\
\label{a21}
    \frac {\partial^4 u}{\partial y^4}: \quad
\frac{h^4}{12h(1+\beta)4!}\big[0+\frac{\beta}{2}(16-8+8-16+16-8+8-16)
\big]=0\,.
\end{gather}

\begin{thebibliography}{10}

\bibitem{Ash}
G.~Ashcroft and X.~Zhang.
\newblock Optimized prefactored compact schemes.
\newblock {\em J. Comput. Phys.}, 190:459--477, 2003.

\bibitem{Har}
J.~C. Hardin, J.~R. Ristorcelli, and C.~K.~W. Tam.
\newblock I{C}{A}{S}{E}/{L}a{R}{C} {Workshop} on {Benchmark} {Problems} in
  {Computational} {Aeroacoustics} ({C}{A}{A}).
\newblock {\em NASA Conference Publication}, 3300, 1995.

\bibitem{Hix}
R.~Hixon.
\newblock Prefactored small-stencil compact schemes.
\newblock {\em J. Comput. Phys.}, 165:522--541, 2000.

\bibitem{Hu}
F.~Q. Hu, M.~Y. Hussaini, and J.~L. Manthey.
\newblock Low-{Dissipation} and �{Dispersion} {Runge}-{Kutta} {Schemes} for
  {Computational} {Acoustics}.
\newblock {\em J. Comput. Phys.}, 124:177--191, 1996.

\bibitem{Ken}
C.~A. Kennedy and M.~H. Carpenter.
\newblock Comparison of {Several} {Numerical} {Methods} for {Simulation} of
  {Compressible} {Shear} {Layers}.
\newblock {\em Appl. Num. Math.}, 14:397--433, 1994.

\bibitem{Kim}
J.~W. Kim and D.~J. Lee.
\newblock Optimized compact finite difference schemes with maximum resolution.
\newblock {\em AIAA Journal}, 34 (5):887--893, 1996.

\bibitem{Lel}
S.~K. Lele.
\newblock Compact finite difference schemes with spectral-like resolution.
\newblock {\em J. Comput. Phys.}, 103:16--42, 1992.

\bibitem{Li}
Y.~Li.
\newblock Wave-{Number} {Extended} {High}-{Order} {Upwind}-{Biased}
  {Finite}-{Difference} {Schemes} for {Convective} {Scalar} {Transport}.
\newblock {\em J. Comput. Phys.}, 133:235--255, 1997.

\bibitem{Lin}
R.~K. Lin and Tony W.~H. Sheu.
\newblock Application of dispersion-relation-preserving theory to develop a
  two-dimensional convection-diffusion scheme.
\newblock {\em J. Comput. Phys.}, 208:493, 2005.

\bibitem{Mah}
K.~Mahesh.
\newblock A family of high order finite difference schemes with good spectral
  resolution.
\newblock {\em J. Comput. Phys.}, 145:332, 1998.

\bibitem{Tam93}
C.~K.~W. Tam and J.~C. Webb.
\newblock Dispersion-{Relation}-{Preserving} {Finite} {Difference} {Schemes}
  for {Computational} {Aeroacoustics}.
\newblock {\em J. Comput. Phys.}, 107:262--281, 1993.

\bibitem{Tam94}
C.~K.~W. Tam and J.~C. Webb.
\newblock Radation {Boundary} {Condition} and {Anisotropy} {Correction} for
  {Finite} {Difference} {Solutions} of the {Helmholtz} {Equation}.
\newblock {\em J. Comput. Phys.}, 113:122--133, 1994.

\bibitem{Tre}
L.~N. Trefethen.
\newblock Group {Velocity} in {Finite} {Difference} {Schemes}.
\newblock {\em SIAM Rev.}, 24:113, 1982.

\bibitem{Vic}
R.~Vichnevetsky and J.~B. Bowles.
\newblock {\em Fourier {Analysis} of {Numerical} {Approximations} of
  {Hyperbolic} {Equations}}.
\newblock SIAM Studies in Applied Mathematics, Philadelphia, 1982.

\bibitem{Whi}
G.~B. Whitham.
\newblock {\em Linear and {Nonlinear} {Waves}}.
\newblock A Wiley-Interscience Publication, New York, 1974.

\bibitem{Zin}
D.~W. Zingg and H.~Lomax.
\newblock Finite {Difference} {Schemes} on {Regular} {Triangular} {Grids}.
\newblock {\em J. Comput. Phys.}, 108:306--313, 1993.

\end{thebibliography}

\end{document}