Module 1: Setting up the problem

Before we bgin, import SimPEG into ipython notebook as follows:


In [26]:
from SimPEG import *
from IPython.html.widgets import interactive


Efficiency Warning: Interpolation will be slow, use setup.py!

            python setup.py build_ext --inplace
    

Introduction

Every geophysical survey consists of a similar basic framework. An energy source is delivered into the earth, and an array of recievers pick up a signal at the surface...


Module summary:

(1) Start with an expression that relates a kernel function with the continuous distribution of a physical property.
(2) Discretize this expression.
(3) Define a mesh that organizes the data.
(4) Build up the matrix equation $d = Gm$.


Step 1: Physical property distribution and the kernel function.
1.1 define kernel function (...add description)
1.2 define the model. (...add description)

Each datum can be expressed as the inner product of the kernel function and the model:

--the integral needs a physical description--

\begin{equation}d_j = \int_a^b g(x) m(x) dx \end{equation}

...each data point is a measure of a physical property of the entire earth at that point. Consequently it is the sum total (and is therefore the integral) of the influence of the material property specified by the model,


Step 2: Discretize the function for each datum.

\begin{equation}d_i = \sum_{j=1}^N g_i (x_j) m_i \Delta x\end{equation}


We can then gather all the data (the $d^i$s) into a column vector (and say we have $M$ data). Similarly we assemble each kernel function as a row in a matrix, $\widetilde{G}$, and our model parameters $m$ into a column vector of length $N$. For the spacing increments, we require that we obtain a single output for each $\Delta x$ in the kernel function. It follows that this is best represented as a diagonal matrix, with a $\Delta x$ on every diagonal entry. Combining all of the aforemention yields the following vector equation:

\begin{equation} \left[ \begin{array}{c} d_1\\ d_2\\ \vdots\\ d_M \end{array} \right] = \begin{bmatrix} g_1(x_j) & g_1(x_j) & g_1(x_j) & \dots & g_1(x_N) \\ g_2(x_j) & g_2(x_j) & g_2(x_j) & \dots & g_2(x_N) \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ g_M(x_j) & g_M(x_j) & g_M(x_j) & \dots & g_M(x_N) \\ \end{bmatrix} \begin{bmatrix} \Delta x_1 & 0 & \dots &0 \\ 0 & \Delta x_2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \dots & \Delta x_N \\ \end{bmatrix} \left[ \begin{array}{c} m_1\\ m_2\\ \vdots\\ m_N \end{array} \right] \end{equation}


Or put more succinctly: $d = \widetilde{G} diag(\Delta x) m$.

Next, set up "toy problem" with $d_1$, $d_2$ and $j=1:5$.


Step 3: Define a mesh.

While the previous section describes how the data is converted from a continuous function in physical space into discreet data "chunks," what it does not address is the manner in which it is to be represented on a computer so that it is available for processing.
3.1 The physical property of interest lies in the cell centers. (...add sketch)
3.2 The kernel functions reside on the nodes (...add sketch)


Step 4: Build up the matrix equation $d=Gm$.

Note again that the kernel functions "live" on the nodes and the model values in the cell centers. To achieve the inner product of the kernel function with the model parameters we need to refomulate the $G$ matrix to obtain values for the kernel functions at the cell centers.

4.1 Define the matrix $G_n$ (as in "n" for "nodes") as the matrix containing the values for the kernel functions on the nodes. Given that we have $M$ data points and seek $N$ model parameters, it follows that the dimensions of $G_n$ will be $M \times (N+1)$. Schematically, $G_n$ will appear like this:

\begin{equation} G_n = \begin{bmatrix} g_{n_1}^1 & g_{n_2}^1 & g_{n_3}^1 & \dots & g_{n_{N+1}}^1 \\ g_{n_1}^2 & g_{n_2}^2 & g_{n_3}^2 & \dots & g_{n_{N+1}}^2 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ g_{n_1}^M & g_{n_2}^M & g_{n_3}^M & \dots & g_{n_{N+1}}^M \end{bmatrix} \end{equation}

In order to evaluate the kernel functions on the cell centers, we will employ the trapezoidal rule, which states:


In [1]:
import numpy as np

In [2]:
m=5
n=6
w=n-1 
s = (m,n)
A = np.zeros(s)

for i in range(m):
    j=i
    k=i+1
    A[i,j] = 0.5  
    A[i,k] = 0.5
print A


[[ 0.5  0.5  0.   0.   0.   0. ]
 [ 0.   0.5  0.5  0.   0.   0. ]
 [ 0.   0.   0.5  0.5  0.   0. ]
 [ 0.   0.   0.   0.5  0.5  0. ]
 [ 0.   0.   0.   0.   0.5  0.5]]

In [25]:
g = lambda x, k, p, q: np.exp(-p*k*x)*np.cos(2*np.pi*q*k*x)
x=np.linspace(0,1,6)
#x = np.array([0., 0.2, 0.4, 0.6, 0.8, 1.])
p = 0.01
q = 0.1
k = np.array([1, 2, 3, 4, 5, 6])
Gn = np.zeros((len(x), len(k)))

for i in range(len(k)):
    f = g(x,k[i],p,q)
    Gn[:,i] = f
#print f
print Gn


[[ 1.          1.          1.          1.          1.          1.        ]
 [ 0.99013245  0.96471657  0.92421453  0.86932419  0.80096714  0.72027328]
 [ 0.96471657  0.86932419  0.72027328  0.52732179  0.30289805  0.06130149]
 [ 0.92421453  0.72027328  0.41818383  0.06130149 -0.29988416 -0.61488486]
 [ 0.86932419  0.52732179  0.06130149 -0.41237005 -0.77729498 -0.94561804]
 [ 0.80096714  0.30289805 -0.29988416 -0.77729498 -0.95122942 -0.76190351]]

In [ ]:


In [77]:
f = lambda x, y : x + y
 f(1,1)


Out[77]:
2

In [16]:
# make the delta x array
deltax = 0.2*np.ones(m)
V = np.diag(deltax)
print V


[[ 0.2  0.   0. ]
 [ 0.   0.2  0. ]
 [ 0.   0.   0.2]]

In [ ]: