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


\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{Infinite Order Finite Differences}

	%\begin{frame}
		%\begin{block}{}
			%\begin{center}
				%\huge{Infinite Order\\
				 %Finite Differences}
			%\end{center}
		%\end{block}
	%\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{Infinite Order Finite Differences}
		\begin{minipage}{0.65\textwidth}
		\small
			\begin{itemize}
				\item Any finite-difference operation can be decribed as a convolution operation implying that the specific finite-difference operator has a spectral representation that can be compared with the exact $-ik$ operator
				\item Using Taylor Series for accuracy improvement (= length of operator)
				\item Using Fourier concepts for calculating exact derivatives
			\end{itemize}
			$\Rightarrow$ Fourier method can be interpreted as an infinite order finite-difference scheme		
		\end{minipage}
		\begin{minipage}{0.3\textwidth}
			\begin{figure}
				\includegraphics[width=0.6\textwidth]{\figpath/fig_ps_scheme.png}
			\end{figure}
		\end{minipage}
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{Infinite Order Finite Differences}
		Let us restate the previous result of the partial derivative as an inverse Fourier transform defined as 
		\begin{equation*}
			\begin{split}
				\partial_x f(x)& = \frac{1}{\sqrt{2\pi}} \int_\infty^\infty \partial_x F(k) e^{-ikx}dk \\
 				&= \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty -ik F(k) e^{-ikx}dk
			\end{split}
		\end{equation*}
		Defining the factors in front of the complex amplitude spectrum $F(k)$ of function $f(x)$ as
		\begin{equation*}
			\partial_x f(x) \ = \int_{-\infty}^\infty D(k) F(k) e^{-ikx}dk, \ \ D(k)\ = \ -ik
		\end{equation*}
	\end{frame}
		
% ---------------------------------------------------------

	\begin{frame}{Infinite Order Finite Differences}
		Two functions $d(x)$ and $f(x)$ with complex spectra $D(k)$ and $F(k)$, are thus linked by 
		\begin{equation*}
			\begin{split}
				D(k) \ &= \ \mathscr{F}[d] \\
				F(k) \ &= \ \mathscr{F}[f] \\
				d \ast f  \ &=\ \mathscr{F}^{-1}[D(k) F(k)] 
			\end{split}
		\end{equation*}
		where  $\mathscr{F}$ represents the Fourier transform, and $\ast$ denotes convolution, defined in the continuous case as 
		\begin{equation*}
			(d \ast f) (x) \  \eqdef \   \int_{-\infty}^{\infty}d(x')f(x-x')dx'
		\end{equation*}
		and in the discrete case with vectors $d_i,\ i=0,1,\dots,m$, and   $f_j,\ j=0,1,\dots,n$
		\begin{equation*}
			(d \ast f)_k \  = \   \sum_{i=0}^{m}d_k f_{k-i} \ \ , \ k=0,1,\dots,m+n 
		\end{equation*}
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{Infinite Order Finite Differences}
		$D(k)$ in general is nothing else but a function defined in the spectral domain acting like a filter on the complex spectrum $F(k)$. \\
		The convolution theorem implies that
		\begin{equation*}
			\partial_x f(x) \ = \  \int_{-\infty}^\infty d(x-x')f(x')dx'
		\end{equation*}
		where d(x) is a real function, the spatial representation of spectrum $D(k)$, in other words
		\begin{equation*}
			d(x) \ =  \mathscr{F}^{-1}[D(k)] \ .
		\end{equation*}
	\end{frame}
		
% ---------------------------------------------------------

	\begin{frame}{Infinite Order Finite Differences}
			Limiting the wavenumber domain to the Nyquist wavenumber $k_{max}=\pi/dx$. Thus $D(k)$ becomes
			\begin{equation*}
				D(k) \ = \ ik[H(k+k_{max}) - H(k_{kmax}) ]
			\end{equation*}
			where $H()$ denotes the Heaviside function, and to obtain $d(x)$ we simply have to inverse transform 
			\begin{equation*}
				d(x) \ =  \  \mathscr{F}^{-1} [ ik[H(k+k_{max}) - H(k-k_{kmax}) ] ] \ 
			\end{equation*}
			leading to
			\begin{equation*}
				d(x) \ =  \frac{1}{\pi x^2}[sin(k_{max}x) \ - \ k_{max}x \cos(k_{max}x)]
			\end{equation*}
			\begin{figure}
				\includegraphics[width=0.25\textwidth]{\figpath/fig_heaviside.png}
			\end{figure}
	\end{frame}
		
% ---------------------------------------------------------

	\begin{frame}{Infinite Order Finite Differences}
			The r.h.s. of the Fourier integral is a multiplication of two spectra:
			\begin{itemize}
				\item derivative operator $ik$
				\item boxcar function and its solution is the $sinc$ function of the form $\sin(x)/x$
			\end{itemize}
			If space is discretized according to
			\begin{equation*}
				x_n=n\ dx, \ \ n \ = \ -N,\dots\ ,0, \dots,N
			\end{equation*} 
			In this case the convolution integral becomes a convolution sum
			\begin{equation*}
				\partial_x f(x) \ \approx \  \sum_{n=-N}^{n=N} d_n f(x-ndx)
			\end{equation*}
			where $d_n$ is the difference operator.
	\end{frame}
			
% ---------------------------------------------------------

	\begin{frame}{Infinite Order Finite Differences}
		\begin{minipage}{0.5\textwidth}	
			Inserting the discretization into 
			\begin{equation*}
				\begin{split}
				d(x) \ \ & =  \frac{1}{\pi x^2}[sin(k_{max}x) \ \\
				\ & - \ k_{max}x \cos(k_{max}x)]
				\end{split}
			\end{equation*}
			we obtain analytically the discrete difference operator
			\begin{eqnarray*}
 			d_n & = & \left\{
				\begin{array}{ll}
					0 & \mbox{for}\ \ n =  0\\
					\frac{(-1)^n}{n dx}  & \mbox{for} \ \ n \neq 0 \ .
				\end{array}
				\right . 
			\end{eqnarray*}	
		\end{minipage}	
	\mbox{}		
	\hfill
		\begin{minipage}{0.45\textwidth}
			\begin{figure}
				\includegraphics[width=0.95\textwidth]{\figpath/fig_conv.png}
			\end{figure}
		\end{minipage}
	\end{frame}
				
% ---------------------------------------------------------

	\begin{frame}{Infinite Order Finite Differences}
		\begin{minipage}{0.45\textwidth}			
			Truncated Fourier operator\\
			\begin{figure}
				\includegraphics[width=0.9\textwidth]{\figpath/fig_conv_gauss.png}
			\end{figure}
		\end{minipage}
	\mbox{}		
	\hfill		
		\begin{minipage}{0.45\textwidth}			
			Convolutional difference operators
			\begin{figure}
				\includegraphics[width=0.8\textwidth]{\figpath/fig_conv_x.png}
			\end{figure}
		\end{minipage}
	\end{frame}
					
% ---------------------------------------------------------

	\begin{frame}{Infinite Order Finite Differences}
		\begin{minipage}{0.55\textwidth}
			How can we conveniently compare the accuracy of such operators?\\
			\vfil
			$\Longrightarrow$ The space representation of the exact difference operator $D(k) = -ik$ in the wavenumber domain\\
			\vfil
			Thus, for a finite-difference operator $d_n^{FD}$ we will obtain
			\begin{equation*}
				D^{FD}(k) \ = i\ \tilde{k}_\nu(k) \ =  \ \mathscr{F}[d^{FD}_n]
			\end{equation*}
		\end{minipage}
	\mbox{}		
	\hfill				
		\begin{minipage}{0.4\textwidth}
			\begin{figure}
				\includegraphics[width=0.95\textwidth]{\figpath/fig_conv_k.png}
			\end{figure}
		\end{minipage}
	\end{frame}
	
\section{Chebyshev Pseudospectral Method}
			

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

	\begin{frame}{Chebyshev Polynomials}
		Let us start with the trigonometric relation
		\begin{equation*}
			\cos \left[ (n+1) \phi \right] + \cos \left[ (n-1) \phi \right]
			  \ = \ 2 \cos (\phi) \cos (n\phi) \qquad n \in \mathbb{N}  \ . 
		\end{equation*}
		Inserting $n=0$ leads to a trivial statement. However, for $n\geq1$ we obtain statements like
		\begin{equation*}
			\begin{split}
				\cos (2\phi) \ &=\ 2 \cos^2 (\phi) - 1 \\
				\cos (3\phi) \ &=\ 4 \cos^3 (\phi) - 3 \cos (\phi) \\
				\cos (4\phi) \ &=\ 8 \cos^4 (\phi) - 8 \cos ^2 (\phi) + 1\\
				&\ \, \vdots
			\end{split}
		\end{equation*}
	\end{frame}		
							
% ---------------------------------------------------------

	\begin{frame}{Chebyshev Polynomials}
		\begin{block}{Chebyshev Polynomials}
			\begin{equation*}
				\cos (n \phi) \ \eqqcolon \ T_n ( \cos (\phi)) \ = \ T_n (x) 
			\end{equation*}
			with
			\begin{equation*}
				x \ = \ \cos (\phi) \qquad x \in \left[ -1, 1 \right] \ ,   \qquad n \in \mathbb{N}_0 
			\end{equation*}
		\end{block}
		\vspace{0.5cm}
		$T_n$ being the n-th order Chebyshev polynomial. Furthermore
		 \begin{equation*}
			| T_n (x) | \leqslant 1 \mbox{  for  } \left[ -1, 1 \right], \qquad n \in \mathbb{N}_0 
		\end{equation*}
	\end{frame}
							
% ---------------------------------------------------------

	\begin{frame}{Chebyshev Polynomials}
		\begin{minipage}{0.55\textwidth}
			Finally, we can write down the first polynomials in $x\in [-1,1]$
			\begin{equation*}
				\begin{split}
					T_0 (x) \ &= \ 1 \\
					T_1 (x) \ &= \ x \\
					T_2 (x) \ &= \ 2x^2 - 1 \\
					T_3 (x) \ &= \ 4x^3 - 3x \\
					T_4 (x) \ &= \ 8x^4 - 8x^2 - 1 \\
  	   			     & \ \, \vdots
				\end{split}
			\end{equation*}
		\end{minipage}
	\mbox{}		
	\hfill			
		\begin{minipage}{0.4\textwidth}
			\begin{figure}
				\includegraphics[width=0.95\textwidth]{\figpath/fig_cheby.png}
			\end{figure}
		\end{minipage}
	\end{frame}	
							
% ---------------------------------------------------------

	\begin{frame}{Chebyshev Polynomials}
			A generating function calculates the Chebyshev polynomials of any order $n$
			\begin{equation*}
				T_{n+1}(x) \ = \ 2x T_n(x) \ - \ T_{n-1}(x), \ \ n\geq 1 
			\end{equation*} 
			The extremal values $x_k^{(e)}$ of these polynomials have a very simple form
			\begin{equation*}
				x_k^{(e)} \ = \ \cos ( \frac{k \pi}{n} ) \qquad k \ =  0, 1, 2, \dotsc , n \ 
			\end{equation*}
			\begin{figure}
				\includegraphics[width=0.55\textwidth]{\figpath/fig_cheby_pts.png}
			\end{figure}
	\end{frame}	
							
% ---------------------------------------------------------

	\begin{frame}{Chebyshev Polynomials}
		The Chebyshev polynomials form an orthogonal set with respect to the weighting function $w(x)=1/\sqrt{1-x^2}$. \\
		$\Rightarrow$ Using them as a basis for function interpolation
		\begin{equation*}
			f(x) \ \approx \ g_n(x) \ = \ \frac{1}{2} c_0 T_0 (x) + \sum _{k = 1} ^{n} c_k T_k (x) 
		\end{equation*}
		where f(x) is an arbitrary function in the interval $[-1,1]$, $T_n(x)$ are the Chebyshev polynomials, and $c_k$ are real coefficients.
	\end{frame}
							
% ---------------------------------------------------------

	\begin{frame}{Chebyshev Polynomials}
		By minimizing the least-squares misfit between $f(x)$ and $g_n (x)$, the coefficients $c_k$ can be found
		\begin{equation*}
			c_k \ = \ \frac{2}{\pi} \int \limits _{-1} ^1 f(x) T_k(x) \frac{\mathrm{d}x}{\sqrt{1-x^2}} 
		   \qquad k  =  0, 1, \dotsc, n  
		\end{equation*}
		which - after substituting $x=\cos(\phi)$ - can be written as 
		\begin{equation*}
			c_k \ = \ \frac{1}{\pi} \int \limits _{-\pi} ^\pi f(\cos (\phi)) \cos (k \phi) \mathrm{d}\phi   \qquad k  =  0, 1, \dotsc , n   
		\end{equation*}
		These coefficients turn out to be the Fourier coefficients for the even $2\pi-$periodic function $f(\cos(\phi))$ with $x=\cos(\phi)$.
	\end{frame}
							
% ---------------------------------------------------------

	\begin{frame}{Chebyshev Polynomials}
		The points we need are the extrema of the Chebyshev polynomials (Chebyshev-) Gauss-Lobatto points defined as
		\begin{equation*}
			x_i \ = \ \cos (\frac{\pi}{N} i ) \qquad i  =  0, 1, \dotsc , n \ .
		\end{equation*}
		With these unevenly distributed grid points we can define the discrete Chebyshev transform as follows. The approximating function  is
		\begin{equation*}
			g_n^* (x) \ = \ \frac{1}{2} c_0^* T_0 + \sum _{k=1} ^{n-1} c_k^* T_k (x)  +  \frac{1}{2} c_n^* T_n 
		\end{equation*}
	\end{frame}
							
% ---------------------------------------------------------

	\begin{frame}{Chebyshev Polynomials}			
		with the coefficients defined as
		\begin{equation*}
			\begin{split}
				c_k^* \ &= \ \frac{2}{m} \left[ \frac {1}{2} (f(1) + (-1)^k f(-1)) +  \sum _{j=1} ^{m-1} f_j \cos (\frac{k j \pi}{m}) \right] \\
 					  & k \ = \ 0, 1, \dotsc, n, \qquad n=m 
			\end{split}
		\end{equation*}
		where $f(1)$ and $f(-1)$ are the function values at the interval boundaries and $f_j$ are the values at the collocation points $f(x=\cos(j\pi/m))$. The fundamental property is
		\begin{equation*}
			g_m^* (x_i) \ = \ f(x_i) 
		\end{equation*}
		where $x_i$ are the collocation points.
	\end{frame}
								
% ---------------------------------------------------------

	\begin{frame}{Example}
		When we have a function $f(x) = x^3$ in the interval $[-1,1]$ using the Chebyshev transform, the function $f(x)$
		\begin{itemize}
			\item[1] can be exactly interpolated at the collocation points
			\item[2] converges very rapidly with just a few polynomials
		\end{itemize}
		\begin{figure}
			\includegraphics[width=0.6\textwidth]{\figpath/fig_discrete_cheby.png}
		\end{figure}
	\end{frame}
									
% ---------------------------------------------------------

	\begin{frame}{Chebyshev polynomials}
		\begin{figure}
			\includegraphics[width=0.35\textwidth]{\figpath/fig_cheby_cardinal.png}
		\end{figure}
		Two cardinal functions with Chebyshev polynomials for grid points $i=n/2$ (solid line) and $i=n$ (dashed line) are shown for $n=8$ 
	\end{frame}
									
% ---------------------------------------------------------

\section{Chebyshev Derivatives, Differentiation matrices}	
	
	\begin{frame}{Chebyshev Derivatives, Differentiation matrices}
		A convolution operation can be formulated as a matrix-vector product involving Toeplitz matrices. Defining a derivative matrix $D_{ij}$
		\begin{eqnarray*}
 			D_{ij} & = & \left\{
			\begin{array}{lll}
				& -\frac{2N^2+1}{6} & \text{for \ i \ = j = N} \\[6pt]
				& -\frac{1}{2} \frac{x_i}{1-x_i^2} & \text{for \ i \ = j  \ = 1,2,...,N-1} \\[6pt]
				& \frac{c_i}{c_j} \frac{(-1)^{i+j}}{x_i-x_j} & \text{for \ i $\neq$ j = 0,1,...,N} 
			\end{array}
			\right . 
		\end{eqnarray*}
		where N+1 is the number of Chebyshev collocation points $x_i=\cos(i\pi/N)$, $i=0,\dots,N$ and the $c_i$ are given as   
		\begin{eqnarray*}
			 c_i & = & \left\{
			\begin{array}{ll}
 				2 & \quad \text{for i=0 or N} \\
				1 & \quad \text{otherwise}
			\end{array}
			\right . 
		\end{eqnarray*}
	\end{frame}
									
% ---------------------------------------------------------
	
	\begin{frame}{Chebyshev Derivatives, Differentiation matrices}	
		This differentiation matrix allows us to write the derivative of function $u_i = u(x_i)$ simply as
		\begin{equation*}
			\partial_x u_i \ = \ D_{ij} \ u_j
		\end{equation*}
		where the right-hand side is a matrix-vector product, and the Einstein summation convention applies.
	\end{frame}
										
% ---------------------------------------------------------
	
	\begin{frame}{Chebyshev Derivatives, Differentiation matrices}	
		\begin{minipage}{0.65\textwidth}
			\begin{figure}
				\includegraphics[width=0.65\textwidth]{\figpath/fig_diffmat.png}
			\end{figure}		
		\end{minipage}
	\mbox{}		
	\hfill	
		\begin{minipage}{0.3\textwidth}
			\scriptsize{Illustration of differentiation matrices (n=64). {\bf Top left:} Exact Fourier differentiation matrix for regular grid (full). {\bf Top right:} Exact Chebyshev differentiation matrix for Chebyshev collocation points. Increasing weights at the corners overshadows interior values.  {\bf Bottom Left:} Standard 2-point finite difference operator (banded). {\bf Bottom Right:} Tapered Fourier operator (12-point). Matrix is banded. For illustration purposes the square root of the absolute values are shown.}
		\end{minipage}				
	\end{frame}
											
% ---------------------------------------------------------
	
	\begin{frame}{Chebyshev Derivatives, Differentiation matrices}	
		\begin{minipage}{0.3\textwidth}
			\begin{figure}
				\includegraphics[width=0.75\textwidth]{\figpath/fig_chebyderivative.png}
			\end{figure}
		\end{minipage}
	\mbox{}		
	\hfill		
		\begin{minipage}{0.65\textwidth}	
			By testing the differentiation, we define a function to seismic wavefield calculations as
			\begin{equation*}
				f(x_i)=\sin(2x_i)-\sin(3x_i)+\sin(4x_i)-\sin(10x_i)
			\end{equation*}
			in the interval $x_i\in[-1,1]$, where the discrete points are the Chebyshev collocation points $x_i=\cos(\pi i/n),\ i=0,\dots,n$ given for $n = 63$.
		\end{minipage}
	\end{frame}
											
% ---------------------------------------------------------

\section{Elastic 1D with Chebyshev Method}

	\begin{frame}{Elastic 1D with Chebyshev Method}
		\begin{block}{Elastic 1D wave equation using the standard 3-point operator}
			\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*}
		\end{block}
		\vspace{0.5cm}
		where the lower index $i$ corresponds to the spatial discretization and the upper index $j$ to the discrete time levels. \\
		The displacement field as well as the geophysical parameters like density $\rho_i$ and shear modulus $\mu_i$ are defined on the irregular Chebyshev collocation points.
	\end{frame}		
											
% ---------------------------------------------------------

	\begin{frame}{Example}	
		\begin{minipage}{0.55\textwidth}
			\begin{itemize}
				\item Distance between grid points is 80 times smaller at the boundaries
				\item The time step for a stable simulation requires $c dt/dx \leq \epsilon$			
				\end{itemize}
				$\Rightarrow$ Grid distance near boundary is responsible for the global simulation time step
		\end{minipage}
	\mbox{}		
	\hfill				
		\begin{minipage}{0.4\textwidth}
			\begin{table}
				\begin{tabular}[t]{ll}
					Parameter & Value \\
					\hline
					nx & 200 \\
					c & 3000 m/s \\
					$\rho$ & 2500 $kg/m^3$ \\
					dt & 6 $\times 10^{-8}$ s \\
					$dx_{min}$ & 1.2 $\times 10^{-4}$ m \\
					$dx_{max}$ & 0.015 m \\
					$f_0$ & 100 kHz \\ 
					$\epsilon$ & 1.4 \\
				\end{tabular}
			\end{table}
		\end{minipage}
	\end{frame}
											
% ---------------------------------------------------------

	\begin{frame}{Result}	
		\begin{minipage}{0.55\textwidth}
			\begin{figure}
				\includegraphics[width=0.35\textwidth]{\figpath/fig_cheby_el1d.png}
			\end{figure}
		\end{minipage}
	\mbox{}		
	\hfill			
		\begin{minipage}{0.4\textwidth}
			\tiny
			\begin{mdframed}[style=MyFrame]
{\sffamily
{\tiny
\% Time extrapolation \\
\% \\
\indent for i = 1:nt, \\
\%  \\
\% (...) \\
\% Space derivatives \\
\indent \indent     du=D*u';
\indent \indent     du=mu./rho.*du; \\
\indent \indent     du=D*du; \\
\% (...) \\
\% Time extrapolation
\indent \indent   unew=2*u-uold+du'*dt*dt; \\
\% Source injection \\
\indent \indent    unew=unew+gauss./rho*src(i)*dt*dt; \\
\% remapping \\
\indent \indent   uold=u;	 \\									
\indent \indent  u=unew; \\
\% (...) \\
\indent end
}}
			\end{mdframed}
		\end{minipage}
	\end{frame}	
										
% ---------------------------------------------------------

	\begin{frame}{Result}	
	\small
		\begin{itemize}
			\item To obtain a stable solution we need a very small time step that is only needed at the boundaries. Mathematically the time step scales with $O(N^{-2})$.
			\item In principle we can stretch the spatial grid such that the grid points close to the boundaries are further apart while the grid point distances at the centre remain basically unchanged. If that stretching function is $\xi(x)$ then the derivative of a function $f(x)$ on the stretched grid is defined as 
\begin{equation*}
\partial_x f(x)\  = \ \frac{\partial f}{\partial \xi}   \frac{d\xi}{dx} \ . 
\end{equation*}
		\end{itemize}
	\end{frame}
	
												
% ---------------------------------------------------------

	\begin{frame}{Summary}	
	\scriptsize
		\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.
\item The Chebyshev method is based on the description of spatial fields using Chebyshev polynomials defined in the interval $[-1,1]$ (easily generalized to arbitrary domain sizes). Exact interpolation is possible when the discrete fields are defined at the Chebyshev collocation points given by $x_i=\cos(\pi i/n),\ i=0,\dots,n$. Therefore, the derivatives can also be evaluated exactly and errors accumulate only due to the finite-difference time extrapolation.		
		\end{itemize}
	\end{frame}
												
% ---------------------------------------------------------

	\begin{frame}{Summary}	
	\scriptsize
		\begin{itemize}
\item Because of the grid densification at the boundaries of the Chebyshev collocation points very small time steps are required for stable simulations when $n$ is large. This can be avoided by stretching the grids by a coordinate transformation.
\item A main advantage of the Chebyshev method is an elegant formulation of boundary conditions (free surface or absorbing) through the definition of so-called characteristic variables.
\item Pseudospectral methods have isotropic errors. Therefore they lend themselves to the study of physical anisotropy.
\item The derivative operations of pseudospectral methods are of a global nature. That means every point on a spatial grid contributes to the gradient. While this is the basis for the high precision, it creates problems when implementing pseudospectral algorithms on parallel computers with distributed memory architectures. As communication is usually the bottleneck, efficient and scalable parallelization of pseudospectral methods is difficult.		
		\end{itemize}
	\end{frame}

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

\end{document}