% beamer version of chapter 4 "The Finite Difference 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{amsmath}

\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}

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

% 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_fd}
\newcommand{\figpatha}{c:/mydata/Heiner/Text/Books/OUP/Submission/Igel-ComputationalSeismology/Part2/Figs_fd}

\newcommand{\fs}{\footnotesize}
\def\be{\begin{equation*}}
\def\ee{\end{equation*}}

\title[Computational Seismology] % (optional, use only with long paper titles)
{The Finite Difference 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.55\textwidth}
			\begin{figure}
				\includegraphics[width=0.95\textwidth]{\figpath/fig_cologne.png}
			\end{figure}
		\end{minipage}
    \mbox{ }
    \hfill		
		\begin{minipage}{0.4\textwidth}
			\begin{block}{}
				\begin{itemize}
					\item Simple concept
					\item Robust
					\item Easy to parallelize
					\item Regular grids
					\item Explicit method
				\end{itemize}	
			\end{block}	
		\end{minipage}
	\end{frame}

% ---------------------------------------------------------
\subsection{History}
	
	\begin{frame}{History}
            \begin{itemize}[<+->]
					\item Several Pioneers of solving PDEs with finite-difference method (Lewis Fry Richardson, Richard Southwell, Richard Courant, Kurt Friedrichs, Hans Lewy, Peter Lax and John von Neumann)
					\item First application to elastic wave propagation (Alterman and Karal, 1968)
					\item Simulating Love waves and was the frst showing snapshots of seimsic wave fields (Boore, 1970)
					\item Concept of staggered-grids by solving the problem of rupture propagation (Madariaga, 1976 and Virieux and Madariaga, 1982)			
            \end{itemize}
	\end{frame}	
	\begin{frame}{History}
            \begin{itemize}[<+->]				
					\item Extension to 3D because of parallel computations (Frankel and Vidale, 1992; Olsen and Archuleta, 1996; etc.)
					\item Application to spherical geometry by Igel and Weber, 1995; Chaljub and Tarantola, 1997 and 3D spherical sections by Igel et al., 2002
					\item Incorporation in the first full waveform inversion schemes initially in 2D, e.g. (Crase et al., 1990) and later in 3D (Chen et al., 2007)
            \end{itemize}
	\end{frame}
% ---------------------------------------------------------
\subsection{Finite Differences in a Nutshell}	
	
	\begin{frame}{Finite Differences in a Nutshell}
  	 	 \begin{minipage}{0.55\textwidth}
			\begin{figure}
				\includegraphics[width=0.45\textwidth]{\figpath/fig_fd_nutshell.png}
			\end{figure}
		\end{minipage}
    \mbox{ }
    \hfill
		\begin{minipage}{0.42\textwidth}
			\begin{block}{}	
				\begin{itemize}
					\item[]
					\item Snapshot in space of the pressure field p
					\item Zoom into the wave field with grid points indicated by +
					\item Exact interpolate using Taylor series
					\item[]
				\end{itemize}
			\end{block}	
		\end{minipage}
	\end{frame}
	
% ---------------------------------------------------------
	\begin{frame}{Scalar wave equation}
	\vskip1truecm
		\begin{block}{1D acoustic wave equation}	
			\begin{equation*}
				\ddot{p}(x,t) = c(x)^{2} \; \partial_x^2 p(x,t) + s(x,t)
			\end{equation*}
		\end{block}	
	\vspace*{0.3cm}
		\begin{itemize}
		\item[] p \hspace{0.5cm}  pressure
		\item[] c \hspace{0.5cm}   acoustic velocity
		\item[] s \hspace{0.5cm}   source term
		\end{itemize}
	\vspace*{0.3cm}
		\begin{block}{Approximation with a difference formula}	
			\begin{equation*}
				\ddot{p}(x,t) \approx \frac{p(x, t + dt) - 2p(x,t) + p(x,t - dt)}{dt^2}
			\end{equation*}
		\end{block}
	\vspace*{0.1cm}
		and equivalently for the space derivative		
	\end{frame}
	
% ---------------------------------------------------------

