Accuracy of the Born Approximation

In this sequence of problems you will construct an experiment to verify the analytical results on the accuracy of the Born approximation that was developed in lecture.

Resources and Code Setup

These exercises are best solved using a single python script file, which I will call You can download a version of this file here, where there is space for you to code your solutions, without having to worry about the 'glue' code.

This problem is based on the 1D wave operators that you developed as part of the first three computer exercises. If you are not confident that your solutions are 100% correct, please use the provided solutions to computer problems 1-3 (;; Please note, the provided assumes that you are using these files.

I suggest you develop your solutions inside this file, and then run the script in an IPython console. This is accomplished by starting an IPython console, from an command prompt, with the command ipython --pylab.

Then, from the IPython console, you can run your script with the command:

%run cp5_sol

Velocity Model

For this problem, the velocity models \(c_0\) and \(c_\varepsilon\) are constant functions and step functions, respectively

\[ c_0(x) = 1\\ c_\varepsilon(x) = \left\{ \begin{array}{cc}c_0(x) & 0 \le x < \frac{1}{2} \\ c_0(x) + \varepsilon & \frac{1}{2} \le x \le 1\end{array} \right. . \]

The function build_step_model has been provided for you in and is called by

C, C0 = build_step_model(eps, config)

where eps is the scalar \(\varepsilon\).

In [1]:
from cp5_sol import forward_operator, linear_forward_operator, build_step_model

# Setup the problem paramters
config = dict()

# Physical domain parameters
config['x_limits'] = [0.0, 1.0]
config['nx'] = 201
config['dx'] = (config['x_limits'][1] - config['x_limits'][0]) / (config['nx']-1)

# Source position
config['x_s'] = 0.1

# Receiver position
config['x_r'] = 0.15

# Set CFL safety constant
config['alpha'] = 1.0/6.0

# Configure wavelet
config['nu0'] = 10 #Hz

# Determine max epsilon
max_eps = 1.0
config['max_eps'] = max_eps

# Define time step parameters
config['T'] = 1.5 #seconds
config['dt'] = config['alpha'] * config['dx'] / (1 + max_eps)
config['nt'] = int(config['T']/config['dt'])


Problem 1: On paper, sketch the perturbed velocity model \(c_\varepsilon\). Where will you see the reflected wave? The transmitted wave?

Problem 2: Assume that the input wavelet has a characteristic freuqency \(\nu_0=10\text{Hz}\).

From the notes, in the reflection regime, the Born approximation will be accurate if \(\varepsilon \ll 1\) and in the transmission regime, it will be accurate if \(\varepsilon \ll \min(1, \frac{\lambda}{\tilde x})\), where \(\tilde x\) is the distance to the reflector. For the purposes of this exercise, assume that \(a \ll b\) is equivalent to \(a < \alpha b\) for \(\alpha = \frac{1}{4}\).

For \(\varepsilon \in \{0.001, 0.01, 0.04, 0.1, 0.25, 0.5, 0.75,1.0\}\), predict whether the Born approximation will be accurate at each position \(x \in \{0, 0.45, 0.55,1\}\).

Problem 3: For each of the above positions and values of \(\varepsilon\), verify your prediction numerically by checking to see if

\[ \frac{||u_{sc}(\cdot,t) - u_1(\cdot,t)||^2_2}{||u_{sc}(\cdot,t)||^2_2} < \alpha, \]

where the \(\cdot\) symbol is the placeholder for a fixed spatial position.

You may use the provided forward_operator and linear_forward_operator routines, found in, which return the wavefields as an \((n_t, n_x)\) array, for this.

How do these results compare to the theoretical prediction? Is there a better measure of accuracy?

In [2]:
from cp5_sol import print_table

eps = [0.001, 0.01, 0.04, 0.1, 0.25, 0.5, 0.75, 1.0]

xs = [0.0, 0.45, 0.55, 1.0]

alpha = 0.25

lam = 1.0 / config['nu0']

print_table(eps, xs, alpha, lam, config)
eps	x=0.0     	x=0.45    	x=0.55    	x=1.0     
0.001	True/True	True/True	True/True	True/True
0.01	True/True	True/True	True/True	True/True
0.04	True/True	True/True	True/False	True/False
0.1	True/True	True/False	True/False	False/False
0.25	True/True	True/False	True/False	False/False
0.5	False/True	False/False	False/False	False/False
0.75	False/True	False/False	False/False	False/False
1.0	False/True	False/False	False/False	False/False

Problem 4: For \(\varepsilon \in \{0, \frac{\lambda}{x}, 1\}\), use the provided plot_space_time_diagrams(u, u0, u1, config) function to compare the four relevant wavefields. Compare the diagrams to the ones in your notes and use the diagrams to understand where inaccuracy arrises.

Problem 5: For \(\varepsilon \in \{0, \frac{\lambda}{x}, 1\}\), plot the traces of \(u_{sc}\) and \(u_1\) for \(x \in \{0, 0.45, 0.5, 0.55,1\}\) to see indeed how different the wavefields are in each regime. How will these differences impact the optimization/full waveform inversion process?

Are there areas that look fairly accurate but your metric says they are inaccurate? What about the opposite situation?

P4-5 Examples: (\(\varepsilon\), \(x\), theoretically accurate[T/F], numerically accurate[T/F])

In [3]:
from cp5_sol import plot_interesting_pair

A) (0.001, 0.0, T, T)

In [4]:
plot_interesting_pair(0.001, 0.0, True, True, config, ipynb=True)
<matplotlib.figure.Figure at 0x6b5eba8>

B) (0.04, 0.55, T, F)

In [5]:
plot_interesting_pair(0.04, 0.55, True, False, config, ipynb=True)
<matplotlib.figure.Figure at 0x7181358>

C) (0.1, 0.45, T, F)

In [6]:
plot_interesting_pair(0.1, 0.45, True, False, config, ipynb=True)
<matplotlib.figure.Figure at 0x6bfa9e8>

D) (1.0, 0.0, F, T)

In [7]:
plot_interesting_pair(1.0, 0.0, False, True, config, ipynb=True)
<matplotlib.figure.Figure at 0x6da4d68>

E) (1.0, 1.0, F, F)

In [8]:
plot_interesting_pair(1.0, 1.0, False, False, config, ipynb=True)
<matplotlib.figure.Figure at 0x7190cf8>