\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