\section{Finite Differences and Taylor Series}


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

	\begin{frame}{Finite Differences}
		\begin{block}{Forward derivative}	
			\begin{equation*}
				d_x f(x) = \underset{dx \: \rightarrow \: 0}{lim} \; \frac{f(x + dx) - f(x)}{dx}
			\end{equation*}
		\end{block}
		\begin{block}{Centered derivative}			
			\begin{equation*}
				d_x f(x) = \underset{dx \: \rightarrow \: 0}{lim} \; \frac{f(x + dx) - f(x - dx)}{2dx}
			\end{equation*}
		\end{block}	
		\begin{block}{Backward derivative}	
			\begin{equation*}
				d_x f(x) = \underset{dx \: \rightarrow \: 0}{lim} \; \frac{f(x) - f(x - dx)}{dx}
			\end{equation*}
		\end{block}	
	\end{frame}

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

	\begin{frame}{Finite Differences}
		\begin{block}{Forward derivative}	
			\begin{equation*}
				d_x f^{+} \approx \frac{f(x + dx) - f(x)}{dx}
			\end{equation*}
		\end{block}
		\begin{block}{Centered derivative}			
			\begin{equation*}
				d_x f^{c} \approx  \frac{f(x + dx) - f(x - dx)}{2dx}
			\end{equation*}
		\end{block}	
		\begin{block}{Backward derivative}	
			\begin{equation*}
				d_x f^{-} \approx  \frac{f(x) - f(x - dx)}{dx}
			\end{equation*}
		\end{block}	
	\end{frame}

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

	\begin{frame}{Finite Differences and Taylor Series}
		\begin{itemize}
			\item[] The approximate sign is important here as the derivatives at point x are not exact. Understanding the accuracy by looking at the definition of Taylor Series:
			\item[]
			\item[] $f(x + dx) = f(x) + f{}' (x) \; dx + \frac{1}{2!} \; f{}''(x) \; dx^2 + O(dx^3)$
			\item[]
			\item[] Subtraction with f(x) and division by dx leads to the definition of the forward derivative:
			\item[]
			\item[] $\frac{f(x + dx) - f(x)}{dx} = f{}' (x) \; + \frac{1}{2!} \; f{}''(x) \; dx + O(dx^2)$
		\end{itemize}			
	\end{frame}

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

	\begin{frame}{Finite Differences and Taylor Series}
		\begin{itemize}
			\item[] Using the same approach - adding the Taylor Series for $f(x + dx)$ and $f(x - dx)$ and dividing by $2dx$ leads to:
			\item[]
			\item[] $\frac{f(x + dx) - f(x - dx)}{2dx} = f{}' (x) + O(dx^2)$
			\item[]
			\item[] This implies a centered finite-difference scheme more rapidly converges to the correct derivative on a regular grid
			\item[] $\Longrightarrow$ It matters which of the approximate formula one chooses
			\item[] $\Longrightarrow$ It does not imply that one or the other finite-difference approximation is always the better one
		\end{itemize}			
	\end{frame}

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

