% beamer version of chapter 5 "The Pseudospectral method Method" in "Computational Seismology, A Practical Introduction" (H. Igel)

%\documentclass{beamer}
%
%\mode<presentation>
%{
  %\usetheme{CambridgeUS}
%
  %\setbeamercovered{transparent}
  %\setbeamertemplate{navigation symbols}{}
%}


\documentclass[aspectratio=169]{beamer}

\mode<presentation>
{
  %\usetheme{metropolis}
	\usetheme{metropolis}
  %\usecolortheme{beaver}
	\usecolortheme{crane}

  \setbeamercovered{transparent}
  \setbeamertemplate{navigation symbols}{}
}


\usepackage[english]{babel}
\usepackage[latin1]{inputenc}

\usepackage{times}
\usepackage[T1]{fontenc}
% Or whatever. Note that the encoding and the font should match. If T1
% does not look nice, try deleting the line with the fontenc.

% Math packages
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{mathtools}
\usepackage{graphicx}
%using a different letter style for capital letters
\usepackage{mathrsfs}

\usepackage{mdframed}
\usepackage{listings}
\usepackage[percent]{overpic}
\usepackage{capt-of}
\usepackage{tikz}
\usepackage{grffile}

% Own definitions ... 
\newcommand{\eqdef}{\mathrel{\mathop:}=} % for symbol :=
\newcommand*\Eval[3]{\left.#1\right\rvert_{#2}^{#3}} % for integration bound |
\newcommand{\figpath}{c:/mydata/Heiner/Text/Books/OUP/TheBook/Part2/Figs_ps}
\newcommand{\fs}{\footnotesize}

\title[Computational Seismology] % (optional, use only with long paper titles)
{The Pseudospectral Method}

%\subtitle{Introduction} 

\author[Heiner Igel]{Heiner Igel}


\institute[] 
  {Department of Earth and Environmental Sciences \\ Ludwig-Maximilians-University Munich} 



\date[] % (optional)
{}

\subject{Talks}


%\beamerdefaultoverlayspecification{<+->}

\begin{document}

% ----------------------------------------------------------

\begin{frame}
  \titlepage
\end{frame}

% ----------------------------------------------------------

%\begin{frame}{Outline}
%\fs
  %\tableofcontents
  %% You might wish to add the option [pausesections]
%\end{frame}

% ---------------------------------------------------------

\section{Introduction}

% ---------------------------------------------------------
\subsection{Motivation}

	\begin{frame}{Motivation}
		\begin{minipage}{0.4\textwidth}
			\begin{figure}
				\includegraphics[width=0.9\textwidth]{\figpath/fig_ssection.png}
			\end{figure}
		\end{minipage}
	\mbox{ }
    \hfill	
		\begin{minipage}{0.55\textwidth}
			\begin{block}{}
				\begin{enumerate}
					\item Orthogonal basis functions, special case of FD
 					\item Spectral accuracy of space derivatives
					\item High memory efficiency
					\item Explicit method
					\item No requirement of grid staggering
					\item Problems with  strongly heterogeneous media
				\end{enumerate}
			\end{block}
		\end{minipage}
	\end{frame}

% ---------------------------------------------------------

\subsection{History}

	\begin{frame}{History}
	\scriptsize
		\begin{itemize}[<+->]
			\item Coining as transform methods as their implementation was based on the Fourier transform (Gazdag, 1981; Kossloff and Bayssal, 1982)
			\item Initial applications to the acoustic wave equation were extended to the elastic case (Kossloff et al., 1984), and to 3D (Reshef et al., 1988)
			\item Developing efficient time integration schemes (Tal-Ezer et al., 1987) that allowed large times steps to be used in the extrapolation procedure
			\item Replacing harmonic functions as bases for the function interpolation by Chebyshev poly-nomials (Kosloff et al., 1990)
			\item To improve the accurate modelling of curved internal interfaces and surface topography grid stretching as coordinate transforms was introduced and applied (Tessmer et al., 1992; Komatitsch et al., 1996)
			\item By mixing finite-difference operators and pseudospectral operators in the different spatial directions, the method was used for interesting seismological problems (Furumura et al., 1998b; Furumura and Kennett, 2005)
		\end{itemize}
	\end{frame}

% ---------------------------------------------------------

\subsection{The Pseudospectral Method in a Nutshell}

	\begin{frame}{The Pseudospectral Method in a Nutshell}
		The Pseudospectral method is:
		\begin{itemize}
			\item a grid point method
			\item a series expansion method (Fourier or Chebyshev)
		\end{itemize}
		\vspace*{1cm}
		Looking at the acoustic wave equation using finite-difference method leaves us with\\
		\begin{equation*}
			\frac{p(x,t+\mathrm{d}t) - 2 p(x,t) + p(x,t-\mathrm{d}t)}{\mathrm{d}t^2} \ = \ c(x)^2 \partial_x^2 p(x,t) + s(x,t)
		\end{equation*}
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{The Pseudospectral Method in a Nutshell}
		The remaining task is to calculate the space derivative on the r.h.s.\\
		\begin{equation*}
			\partial_x^{(n)} p(x,t) \ = \ \mathscr{F}^{-1}[(-ik)^n P(k,t)] 
		\end{equation*} 
	\vfill
		where $i$ is the imaginary unit, $\mathscr{F}^{-1}$ is the inverse Fourier transform, and $P(k,t)$ is the spatial Fourier transform of the pressure field $p(x,t)$, $k$ being the wavenumber.
	\vfill
		Using discrete Fourier transform of functions defined on a regular grid, we obtain exact derivatives up to the Nyquist wavenumber $k_{N}=\pi/dx$.
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{The Pseudospectral Method in a Nutshell}
		\begin{minipage}{0.4\textwidth}
			\includegraphics[width=0.6\textwidth]{\figpath/fig_ps_nutshell.png}
		\end{minipage}
	\mbox{ }
    \hfill	
    \begin{minipage}{0.55\textwidth}
		\vspace{-1truecm}
    	\begin{block}{Principle of the pseudospectral method based on the Fourier series}
			\begin{itemize}
				\item Use of sine and cosine functions for the expansions implies periodicity
				\item Using Chebyshev polynomials similar accuracy of common boundary conditions (free surface, absorbing) can be achieved
			\end{itemize}
    	\end{block}
    	\vspace{0.5cm}
	\end{minipage}    		
	\end{frame}
	
% ---------------------------------------------------------

\section{The Pseudospectral Method: Ingredients}
	
% ---------------------------------------------------------

	\begin{frame}{Orthogonal Functions, Interpolation, Derivative}
		\begin{block}{In many situations we either...}
			\begin{itemize}
				\item seek to approximate a known analytic function by an approximation
				\item know a function only at a discrete set of points and we would like to interpolate in between those points
			\end{itemize}
		\end{block}
		\vspace*{0.5cm}
		Let us start with the first problem such that our known function is approximated by a finite sum over some N basis functions $\Phi_i$ 
		\begin{equation*}
			f(x) \ \approx \ g_N (x) \ = \ \sum_{i=1}^N a_i \Phi_i (x)
		\end{equation*}
		and assume that the basis functions form an orthogonal set
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{Orthogonal Functions, Interpolation, Derivative}
		Why would one want to replace a known function by something else?
		\vfill
		\begin{minipage}{0.4\textwidth}
			\includegraphics[width=0.9\textwidth]{\figpath/fig_nondiff.png}
		\end{minipage}
	\mbox{ }
    \hfill	
    	\begin{minipage}{0.55\textwidth}
    		\begin{block}{}
    			Dynamic phenomena are mostly expressed by PDEs
    			\begin{itemize}
    				\item[] Either nature is not smooth and differentiable
    				\item[] mathemical functions are non-differentiable
    			\end{itemize}
    		\end{block}
		\end{minipage}    		
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{Orthogonal Functions, Interpolation, Derivative}
		With the right choice of differentiable basis functions $\Phi_i$ the calculation becomes
		\begin{equation*}
			\partial_x f(x) \ \approx \ \partial_x g_N (x) \ = \ \sum_{i=1}^N a_i \partial_x \Phi_i (x) 
		\end{equation*}
		Consider the set of (trigonometric) basis functions
		\begin{equation*}
			\begin{split}
				\cos (nx) &\qquad n \ = \ 0, 1, \dotsc , \infty \\
				\sin (nx) &\qquad n \ = \ 0, 1, \dotsc , \infty 
			\end{split}
		\end{equation*}
		with
		\begin{equation*}
			\begin{split}
				1, &\cos (x), \cos (2x), \cos (3x), \dotsc  \\
				0, &\sin (x), \sin (2x), \sin (3x), \dotsc
			\end{split}
		\end{equation*}
		in the interval $[-\pi, \pi]$
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{Orthogonal Functions, Interpolation, Derivative}
		\begin{minipage}{0.25\textwidth}
			\begin{figure}
				\includegraphics[width=0.85\textwidth]{\figpath/fig_ortho.png}
			\end{figure}
		\end{minipage}
	\mbox{ }
    \hfill			
		\begin{minipage}{0.7\textwidth}
		Checking whether these functions are orthogonal by evaluating integrals with all possible combinations
		\begin{equation*}
			\begin{split}
				&\int \limits_{-\pi}^{\pi} \cos(jx) \cos(kx) \mathrm{d}x \ = \
   				 \begin{cases} 0    &\mbox{for } j \ \neq \ k \\
    	        		  2\pi &\mbox{for } j \ = \ k \ = \ 0 \\
    	       		  	 \pi  &\mbox{for } j \ = \ k \ > \ 0
				\end{cases} \\
				&\int \limits_{-\pi}^{\pi} \sin(jx) \sin(kx) \mathrm{d}x \ = \
    			\begin{cases} 0    &\mbox{for } j \ \neq \ k \ ; \ j, k \ > \ 0 \\
    	      	  	  \pi  &\mbox{for } j \ = \ k \ > \ 0
				\end{cases} \\
				&\int \limits_{-\pi}^{\pi} \cos(jx) \sin(kx) \mathrm{d}x \ = \
   				 0    \mbox{  for } j \ \ \geqslant \ 0 \ , \ k \ >  \ 0
			\end{split}
		\end{equation*}
		\end{minipage}
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{Orthogonal Functions, Interpolation, Derivative}
		The approximate function $g_N (x)$ can be stated as
		\begin{equation*}
			f(x) \ \approx \ g_N (x) \ = \ \sum_{k=0}^N a_k \cos(kx) \ + \ b_k \sin(kx) 
		\end{equation*}  
		By minimizing the difference between approximation $g_N (x)$ and the original function $f(x)$, the so-called $l_2$-norm, the coefficients $a_k, b_k$ can be found
		\begin{equation*}
			\| f(x) - g_N(x) \|_{l_2} \ = \ \left[ \int \limits_a^b 
   			 \lbrace f(x) - g_N(x) \rbrace ^2 \, \mathrm{d}x \right]
  			  ^{ \tfrac{1}{2} } \ = \ \text{Min} 
		\end{equation*}
		$\Longrightarrow$ independent of the choice of basis functions
	\end{frame}
		
% ---------------------------------------------------------

\subsection{Fourier Series and Transforms}

	\begin{frame}{Fourier Series and Transforms}
		The most important concept of this section will consist of the properties of Fourier series on regular grids.\\
		\vspace{0.3cm}
		The approximate function $g_{N} (x)$ has the following form
		\begin{equation*}
			g_N (x) \ = \ \frac{1}{2} a_0 + \sum _{k = 1} ^n a_k \cos (kx) + b_k \sin (kx) 
		\end{equation*}
		and leads to the coefficients
		\begin{equation*}
			\begin{split}
				a_k \ &= \ \frac{1}{\pi} \int \limits _{-\pi}^\pi f(x) \cos (kx) \mathrm{d}x 
			    \qquad k \ = \ 0, 1, \dotsc , n \\
				b_k \ &= \ \frac{1}{\pi} \int \limits _{-\pi}^\pi f(x) \sin (kx) \mathrm{d}x 
 			   \qquad k \ = \ 0, 1, \dotsc , n \ \ .
			\end{split}
		\end{equation*}
	\end{frame}
			
% ---------------------------------------------------------

	\begin{frame}{Fourier Series and Transforms}
		Using Euler's formulae, yields to
		\begin{equation*}
			g_N (x) \ = \sum _{k = -n} ^{k=n} c_k e^{ikx}
		\end{equation*}
		with complex coefficients $c_k$ given by		
		\begin{equation*}
			\begin{split}
				c_k \   &= \frac{1}{2} \ (a_k - i b_k) \\
				c_{-k} \  &= \frac{1}{2} \ (a_k + i b_k) \ \ \ k \ > \ 0 \\
				c_0 \  &=  \frac{1}{2} a_0 \ .
			\end{split}
		\end{equation*}
	\end{frame}
				
% ---------------------------------------------------------

	\begin{frame}{Fourier Series and Transforms}
		\begin{minipage}{0.3\textwidth}
			\begin{figure}
				\includegraphics[width=0.75\textwidth]{\figpath/fig_continuous.png}
			\end{figure}
		\end{minipage}
	\mbox{ }
    \hfill	
    	\begin{minipage}{0.65\textwidth}
    		Finding the interpolating trigonometric polynomial for the periodic function
    		\begin{equation*}
				f(x + 2\pi x ) \ = \ f(x) \ = \ x^2 \qquad x \in \left[ 0, 2\pi \right]
			\end{equation*}
			The approximation $g_N (x)$ can be obtained with
			\begin{equation*}
				g_N (x) \ = \ \frac{4\pi^2}{3} + \sum _{k=1} ^N \lbrace \frac{4}{k^2} \cos (kx)
  				 - \frac{4\pi}{k} \sin (kx) \rbrace
			\end{equation*}
		\end{minipage}    			
	\end{frame}
					
% ---------------------------------------------------------

	\begin{frame}{Fourier Series and Transforms}
		We assume that we know our function f(x) at a discrete set of points $x_i$ given by
		\begin{equation*}
			x_i \ = \ \frac{2\pi}{N} i \ \ \ i=0, \dots, N \ .
		\end{equation*}
		Using the "trapezoidal rule" for the integration of a definite integral we obtain for the Fourier coefficients
		\begin{align*}
			a_k^* \ &= \ \frac{2}{N} \sum _{j=1}^N f(x_j) \cos (kx_j) 
  			  \qquad &k \ = \ 0, 1, \dotsc , n \\
			b_k^* \ &= \ \frac{2}{N} \sum _{j=1}^N f(x_j) \sin (kx_j)  
   			 \qquad &k \ = \ 0, 1, \dotsc , n
		\end{align*}
	\end{frame}
						
% ---------------------------------------------------------

	\begin{frame}{Fourier Series and Transforms}
		\begin{minipage}{0.3\textwidth}
			\begin{figure}
				\includegraphics[width=0.75\textwidth]{\figpath/fig_discrete.png}
			\end{figure}
		\end{minipage}
	\mbox{ }
    \hfill			
		\begin{minipage}{0.65\textwidth}
		We thus obtain the specific Fourier polynomial with $N = 2n$
		\begin{equation*}
		\begin{split}
			g_n^* \ \coloneqq \ \ &  \frac{1}{2} a_0^* + \sum _{k=1} ^{n-1} \lbrace a_k^* \cos (kx)
   			 - b_k^* \sin (kx) \rbrace\\
   			  \ & + \frac{1}{2} a_n^* cos(nx)
   		\end{split}
		\end{equation*}
		with the tremendously important property
		\begin{equation*}
			g_n^* (x_i) \ = \ f(x_i)    \ .
		\end{equation*}
		\end{minipage}
	\end{frame}
							
% ---------------------------------------------------------

	\begin{frame}{Cardinal functions}
		\begin{minipage}{0.4\textwidth}
			\begin{figure}
				\includegraphics[width=0.95\textwidth]{\figpath/fig_cardinal.png}
			\end{figure}
		\end{minipage}
	\mbox{ }
    \hfill	
    	\begin{minipage}{0.55\textwidth}
    		\begin{itemize}
    			\item[] Discrete interpolation and derivative operations can also be formulated in terms of convolutions
    			\item[] It is unity at grid point $x_i$ and zero at all other points on the discrete grid
    			\item[] It has the form of a $sinc$-function
    		\end{itemize}
		\end{minipage}    			
	\end{frame}
						
% ---------------------------------------------------------

	\begin{frame}{Fourier Series and Transforms}
		\begin{block}{Forward Transform}
			\begin{equation*}
				F(k) \ = \mathscr{F}[f(x)] \ = \ \frac{1}{\sqrt{2\pi}}  \ \int^\infty_{-\infty} f(x) e^{ikx} dx 
			\end{equation*}
			\end{block}
			\vspace{1cm}
			\begin{block}{Inverse Transform}
			\begin{equation*}
				f(x) \ = \ \mathscr{F}^{-1}[F(k)] \ = \frac{1}{\sqrt{2\pi}} \ \int^\infty_{-\infty} F(k) e^{-ikx} dk 
			\end{equation*}
		\end{block}
	\end{frame}
						
% ---------------------------------------------------------

	\begin{frame}{Fourier Series and Transforms}
		Taking the formulation of the inverse transform to obtain the derivative of function f(x)
		\begin{equation*}
			\begin{split}
				\frac{d}{dx}f(x) \ &= \ \frac{d}{dx} \frac{1}{\sqrt{2\pi}} \ \int^\infty_{-\infty} F(k) e^{-ikx} dk \\
				&= \  \frac{1}{\sqrt{2\pi}} \ \int^\infty_{-\infty} - ik \ F(k) e^{-ikx} dk \\
				&= \  \frac{1}{\sqrt{2\pi}} \  \int^\infty_{-\infty} D(k) \ F(k) e^{-ikx} dk
			\end{split}
		\end{equation*}
		with $D(k) = -ik$\\
		We can extend this formulation to the calculation of the $n-th$ derivative of $f(x)$  to achieve
		\begin{equation*}
			F^{(n)}(k) \ = \ D(k)^n \ F(k) \ = \ (-ik)^n \ F(k)
		\end{equation*}
	\end{frame}
						
% ---------------------------------------------------------

	\begin{frame}{Fourier Series and Transforms}
		Thus using the condense Fourier transform operator $\mathscr{F}$ we can obtain an exact $n-th$ derivative using
		\begin{equation*}
			\begin{split}
				f^{(n)}(x) \ &= \ \mathscr{F}^{-1}[ (-ik)^n \ F(k) ] \\
             &= \ \mathscr{F}^{-1}[ (-ik)^n \ \mathscr{F}[f(x)] ] \ .
			\end{split}
		\end{equation*}
		Adopting the complex notation of the forward transform we gain
		\begin{equation*}
			F_k \ = \ \sum_{j=0}^{N-1} f_j \ e^{i\ 2\pi j k/N} \ \ k=0,\dots,N
		\end{equation*}
		and the inverse transform 
		\begin{equation*}
			f_j \ = \ \frac{1}{N} \ \sum_{k=0}^{N-1} F_k \ e^{-i\ 2\pi j k/N} \ \ j=0,\dots,N
		\end{equation*}
	\end{frame}
						
% ---------------------------------------------------------

	\begin{frame}{Fourier Series and Transforms}
		We are able to gain exact $n - th$ derivatives on our regular grid by performing the following operations on vector $f_j$ defined at grid points $x_j$
		\vfill
		\begin{equation*}
			\partial^{(n)}_x f_j  \  = \mathscr{F}^{-1} [ (-ik)^n \ F_k]
		\end{equation*}
		\vfill
		where
		\begin{equation*}
			F_k  = \mathscr{F}[f_j]
		\end{equation*} 
	\end{frame}
						
% ---------------------------------------------------------

	\begin{frame}{Example}
		\begin{minipage}{0.55\textwidth}
			We initialize a $2 \pi$-periodic Gauss-function in the interval $x\in[0,2\pi]$ as
			\begin{equation*}
				f(x)\ = \ e^{-1/\sigma^2 \ (x-x_0)^2}  
			\end{equation*}
			with $x_0=\pi$ and the  derivative
			\begin{equation*}
				f'(x)\ = \ - 2 \frac{(x-x_0)}{\sigma^2} \ e^{-1/\sigma^2 \ (x-x_0)^2}  
			\end{equation*}
			 Grid spacing of \small{$dx=\frac{2\pi}{N}$ with $N=127$ and $x_j=j \frac{2\pi}{N}, \ j=0,\dots,N$}. 
		\end{minipage}
	\mbox{ }
    \hfill	
    	\begin{minipage}{0.39\textwidth}
    	\tiny
    		\begin{mdframed}[style=MyFrame]
{\sffamily
{\tiny
\% Main program \\
\% \\
\% Initialize function vector f \\
\%(...) \\
\% Calculate derivative of vector f in interval [a,b]\\
\indent df =  fder(f,a,b) \\
\%  \\
\% (...) \\
\% Subroutines/Functions \\
\indent function df = fder(f,a,b) \\
\indent \% Fourier Derivative of vector f \\
\indent \% (...) \\
\indent \% length of vector f \\
\indent n = max(size(f)); \\
\indent \% initialization of k vector (wavenumber axis) \\
\indent k = 2*pi/(b-a)*[0:n/2-1 \\
\indent 0 -n/2+1:-1]; \\
\indent \% Fourier derivative \\
\indent df = ifft(-i*k.*fft(f)); \\
(...)\\
}}
\end{mdframed}
		\end{minipage}    			
	\end{frame}
						
% ---------------------------------------------------------

	\begin{frame}{Result}
			\begin{figure}
				\includegraphics[width=0.75\textwidth]{\figpath/fig_fourierderivative_gauss3.png}
			\end{figure}
	\end{frame}

% ---------------------------------------------------------

\section{The Fourier Pseudospectral Method}
\subsection{Acoustic 1D}

	\begin{frame}
		\begin{block}{}
			\begin{center}
				\huge 	The Fourier Pseudospectral Method	
			\end{center}
		\end{block}
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{Acoustic 1D}
		\begin{block}{Constant-density acoustic wave equation in 1D}
			\begin{equation*}
				\ddot{p} \ = \ c^2 \partial_x^2 p + s  \ 
			\end{equation*}
		\end{block}
		\vfill
		The time-dependent part is solved using a standard 3-point finite-difference operator leading to  
		\begin{equation*}
			\frac{p_j^{n+1} - 2 p_j^n + p_j^{n-1}}{\mathrm{d}t^2} \ = \ c_j^2 \partial_x^2 p_j^n
    		 + s^n_j  \ 
		\end{equation*}
		where upper indices represent time and lower indices space.
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{Acoustic 1D}
		\begin{minipage}{0.55\textwidth}
			Calculating the $2^{nd}$ derivatives using the Fourier transform
			\begin{equation*}
				\begin{split}
 					\partial_x^2 p_j^n \ \ & =  \ \mathscr{F}^{-1} [ (-ik)^2 \ P^n_\nu] \ \\
 					\ & = \ \mathscr{F}^{-1} [ -k^2 \ P^n_\nu]
 				\end{split}
			\end{equation*}	
			where $P^n_\nu$ is the discrete complex wavenumber spectrum at time $n$ \\
			leading to an exact derivative with only numerical rounding errors.	
		\end{minipage}
	\mbox{ }
    \hfill		
    \begin{minipage}{0.4\textwidth}
    \tiny		
    \hfill
			\begin{mdframed}[style=MyFrame]
{\sffamily
{\tiny
\% Main program \\
\%(...) \\
\% Time exploration\\
for i=1:nt, \\
\% (...) \\
\% 2nd space derivate  \\
\indent d2p=s2der1d(p,dx);   \\
\% Extrapolation \\
\indent pnew=2*p-pold+c.*c.*d2p*dt*dt;  \\
\% Add sources \\
\indent pnew=pnew+sg*src(i)*dt*dt;  \\
\% Remap pressure field \\
\indent pold=p; \\
\indent p=pnew; \\
\% (...) \\
end \\
\% (...) \\
\% Subroutines \\
function df=s2der1d(f,dx) \\
\% (...) \\
\% 2nd Fourier derivative \\
ff=fft(f); ff=k.*k.*ff; df=real(ifft(ff)); \\
}}
\end{mdframed}
		\end{minipage}
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{Example}
		\begin{minipage}{0.65\textwidth}
			In FD method possible to initiate a point-like source at one grid point \\
			\vspace{0.3cm}
			
			In PS method not possible because Fourier transform of a spike-like function creates oscillations \\
			\vspace{0.3cm}
			
			$\Longrightarrow$ Defining a space-dependent part of the source using a Gaussian function $e^{-1/\sigma^2(x-x_0)^2}$ with $\sigma = 2dx$, $dx$ being the grid interval and $x_0$ the source location
		\end{minipage}
	\mbox{ }
    \hfill			
		\begin{minipage}{0.3\textwidth}
			\begin{table}
				\begin{tabular}[t]{ll}
					Parameter & Value \\
					\hline
					$x_{max}$ & 1250 m \\
					nx & 2048 \\
					c & 343 m/s \\
					dt & 0.00036 s \\
					dx & 0.62 m \\
					$f_0$ & 60 Hz \\ 
					$\epsilon$ & 0.2 \\
				\end{tabular}
			\end{table}
		\end{minipage}
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{Result}
		\begin{figure}
			\includegraphics[width=0.6\textwidth]{\figpath/fig_ac_fourier.png}
		\end{figure}
	\end{frame}
	
% ---------------------------------------------------------

\subsection{Stability, Convergence, Dispersion}

	\begin{frame}{Stability, Convergence, Dispersion}
		To understand the behaviour of numerical approximations using discrete plane waves of the form
		\begin{equation*}
			\begin{split}
				p_j^n \ & = \ e^{i(k j \mathrm{d} x - \omega n \mathrm{d}t)} \\
				\partial_x^2 p_j^n \ & = \ -k^2 e^{i(k j \mathrm{d} x - \omega n \mathrm{d}t)} \\
			\end{split}
		\end{equation*}
		The time-dependent part can be expressed as
		\begin{equation*}
			\partial_t^2 p_j^n \ = \ -\frac{4}{dt^2}\ \sin^2(\frac{\omega dt}{2}) \  e^{i(k j \mathrm{d} x - \omega n \mathrm{d}t)} 
		\end{equation*}
		where we made use of Euler's formula and that $2\sin^2 x=1/2(1-\cos2x)$
	\end{frame}
		
% ---------------------------------------------------------

	\begin{frame}{Stability, Convergence, Dispersion}
		\begin{block}{Phase velocity c(k)}
			\begin{equation*}
				c(k) \ = \  \frac{\omega}{k} \ = \frac{2}{k dt}\ \sin^{-1}(\frac{k c dt}{2})  \ .
			\end{equation*}		
		\end{block}
\vfill
		\begin{itemize}
			\item When $dt$ becomes small  $\sin^{-1}(k c dt/2)\approx k c dt/2$
			\item$dx$ does not appear in this equation
			\item The inverse sine must be smaller than one the stability limit requires $k_{max}(c dt/2) \leq 1$. As $k_{max}=\pi/dx$ the stability criterion for the 1D case is $\epsilon=c dt/dx = 2/\pi \approx 0.64$
		\end{itemize}
	\end{frame}
			
% ---------------------------------------------------------

	\begin{frame}{Stability, Convergence, Dispersion}
		\begin{figure}
			\includegraphics[width=0.5\textwidth]{\figpath/fig_fourier_disp.png}
		\end{figure}
	\end{frame}
				
% ---------------------------------------------------------

\subsection{Acoustic 2D}

	\begin{frame}{Acoustic 2D}
		\begin{block}{Acoustic wave equation in 2D}
			\begin{equation*}
				\ddot{p} \ = \ c^2 (\partial_x^2 p + \partial_z^2 p) \ + s
			\end{equation*}
		\end{block}
		\vfill
		The time-dependent part is replaced by a standard 3-point finite-difference approximation
		\begin{equation*}
 			\frac{p_{j,k}^{n+1} - 2 p_{j,k}^n + p_{j,k}^{n-1}}{\mathrm{d}t^2} \ = \ c_{j,k}^2 ( \partial_x^2 p + \partial_z^2 p)_{j,k} \ + s_{j,k}^n
		\end{equation*}
		Using Fourier approach for approximating $2^{nd}$ partial derivatives
		\begin{equation*}
 			\partial_x^2 p + \partial_z^2 p \ = \mathscr{F}^{-1} [ -k_x^2 \ \mathscr{F}[p] ] \ + \  \mathscr{F}^{-1} [ -k_z^2 \ \mathscr{F}[p] ]
		\end{equation*}
	\end{frame}
				
% ---------------------------------------------------------

	\begin{frame}{Acoustic 2D}
		\begin{minipage}{0.35\textwidth}
			\begin{table}
				\begin{tabular}[t]{ll}
					Parameter & Value \\
					\hline
					$x_{max}$ & 200 m \\
					nx & 256 \\
					c & 343 m/s \\
					dt & 0.00046 s \\
					dx & 0.78 m \\
					$f_0$ & 200 Hz \\ 
					$\epsilon$ & 0.2 \\
				\end{tabular}
			\end{table}
		\end{minipage}
	\mbox{ }
    \hfill				
		\begin{minipage}{0.45\textwidth}
		\tiny
			\begin{mdframed}[style=MyFrame]
{\sffamily
{\tiny
\% (...) \\
\% 2nd space derivates  \\
\indent for  j=1:nz, \\
\indent d2xp(:,j)=s2der1d(p(:,j)',dx); \\
\indent end \\
\indent for  i=1:nx, \\
\indent d2zp(i,:)=s2der1d(p(i,:),dx); \\ 
\indent end \\
\% Extrapolation \\
\indent pnew=2*p-pold+c.*c.*(d2xp+d2zp)*dt*dt2;  \\
\% (...) \\
}}
\end{mdframed}
		\end{minipage}
	\end{frame}
					
% ---------------------------------------------------------

	\begin{frame}{Acoustic 2D}
		\begin{figure}
			\includegraphics[width=0.9\textwidth]{\figpath/fig_ps_2D.png}
		\end{figure}
	\end{frame}
						
% ---------------------------------------------------------

\subsection{Numerical anisotropy}

	\begin{frame}{Numerical anisotropy}
		Investigating the dispersion behaviour by finding solutions to monochromatic plane waves propagating in the direction ${ \bf k} =(k_x,k_z)$   
		\begin{equation*}
			p^n_{j,k} \ = \  e^{i (k_x j \mathrm{d}x + k_z k \mathrm{d}x  - \omega n \mathrm{d}t )} 
		\end{equation*}
		With Fourier method the derivatives can be calculated by
		\begin{equation*}
			\begin{split}
				\partial_x p^n_{j,k} \ &= \ -k_x^2 \ e^{i (k_x j \mathrm{d}x + k_z k \mathrm{d}x  - \omega n \mathrm{d}t )} \\
				\partial_z p^n_{j,k} \ &= \ -k_z^2 \  e^{i (k_x j \mathrm{d}x + k_z k \mathrm{d}x  - \omega n \mathrm{d}t )} 
			\end{split}
		\end{equation*}
		Combining this with the 3-point-operator for the time derivative
		\begin{equation*}
			\partial_t^2 p_{j,k}^n \ = \ -\frac{4}{dt^2}\ \sin^2(\frac{\omega dt}{2}) \   e^{i (k_x j \mathrm{d}x + k_z k \mathrm{d}x
			  - \omega n \mathrm{d}t )}
		\end{equation*}
	\end{frame}
							
% ---------------------------------------------------------

	\begin{frame}{Numerical anisotropy}
		We obtain the numerical dispersion relation in 2D for arbitrary wave number vectors (i.e., propagation directions) ${\bf k}$ as
		\begin{equation*}
			c({\bf k}) \ = \ \frac{\omega}{|{\bf k}|} \ = \ \frac{2}{|{\bf k}| dt}\ \sin^{-1}(\frac{c dt \sqrt{k_x^2+k_z^2}}{2})  \ .
		\end{equation*}
		\begin{figure}
			\includegraphics[width=0.3\textwidth]{\figpath/fig_fourier_anis.png}
		\end{figure}
	\end{frame}
								
% ---------------------------------------------------------

\subsection{Elastic 1D}

	\begin{frame}{Elastic 1D}
		\begin{block}{1D Elastic wave equation}
			\begin{equation*}
				\rho(x) \ddot{u}(x,t) \ = \  \partial_x \left[ \mu(x) \partial_x u(x,t) \right] \ + \  f(x,t)
			\end{equation*}
		\end{block}
		\vfill
		$u \; \; \; \; \; \; \; $ displacement field\\
		$\mu  \; \; \; \; \; \; \; $ space-dependent shear modulus\\
		\vfill
		The finite-difference approximation of the extrapolation part leads to
		\begin{equation*}
			\rho_i \frac{ u_i^{j+1} - 2u_i^j + u_i^{j-1} }{\mathrm{d}t^2} \ = \ (\partial_x \left[ \mu(x) \partial_x u(x,t) \right])_i^j  + f_i^j
		\end{equation*}
		with space derivatives to be calculated using the Fourier method.\\
	\end{frame}
								
% ---------------------------------------------------------

	\begin{frame}{Elastic 1D}			
		The sequence of operations required to obtain the r.h.s. reads\\
		\begin{equation*}
			\begin{split}
				u_i^j  \ & \rightarrow   \mathscr{F}[u_i^j]  \rightarrow \ U^j_\nu \rightarrow -ik U^j_\nu \rightarrow \mathscr{F}^{-1}[-ik U^j_\nu]  \rightarrow \partial_x u_i^j  \\
 \partial_x u_i^j   &\rightarrow  \mathscr{F}[\mu_i  \partial_x u_i^j] \rightarrow \tilde{U}^j_\nu  \rightarrow \mathscr{F}^{-1}[-ik \tilde{U}^j_\nu] \rightarrow \partial_x \left[ \mu(x) \partial_x u(x,t) \right] \\
			\end{split}
		\end{equation*}\\
		
		where capital letters denote fields in the spectral domain, lower indices with Greek letters indicate discrete frquencies, and $\tilde{U}^j_\nu=\mu_i  \partial_x u_i^j$ was introduced as an intermediate result to facilitate notation.
	\end{frame}
								
% ---------------------------------------------------------

	\begin{frame}{Elastic 1D}	
		\begin{minipage}{0.5\textwidth}
			Finding a setup for a classic staggered-grid finite-difference solution to the elastic 1D problem, leads to an energy misfit to the analytical solution $u_a$ of 1\%. The energy misfit is simply calculated by $(u_{FD}-u_{a})^2/u_a^2$
		\end{minipage}
	\mbox{ }
    \hfill			
		\begin{minipage}{0.4\textwidth}
			\begin{table}
				\begin{tabular}[t]{lll}
					 & FD & PS \\
					\hline
					nx & 3000 & 1000 \\
					nt & 2699 & 3211 \\
					c & 3000 m/s & 3000 m/s\\
					dx & 0.33 m & 1.0 m\\
					dt & 5.5e-5 s & 4.7e-5 s\\
					$f_0$ & 260 Hz & 260 Hz \\ 
					$\epsilon$ & 0.5  & 0.14 \\
					$n/\lambda$ & 34  &11  \\
				\end{tabular}
			\end{table}
		\end{minipage}
	\end{frame}
									
% ---------------------------------------------------------

	\begin{frame}{Elastic 1D}	
			\begin{figure}
				\includegraphics[width=0.7\textwidth]{\figpath/fig_el1d_fdps3.png}
			\end{figure}	
			\scriptsize Comparing memory requirements and computation speed between the Fourier method (\textbf{right}) and a 4th-order finite-difference scheme (\textbf{left}). In both cases the relative error compared to the analytical solution (misfit energy calculated by $\frac{u_{FD} - u_{a}}{u_{a}^2}$) is approximately 1\%. The big difference is the number of grid points along the x dimension. The ratio is 3:1 (FD:Fourier)	
	\end{frame}
					
% ---------------------------------------------------------

	\begin{frame}{Summary}	
		\small
		\begin{itemize}
		\item Pseudospectral methods are based on discrete function approximations that allow exact interpolation at so-called collocation points. The most prominent examples are the Fourier method based on trigonometric basis functions and the Chebyshev method based on Chebyshev polynomials.
\item The Fourier method can be interpreted as an application of discrete Fourier series on a regular-spaced grid. The space derivatives can be obtained exactly (except for rounding errors). Derivatives can be efficiently calculated with the discrete Fourier transform requiring $n \log n$ operations.
\item The Fourier method implicitly assumes periodic behavior. Boundary conditions like the free surface or absorbing behaviour are difficult to implement.
		\end{itemize}
	\end{frame}

% ---------------------------------------------------------
						
\end{document}