\subsection{Higher Derivatives}

	\begin{frame}{Higher Derivatives}
		\begin{itemize}
			\item[] The partial differential equations have often 2nd (seldom higher) derivatives
			\item[] Developing from first derivatives by mixing a forward and a backward definition yields
			\item[]
		\end{itemize}
		\begin{equation*}
			\partial_x^2 f \approx \frac{\frac{f(x + dx) - f(x)}{dx} - \frac{f(x) - f(x - dx)}{dx}}{dx} = \frac{f(x + dx) - 2 f(x) + f(x - dx)}{dx^2}
		\end{equation*}
	\end{frame}

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

	\begin{frame}{Higher Derivatives - Alternative derivation}
		\begin{itemize}
			\item[] Determining the weights with which the function values have to be multiplied to obtain derivative approximations ...
			\item[]
			\item[] $ \textcolor[rgb]{1,0,0}{a} \ f(x + dx) = \textcolor[rgb]{1,0,0}{a} \left[ f(x) + f{}' (x) \; dx + \frac{1}{2!} \; f{}''(x) \; dx^2 + \ldots \right] $
			\item[] \hspace{0.8cm} $ \textcolor[rgb]{1,0,0}{b} \ f(x) = \textcolor[rgb]{1,0,0}{b} \left[ f(x) \right] $
			\item[] $ \textcolor[rgb]{1,0,0}{c} \ f(x - dx) = \textcolor[rgb]{1,0,0}{c} \left[ f(x) - f{}' (x) \; dx + \frac{1}{2!} \; f{}''(x) \; dx^2 - \ldots \right] $
		\end{itemize}
	\end{frame}
	
	
	
% ---------------------------------------------------------

\begin{frame}{Higher Derivatives - Alternative derivation}
	
	
How to determine a, b, and c?
	
\be
\begin{split}
\textcolor[rgb]{1,0,0}{a} f(x + dx) &+ \textcolor[rgb]{1,0,0}{b} f(x) + \textcolor[rgb]{1,0,0}{c} f(x - dx) \ \approx \\
  f(x) &\left[ \textcolor[rgb]{1,0,0}{a} \ + \ \textcolor[rgb]{1,0,0}{b} \ + \ \textcolor[rgb]{1,0,0}{c} \right] \\
+   dx  f^{\hspace{1pt}'}   &\left[ \textcolor[rgb]{1,0,0}{a} \ \ \ \ \ \ \ \  - \ \textcolor[rgb]{1,0,0}{c} \right] \\
+ \tfrac{1}{2!}  dx^2  f^{\hspace{1pt}''} &\left[ \textcolor[rgb]{1,0,0}{a} \ \ \ \ \ \ \ \  + \ \textcolor[rgb]{1,0,0}{c} \right] \ \ 
\end{split}
\ee
	 
\end{frame}	


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

\begin{frame}{Higher Derivatives - Alternative derivation}
	
	
To obtain a 2nd derivative we require 

\begin{equation*}
\begin{split}
a \ + \ b \ + \  c \ & = \ 0 \\
a \ \ \phantom{+} \ \phantom{b} \ \ - c \ & = \ 0 \\
a \ \ \phantom{+} \ \phantom{b} \ \ + c \ & = \ \frac{2!}{dx^2} \ 
\end{split}
\end{equation*}
	 
\end{frame}	

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

\begin{frame}{Higher Derivatives - Alternative derivation}
	
	
\be
\begin{split}
\left(\begin{array}{rrr}
1 & 1 & 1 \\
1 & 0 & -1 \\
1 & 0 & 1 \\
\end{array}\right)
\left(
\begin{array}{c}
a  \\
b  \\
c \\
\end{array}
\right)
&=
\left(\begin{array}{c}
0  \\
0  \\
\frac{2!}{dx} \\  
\end{array}
\right) \\
\bf{A} \hspace{1.6cm} \bf{w}  \hspace{.4cm}  & =  \hspace{.6cm}\bf{s} 
\end{split}
\label{eq_A}
\ee

\begin{center}
with solution
\end{center}

\be
 \bf{w}  =  \bf{A}^{-1}  \bf{s} 
\ee
	 
\end{frame}	



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

	\begin{frame}{Higher Derivatives - Solution }
\begin{equation*}
\begin{split}
a \ & = \ \frac{1}{dx^2} \\
b \ & = \ - \frac{2}{dx^2} \\
c \ & = \  \frac{1}{dx^2}  \ \ 
\end{split}
\end{equation*}
	\end{frame}

% ---------------------------------------------------------
\subsection{High-Order Operators}
	\begin{frame}{High-Order Operators}
		\begin{itemize}
			\item[] What happens if we extend the \textit{domain of influence} for the derivative(s) of our function f(x)? 
			\item[] Let us search for a 5-point operator for the second derivative
		\end{itemize}
		\begin{eqnarray*}
			f{}'' \approx a f(x + 2dx) + b f(x + dx) + c f(x) + d f(x - dx) + e f(x - 2dx)		
		\end{eqnarray*}
		\begin{equation*}
			\begin{split}
				a \ + \ b \ + \  c \ + d \ + e \ & = \ 0 \\
				2a \ + \ b \ \phantom{+} \ -d \ \ - 2e \ & = \ 0 \\
				4a \ + \ b \ \phantom{+} \ +d \ \ + 4e \ & = \ \frac{1}{2\mathrm{d}x^2} \\
				8a \ + \ b \ \phantom{+} \ -d \ \ - 8e \ & = \ 0 \\
				16a \ + \ b \ \phantom{+} \ +d \ \ +16e \ & = \ 0 
			\end{split}
		\end{equation*}
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{High-Order Operators}
		\begin{itemize}
		\item[] Using matrix inversion we obtain a unique solution
		\end{itemize}
		\begin{equation*}
			\begin{split}
				a \ & = \ -\frac{1}{12\mathrm{d}x^2} \\
				b \ & = \ \ \frac{4}{3\mathrm{d}x^2} \\
				c \ & = \ -\frac{5}{2\mathrm{d}x^2} \\
				d\ & = \  \ \frac{4}{3\mathrm{d}x^2} \\
				e\ & = \ -\frac{1}{12\mathrm{d}x^2} \ .
			\end{split}
		\end{equation*}
		\begin{itemize}
			\item[] with a leading error term for the 2nd derivative is $O(dx^4)$
			\item[] $ \Longrightarrow $ Accuracy improvement
		\end{itemize}
	\end{frame}

% ---------------------------------------------------------
	
	\begin{frame}{High-Order Operators}
		\begin{minipage}{0.5\textwidth}
			\begin{figure}
				\includegraphics[width=0.4\textwidth]{\figpath/fig_taylor.png}
			\end{figure}
		\end{minipage}
    \mbox{ }
    \hfill	
    	\begin{minipage}{0.45\textwidth}
    		\begin{block}{}
    			\begin{itemize}
    				\item[] Graphical illustration of the Taylor Operators for the first derivative for higher orders
    				\item[] The weights rapidly become small for increasing distance to central point of evaluation
    			\end{itemize}
    		\end{block}
		\end{minipage}    			
	\end{frame}
	
% ---------------------------------------------------------

\section{Finite-Difference Approximation of Wave Equations}


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

	\begin{frame}{Acoustic waves in 1D}
		\begin{minipage}{0.7\textwidth}
			\begin{itemize}
				\item[] To solve the wave equation, we start with the simplemost wave equation:
				\item[] The constant density acoustic wave equation in 1D		
			\end{itemize}
			\begin{block}{}
			\begin{equation*}
				\ddot{p} = c^2 \; \partial_x^2 p + s 
			\end{equation*}
			\end{block}
			impossing pressure-free conditions at the two boundaries as
			\begin{equation*}
				p(x) \mid_{x=0,L} = 0
			\end{equation*}
		\end{minipage}
    \mbox{ }
    \hfill	
		\begin{minipage}{0.25\textwidth}
			\begin{figure}
				\includegraphics[width=0.15\textwidth]{\figpath/fig_oboe.jpg}
			\end{figure}					
		\end{minipage}
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{Acoustic waves in 1D}
		The following dependencies apply:
		\begin{alignat*}{4}\label{definitions}
			p  \ &\rightarrow \ &&p(\textbf{x, t}) \qquad &&\text{pressure} \\
			c  \ &\rightarrow \ &&c(\textbf{x}) \qquad &&\text{P-velocity} \\
			s \ &\rightarrow \ &&s(\textbf{x, t}) \qquad &&\text{source term} \\
		\end{alignat*}
		As a first step we need to discretize space and time and we do that with a constant increment that we denote dx and dt.
		\begin{equation*}
			\begin{split}
				x_j \ = \ j \mathrm{d}x, \qquad j \ = \ 0, \ j_{max}  \\
				t_n \ = \ n \mathrm{d}t, \qquad n \ = \ 0, \ n_{max} 
			\end{split}
		\end{equation*}
	\end{frame}

% ---------------------------------------------------------
	
	\begin{frame}{Acoustic waves in 1D}
		Starting from the continuous description of the partial differential equation to a discrete description. The upper index will correspond to the time discretization, the lower index will correspond to the spatial discretization
		\begin{equation*}
			\begin{split}
				p_j^{n+1} \ &\rightarrow \ p(x_j , t_n + \mathrm{d}t ) \\
				p_j^n \ &\rightarrow \ p(x_j , t_n ) \\
				p_j^{n-1} \ &\rightarrow \ p(x_j , t_n - \mathrm{d}t ) \\
				p_{j+1}^n \ &\rightarrow \ p(x_j  + \mathrm{d}x, t_n) \\
				p_j^n \ &\rightarrow \ p(x_j, t_n ) \\
				p_{j-1}^n \ &\rightarrow \ p(x_j  - \mathrm{d}x, t_n) \\  \ .
			\end{split}
		\end{equation*}
	\end{frame}
	
% ---------------------------------------------------------
	
	\begin{frame}{Acoustic waves in 1D}		
		\begin{block}{}
		\begin{equation*}
			\frac{p_j^{n+1} - 2 p_j^n + p_j^{n-1}}{\mathrm{d}t^2} \ = \ c_j^2 \left[ 
   			 \frac{p_{j+1}^n - 2 p_j^n + p^n_{j-1}}{\mathrm{d} x^2} \right] + s^n_j  \ .
		\end{equation*}
		\end{block}
	\begin{minipage}{0.5\textwidth}
		\begin{itemize}
			\item[] the r.h.s. is defined at same time level n
			\item[] the l.h.s. requires information from three different time levels
		\end{itemize}
	\end{minipage}
    \mbox{ }
    \hfill	
	\begin{minipage}{0.45\textwidth}
\begin{figure}
\begin{tikzpicture}[scale = .7]
\draw [->] (0, 0) -- (0,6);
\node [left]  at (0,6) {$t$};
\node [left]  at (0,2) {$n-1$};
\node [left]  at (0,3) {$n$};
\node [left]  at (0,4) {$n+1$};
\draw [->] (0, 0) -- (6,0);
\node [below]  at (2,0) {$j-1$};
\node [below]  at (3,0) {$j$};
\node [below]  at (4,0) {$j+1$};
\node [below]  at (6,0) {$x$};
\draw [fill] (1,1) circle [radius = 0.1];
\draw [fill] (1,2) circle [radius = 0.1];
\draw [fill] (1,3) circle [radius = 0.1];
\draw [fill] (1,4) circle [radius = 0.1];
\draw [fill] (1,5) circle [radius = 0.1];
\draw [fill] (2,1) circle [radius = 0.1];
\draw [fill] (2,2) circle [radius = 0.1];
\draw [fill = gray, gray] (2,3) circle [radius = 0.1];
\draw [fill] (2,4) circle [radius = 0.1];
\draw [fill] (2,5) circle [radius = 0.1];
\draw [fill] (3,1) circle [radius = 0.1];
\draw [fill = gray, gray] (3,2) circle [radius = 0.1];
\draw [fill = gray, gray] (3,3) circle [radius = 0.1];
%\draw [fill = gray, gray] (3,4) circle [radius = 0.1];
\draw  (3,4) circle [radius = 0.1];
\draw [fill] (3,5) circle [radius = 0.1];
\draw [fill] (4,1) circle [radius = 0.1];
\draw [fill] (4,2) circle [radius = 0.1];
\draw [fill = gray, gray] (4,3) circle [radius = 0.1];
\draw [fill] (4,4) circle [radius = 0.1];
\draw [fill] (4,5) circle [radius = 0.1];
\draw [fill] (5,1) circle [radius = 0.1];
\draw [fill] (5,2) circle [radius = 0.1];
\draw [fill] (5,3) circle [radius = 0.1];
\draw [fill] (5,4) circle [radius = 0.1];
\draw [fill] (5,5) circle [radius = 0.1];
\end{tikzpicture}
\end{figure}	
	\end{minipage}
	\end{frame}

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

	\begin{frame}{Acoustic waves in 1D}	
		Assuming that information at time level n (the presence) and $n-1$ (the past) is known, we can solve for the unknown field $p_j^{n+1}$:
		\begin{equation*}
			p_j^{n+1} \ = \ c_j^2 \frac{\mathrm{d}t^2}{\mathrm{d}x^2} \left[ p_{j+1}^n - 2 p_j^n + p_{j-1}^n \right]
			+ 2p_j^n - p_j^{n-1} + \mathrm{d}t^2 s_j^n
		\end{equation*}
		The initial condition of our wave simulation problem is such that everything is at rest at time t = 0:
		\begin{equation*}
			\Eval{p (x,t)}{t = 0}{} \ = \ 0 \ , \ \Eval{\dot{p} (x,t)}{t = 0}{} \ = \ 0 .
		\end{equation*}
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{Acoustic waves in 1D}	
		\begin{itemize}
			\item[] Waves begin to radiate as soon as the source term $s(x,t)$ starts to act
			\item[] For simplicity: the source acts directly at a grid point with index $j_{s}$
			\item[] Temporal behaviour of the source can be calculated by Green's function
			\item[]
		\end{itemize}
		\begin{equation*}
			s(x,t)=\delta(x-x_s) \ \delta(t-t_s)
		\end{equation*}
		\begin{itemize}
		\item[] where $x_{s}$ and $t_{s}$ are source location and source time and $\delta()$ corresponds to the delta function
		\end{itemize}
	\end{frame}

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

	\begin{frame}{Acoustic waves in 1D}	
		A delta function contains all frequencies and we cannot expect that our numerical algorithm is capable of providing accurate solutions
		Operating with a band-limited source-time function:
		\begin{equation*}
			s(x,t)=\delta(x-x_s) \ f(t)
		\end{equation*}
		where the temporal behaviour f(t) is chosen according to our specific physical problem
	\end{frame}

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

	\begin{frame}{Example}	
			Simulating acoustic wave propagation in a 10km column (e.g. the atmosphere) and assume an air sound speed of $c=0.343 m/s$. 
			We would like to $\it hear$ the sound wave so it would need a dominant frequency of at least 20 Hz. 
			For the purpose of this exercise we initialize the source time function $f(t)$ using the first derivative of a Gauss function.
			\begin{equation*}
				f(t) \ = \ -8 \ f_0 \ (t-t_0) \ e^{-\frac{1}{(4f_0)^2} \ (t-t_0)^2}
			\end{equation*}
			where $t_0$ corresponds to the time of the zero-crossing, $f_0$ is the dominant frequency
	\end{frame}

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

	\begin{frame}{Example}	
		\begin{minipage}{0.3\textwidth}
			\includegraphics[width=0.7\textwidth]{\figpath/fig_source.png}		
		\end{minipage}
    \mbox{ }
    \hfill			
		\begin{minipage}{0.6\textwidth}
			\begin{itemize}
			\item What is the minimum spatial wavelength that propagates inside the medium?
           \item  What is the maximum velocity inside the medium?
           \item What is the propagation distance of the wavefield (e.g., in dominant wavelengths)?
			\end{itemize}
		\end{minipage}
	\end{frame}

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

	\begin{frame}{Example}	
		Sufficient to look at the relation between frequency and wavenumber:
		\begin{equation*}
			c \ = \frac{\omega}{k} \ = \ \frac{\lambda}{T} \ = \ \lambda \ f
		\end{equation*}
		where $c$ is velocity, $T$ is period, $\lambda$ is wavelength, $f$ is frequency, and $\omega=2\pi f$ is angular frequency
		\begin{itemize}
			\item[] dominant wavelength of $f_0 = 20 Hz$
			\item[] substantial amount of energy in the wavelet is at frequencies above 20 Hz
			\item[] $\Longrightarrow$ $\lambda = 17m$ and $\lambda = 7m$ for frequencies $20Hz$  and $50Hz$, respectively
		\end{itemize}
	\end{frame}

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


	\begin{frame}[fragile=singleslide]{Python code fragment}		

\begin{minted}[fontfamily=courier, fontsize=\footnotesize]{python}
# Time extrapolation
for it in range(nt):
    # calculate partial derivatives (omit boundaries)
    for i in range(1, nx - 1):
        d2p[i] = (p[i + 1] - 2 * p[i] + p[i - 1]) / dx ** 2
    # Time extrapolation
    pnew = 2 * p - pold + dt ** 2 * c ** 2 * d2p
    # Add source term at isrc
    pnew[isrc] = pnew[isrc] +  dt ** 2 * src[it] / dx
    # Remap time levels
    pold, p = p, pnew
\end{minted}

	\end{frame}
	

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

	\begin{frame}{Result}	
		\begin{minipage}{0.3\textwidth}
			\begin{figure}
				\includegraphics[width=0.55\textwidth]{\figpatha/fig_ac1d.png}
			\end{figure}
		\end{minipage}
    \mbox{ }
    \hfill	
    	\begin{minipage}{0.65\textwidth}
		\begin{itemize}
			\item[] Choosing a grid increment of $dx = 0.5m \longrightarrow$ about 24 points per spatial wavelength for the dominant frequency
			\item[] Setting time increment $dt = 0.0012 \longrightarrow$ around 40 points per dominant period
		\end{itemize}    	
		\end{minipage}    
	\end{frame}

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

	\begin{frame}{Summary}
	\small
		\begin{itemize}
			\item Replacing the partial derivatives by finite differences allows partial differential
equations such as the wave equation to be solved directly for (in principle)
arbitrarily heterogeneous media
			\item The accuracy of finite-difference operators can be improved by using information from more grid points (i.e., longer operators). The weights for the grid points can be obtained using Taylor series
		\end{itemize}
	\end{frame}

\end{document}