

In [1]:

    
from sympy import *
from sympy.abc import *
from sympy.galgebra.ga import *
import numpy as np
from numpy import linalg as LA
from __future__ import print_function
init_printing()

PDE

The acoustic wave equation for the square slowness m and a source q is given in 3D by :

\begin{cases} &m \frac{d^2 u(x,t)}{dt^2} - \rho \nabla(\frac{1}{\rho} \text{grad}(u(x,t)) =q \\ &u(.,0) = 0 \\ &\frac{d u(x,t)}{dt}|_{t=0} = 0 \end{cases}

with the zero initial conditons to guaranty unicity of the solution



In [3]:

    
p=Function('p')
b=Function('b')
m,s,h,r = symbols('m s h r')
m=M(x,y)
q=Q(x,y,t)
d=D(x,y,t)
e=E(x,y)
r=rho(x,y)

Time and space discretization as a Taylor expansion.

The time discretization is define as a second order ( $ O (dt^2)) $) centered finite difference to get an explicit Euler scheme easy to solve by steping in time.

$ \frac{d^2 u(x,t)}{dt^2} \simeq \frac{u(x,t+dt) - 2 u(x,t) + u(x,t-dt)}{dt^2} + O(dt^2) $

And we define the space discretization also as a Taylor serie, with oder chosen by the user. This can either be a direct expansion of the second derivative bulding the laplacian, or a combination of first oder space derivative. The second option can be a better choice in case you would want to extand the method to more complex wave equations involving first order derivatives in chain only.

$ \frac{d^2 u(x,t)}{dt^2} \simeq \frac{1}{dx^2} \sum_k \alpha_k (u(x+k dx,t)+u(x-k dx,t)) + O(dx^k) $



In [4]:

    
dtt=as_finite_diff(p(x,y,t).diff(t,t), [t-s,t, t+s])
dt=as_finite_diff(p(x,y,t).diff(t), [t-s, t+s])
# Spacial finite differences can easily be extended to higher order by increasing the list of sampling point in the next expression. 
# Be sure to keep this stencil symmetric and everything else in the notebook will follow.
dxx=as_finite_diff(p(x,y,t).diff(x,x), [x-h,x, x+h]) 
dyy=as_finite_diff(p(x,y,t).diff(y,y), [y-h,y, y+h])
dy=as_finite_diff(p(x,y,t).diff(y), [y-h, y+h])
dx=as_finite_diff(p(x,y,t).diff(x), [x-h, x+h])
dyr=as_finite_diff(b(x,y).diff(y), [y-h, y+h])
dxr=as_finite_diff(b(x,y).diff(x), [x-h, x+h])
dtt,dxx,dyy,dt,dx,dy









    Out[4]:





$$\left ( - \frac{2}{s^{2}} p{\left (x,y,t \right )} + \frac{1}{s^{2}} p{\left (x,y,- s + t \right )} + \frac{1}{s^{2}} p{\left (x,y,s + t \right )}, \quad - \frac{2}{h^{2}} p{\left (x,y,t \right )} + \frac{1}{h^{2}} p{\left (- h + x,y,t \right )} + \frac{1}{h^{2}} p{\left (h + x,y,t \right )}, \quad - \frac{2}{h^{2}} p{\left (x,y,t \right )} + \frac{1}{h^{2}} p{\left (x,- h + y,t \right )} + \frac{1}{h^{2}} p{\left (x,h + y,t \right )}, \quad - \frac{1}{2 s} p{\left (x,y,- s + t \right )} + \frac{1}{2 s} p{\left (x,y,s + t \right )}, \quad - \frac{1}{2 h} p{\left (- h + x,y,t \right )} + \frac{1}{2 h} p{\left (h + x,y,t \right )}, \quad - \frac{1}{2 h} p{\left (x,- h + y,t \right )} + \frac{1}{2 h} p{\left (x,h + y,t \right )}\right )$$



In [15]:

    
# In practice, compute rho. grad(rho) upfront (this is a constant) and just compute the following expression, it avoids 
# the recomputation of the gradient at every time step. int the following X and Y are the component of rho.grad(rho) in X and Y respectively

gradgrad=X*dx+Y*dy
expand(gradgrad)









    Out[15]:





$$- \frac{X}{2 h} p{\left (- h + x,y,t \right )} + \frac{X}{2 h} p{\left (h + x,y,t \right )} - \frac{Y}{2 h} p{\left (x,- h + y,t \right )} + \frac{Y}{2 h} p{\left (x,h + y,t \right )}$$



In [14]:

    
lap=(dxx+dyy)
expand(lap + gradgrad)









    Out[14]:





$$- \frac{X}{2 h} p{\left (x,- h + y,t \right )} + \frac{X}{2 h} p{\left (x,h + y,t \right )} - \frac{Y}{2 h} p{\left (- h + x,y,t \right )} + \frac{Y}{2 h} p{\left (h + x,y,t \right )} - \frac{4}{h^{2}} p{\left (x,y,t \right )} + \frac{1}{h^{2}} p{\left (x,- h + y,t \right )} + \frac{1}{h^{2}} p{\left (x,h + y,t \right )} + \frac{1}{h^{2}} p{\left (- h + x,y,t \right )} + \frac{1}{h^{2}} p{\left (h + x,y,t \right )}$$

Solve forward in time

The wave equation with absorbing boundary conditions writes

$ \eta \frac{d u(x,t)}{dt} + m \frac{d^2 u(x,t)}{dt^2} - \rho \nabla(\frac{1}{\rho} \text{grad}(u(x,t)) =q $

and the adjont wave equation

$ -\eta \frac{d u(x,t)}{dt} + m \frac{d^2 u(x,t)}{dt^2} - \rho \nabla(\frac{1}{\rho} \text{grad}(u(x,t)) =q $

where $ \eta$ is a damping factor equal to zero inside the physical domain and decreasing inside the absorbing layer from the pysical domain to the border

And in order to simplify (kinda, I just like it that way) we will rewrite

$ \rho \nabla(\frac{1}{\rho} \text{grad}(u(x,t)) = \nabla^2 u(x,t) + \rho \text{grad}(\frac{1}{\rho}) . \text{grad}(u(x,t)) $



In [16]:

    
# Forward wave equation
wave_equation = m*dtt- (dxx+dyy)- r*(dxr*dx+dyr*dy) - q  + e*dt
stencil = solve(wave_equation,p(x,y,t+s))[0]
ts=lambdify((p(x,y,t-s),p(x-h,y,t), p(x,y,t), p(x+h,y,t),p(x,y-h,t), p(x,y+h,t),b(x-h,y), b(x,y), b(x+h,y),b(x,y-h), b(x,y+h), q , m,r, s, h,e),stencil,"numpy")
stencil









    Out[16]:





$$\frac{1}{2 h^{2} \left(s E{\left (x,y \right )} + 2 M{\left (x,y \right )}\right)} \left(4 h^{2} s^{2} Q{\left (x,y,t \right )} + 2 h^{2} s E{\left (x,y \right )} p{\left (x,y,- s + t \right )} + 8 h^{2} M{\left (x,y \right )} p{\left (x,y,t \right )} - 4 h^{2} M{\left (x,y \right )} p{\left (x,y,- s + t \right )} + s^{2} b{\left (x,- h + y \right )} p{\left (x,- h + y,t \right )} \rho{\left (x,y \right )} - s^{2} b{\left (x,- h + y \right )} p{\left (x,h + y,t \right )} \rho{\left (x,y \right )} - s^{2} b{\left (x,h + y \right )} p{\left (x,- h + y,t \right )} \rho{\left (x,y \right )} + s^{2} b{\left (x,h + y \right )} p{\left (x,h + y,t \right )} \rho{\left (x,y \right )} + s^{2} b{\left (- h + x,y \right )} p{\left (- h + x,y,t \right )} \rho{\left (x,y \right )} - s^{2} b{\left (- h + x,y \right )} p{\left (h + x,y,t \right )} \rho{\left (x,y \right )} - s^{2} b{\left (h + x,y \right )} p{\left (- h + x,y,t \right )} \rho{\left (x,y \right )} + s^{2} b{\left (h + x,y \right )} p{\left (h + x,y,t \right )} \rho{\left (x,y \right )} - 16 s^{2} p{\left (x,y,t \right )} + 4 s^{2} p{\left (x,- h + y,t \right )} + 4 s^{2} p{\left (x,h + y,t \right )} + 4 s^{2} p{\left (- h + x,y,t \right )} + 4 s^{2} p{\left (h + x,y,t \right )}\right)$$



In [ ]:

    
# Adjoint wave equation
wave_equationA =  m*dtt- (dxx+dyy) - r*(dxr*dx+dyr*dy)  - d  - e*dt
stencilA = solve(wave_equationA,p(x,y,t-s))[0]
tsA=lambdify((p(x,y,t+s),p(x-h,y,t), p(x,y,t), p(x+h,y,t),p(x,y-h,t), p(x,y+h,t),b(x-h,y), b(x,y), b(x+h,y),b(x,y-h), b(x,y+h), d, m,r, s, h,e),stencilA,"numpy")
stencilA

Define the discrete model



In [ ]:

    
import matplotlib.pyplot as plt
from matplotlib import animation

hstep=25 #space increment d  = minv/(10*f0);
tstep=2 #time increment dt < .5 * hstep /maxv;
tmin=0.0 #initial time
tmax=600 #simulate until
xmin=-500.0 - 10*hstep #left bound
xmax=500.0 + 10*hstep #right bound...assume packet never reaches boundary
ymin=-500.0 - 10*hstep #left bound
ymax=500.0 + 10*hstep #right bound...assume packet never reaches boundary
f0=.010
t0=1/.010
nbpml=10
nx = int((xmax-xmin)/hstep) + 1 #number of points on x grid
ny = int((ymax-ymin)/hstep) + 1 #number of points on x grid
nt = int((tmax-tmin)/tstep) + 2 #number of points on t grid
xsrc=-400
ysrc=0.0
xrec = nbpml+4
#set source as Ricker wavelet for f0 
def source(x,y,t):
    r = (np.pi*f0*(t-t0))
    val = (1-2.*r**2)*np.exp(-r**2)
    if abs(x-xsrc)<hstep/2 and abs(y-ysrc)<hstep/2:
        return val
    else:
        return 0.0
    
def dampx(x):
    dampcoeff=1.5*np.log(1.0/0.001)/(5.0*hstep);
    if x<nbpml:
        return dampcoeff*((nbpml-x)/nbpml)**2
    elif x>nx-nbpml-1:
        return dampcoeff*((x-nx+nbpml)/nbpml)**2
    else:
        return 0.0
def dampy(y):
    dampcoeff=1.5*np.log(1.0/0.001)/(5.0*hstep);
    if y<nbpml:
        return dampcoeff*((nbpml-y)/nbpml)**2
    elif y>ny-nbpml-1:
        return dampcoeff*((y-ny+nbpml)/nbpml)**2
    else:
        return 0.0



In [ ]:

    
# Velocity and density models
def smooth10(vel,nx,ny):
    out=np.ones((nx,ny))
    out[:,:]=vel[:,:]
    for a in range(5,nx-6):
        out[a,:]=np.sum(vel[a-5:a+5,:], axis=0) /10
    return out

# True velocity
vel=np.ones((nx,ny)) + 2.0
rho=np.ones((nx,ny)) 
vel[floor(nx/2):nx,:]=4.5
rho[floor(nx/2):nx,:]=2.0
rho=rho**-1
mt=vel**-2
# Smooth velocity
v0=smooth10(vel,nx,ny)
m0=v0**-2

dm=m0-mt

Create functions for the PDE

The Gradient/Born are here so that everything is at the correct place, it is described later



In [ ]:

    
def Forward(nt,nx,ny,m,rho):
    u=np.zeros((nt,nx,ny))
    rec=np.zeros((nt,ny-2))
    for ti in range(0,nt):
        for a in range(1,nx-1):
            for b in range(1,ny-1):
                src = source(xmin+a*hstep,ymin+b*hstep,tstep*ti)
                damp=dampx(a)+dampy(b)
                if ti==0:
                    u[ti,a,b]=ts(0,0,0,0,0,0,rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
src,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                elif ti==1:
                    u[ti,a,b]=ts(0,u[ti-1,a-1,b],u[ti-1,a,b],u[ti-1,a+1,b],u[ti-1,a,b-1],u[ti-1,a,b+1],rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
src,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                else:
                    u[ti,a,b]=ts(u[ti-2,a,b],u[ti-1,a-1,b],u[ti-1,a,b],u[ti-1,a+1,b],u[ti-1,a,b-1],u[ti-1,a,b+1],rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
src,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                if a==xrec :
                    rec[ti,b-1]=u[ti,a,b] 
    return rec,u

def Adjoint(nt,nx,ny,m,rho,rec):
    v=np.zeros((nt,nx,ny))
    srca=np.zeros((nt))
    for ti in  range(nt-1, -1, -1):
        for a in range(1,nx-1):
            for b in range(1,ny-1):
                if a==xrec:
                    resid=rec[ti,b-1]
                else:
                    resid=0
                damp=dampx(a)+dampy(b)
                if ti==nt-1:
                    v[ti,a,b]=tsA(0,0,0,0,0,0,rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
resid,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                elif ti==nt-2:
                    v[ti,a,b]=tsA(0,v[ti+1,a-1,b],v[ti+1,a,b],v[ti+1,a+1,b],v[ti+1,a,b-1],v[ti+1,a,b+1],rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
resid,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                else:
                    v[ti,a,b]=tsA(v[ti+2,a,b],v[ti+1,a-1,b],v[ti+1,a,b],v[ti+1,a+1,b],v[ti+1,a,b-1],v[ti+1,a,b+1],rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
resid,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                if abs(xmin+a*hstep-xsrc)<hstep/2 and abs(ymin+b*hstep-ysrc)<hstep/2:
                    srca[ti]=v[ti,a,b]
    return srca,v

def Gradient(nt,nx,ny,m,rho,rec,u):
    v1=np.zeros((nx,ny))
    v2=np.zeros((nx,ny))
    v3=np.zeros((nx,ny))
    grad=np.zeros((nx,ny))
    for ti in range(nt-1,-1,-1):
        for a in range(1,nx-1):
            for b in range(1,ny-1):
                if a==xrec:
                    resid=rec[ti,b-1]
                else:
                    resid=0
                damp=dampx(a)+dampy(b)
                v3[a,b]=tsA(v1[a,b],v2[a-1,b],v2[a,b],v2[a+1,b],v2[a,b-1],v2[a,b+1],rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
resid,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                grad[a,b]=grad[a,b]-(v3[a,b]-2*v2[a,b]+v1[a,b])*(u[ti,a,b])
        v1,v2,v3=v2,v3,v1
    # No update inside the pml, only in the physical domain
    # grad[0:nbpml-1,:]=0
    # grad[nx-nbpml-1:nx-1,:]=0
    # grad[:,0:nbpml-1]=0
    # grad[:,ny-nbpml-1:ny-1]=0
    return tstep**-2*grad

def Born(nt,nx,ny,m,rho,dm):
    u1=np.zeros((nx,ny))
    U1=np.zeros((nx,ny))
    u2=np.zeros((nx,ny))
    U2=np.zeros((nx,ny))
    u3=np.zeros((nx,ny))
    U3=np.zeros((nx,ny))
    rec=np.zeros((nt,ny-2))
    src2=0
    for ti in range(0,nt):
        for a in range(1,nx-1):
            for b in range(1,ny-1):
                damp=dampx(a)+dampy(b)
                src = source(xmin+a*hstep,ymin+b*hstep,tstep*ti)
                u3[a,b]=ts(u1[a,b],u2[a-1,b],u2[a,b],u2[a+1,b],u2[a,b-1],u2[a,b+1],rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
src,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                src2 = -tstep**-2*(u3[a,b]-2*u2[a,b]+u1[a,b])*dm[a,b]
                U3[a,b]=ts(U1[a,b],U2[a-1,b],U2[a,b],U2[a+1,b],U2[a,b-1],U2[a,b+1],rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
src2,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                if a==xrec :
                    rec[ti,b-1]=U3[a,b]
        u1,u2,u3=u2,u3,u1
        U1,U2,U3=U2,U3,U1
    return rec

A Forward propagation example



In [ ]:

    
(rect,ut)=Forward(nt,nx,ny,mt,rho)



In [ ]:

    
fig = plt.figure()
plts = []             # get ready to populate this list the Line artists to be plotted
plt.hold("off")
for i in range(nt):
    r = plt.imshow(ut[i,:,:])   # this is how you'd plot a single line...
    plts.append( [r] )  
ani = animation.ArtistAnimation(fig, plts, interval=50,  repeat = False)   # run the animation
#plt.show()



In [ ]:

    
fig2 = plt.figure()
plt.hold("off")
shotrec = plt.imshow(rect,vmin=-10,vmax=10)   # this is how you'd plot a single line...
#plt.show()

Adjoint test

In ordr to guaranty we have the gradient we need to make sure that the solution of the adjoint wave equation is indeed the true adjoint. Tod os so one should check that

$ <Ax,y> - <x,A^Ty> = 0$

where $A$ is the wave_equation, $A^T$ is wave_equationA and $x,y$ are any random vectors in the range of each operator. This can however be expensive as this two vector would be of size $N * n_t$. To test our operator we will the relax this test by

$ <P_r A P_s^T x,y> - <x,P_SA^TP_r^Ty> = 0$

where $P_r , P_s^T$ are the source and recevier projection operator mapping the source and receiver locations and times onto the full domain. This allow to have only a random source of size $n_t$ at a random postion.



In [ ]:

    
(rec0,u0)=Forward(nt,nx,ny,m0,rho)



In [ ]:

    
(srca,v)=Adjoint(nt,nx,ny,m0,rho,rec0)



In [ ]:

    
plts = []             # get ready to populate this list the Line artists to be plotted
plt.hold("off")
for i in range(0,nt):
    r = plt.imshow(v[i,:,:])   # this is how you'd plot a single line...
    plts.append( [r] )  
ani = animation.ArtistAnimation(fig, plts, interval=50,  repeat = False)   # run the animation
#plt.show()



In [ ]:

    
shotrec = plt.plot(srca)   # this is how you'd plot a single line...
plt.show()



In [ ]:

    
# Actual adjoint test
term1=0
for ti in range(0,nt):
    term1=term1+srca[ti]*source(xsrc,ysrc,(ti)*tstep)

term2=LA.norm(rec0)**2

term1,term2,term1-term2,term1/term2

#if abs(term11/term21-1)<1e-9
#    print('Adjoint test passed')

Least square objective Gradient

We will consider here the least square objective, as this is the one in need of an adjoint. The test that will follow are however necessary for any objective and associated gradient in a optimization framework. The objective function can be written

$ min_m \Phi(m) := \frac{1}{2} \| P_r A^{-1}(m) q - d\|_2^2$

And it's gradient becomes

$ \nabla_m \Phi(m) = - (\frac{dA(m)u}{dm})^T v $

where v is the soltuion if the adjoint wave equation. For the simple acoustic case the gradient can be rewritten as

$ \nabla_m \Phi(m) = - \sum_{t=1}^{nt} \frac{d^2u(t)}{dt^2} v(t) $



In [ ]:

    
# Misfit

F0=.5*LA.norm(rec0-rect)**2
F0



In [ ]:

    
Im1=Gradient(nt,nx,ny,m0,rho,rec0-rect,u0)



In [ ]:

    
shotrec = plt.imshow(rec0-rect,vmin=-100,vmax=100)   # this is how you'd plot a single line...
#plt.show()

shotrec = plt.imshow(Im1,vmin=-100000,vmax=100000)   # this is how you'd plot a single line...
#plt.show()

Adjoint test for the gradient

The adjoint of the FWI Gradient is the Born modelling operator, implementing a double propagation forward in time with a wavefield scaled by the model perturbation for the second propagation

$ J dm = - A^{-1}(\frac{d A^{-1}q}{dt^2}) dm $



In [ ]:

    
Im2=Gradient(nt,nx,ny,m0,rho,rec0-rect,u0)



In [ ]:

    
du1=Born(nt,nx,ny,m0,rho,dm)



In [ ]:

    
term11=np.dot((rec0-rect).reshape(-1),du1.reshape(-1))
term21=np.dot(Im2.reshape(-1),dm.reshape(-1))

term11,term21,term11-term21,term11/term21

#if abs(term11/term21-1)<1e-9
#    print('Adjoint test passed')



In [ ]:

    
shotrec = plt.imshow(Im2,vmin=-10000,vmax=10000)   # this is how you'd plot a single line...
#plt.show()

Jacobian test

The last part is to check that the operators are consistent with the problem. There is then two properties to be satisfied

$ U(m + hdm) = U(m) + \mathcal{O} (h) \\ U(m + h dm) = U(m) + h J[m]dm + \mathcal{O} (h^2) $

which are the linearization conditions for the objective. This is a bit slow to run here but here is the way to test it.

1 - Genrate data for the true model m
2 - Define a smooth initial model $m_0$ and comput the data $d_0$ for this model
3 - You now have $U(m_0)$
4 - Define $ dm = m-m_0$ and $ h = {1,.1,.01,.001,...}$
5 - For each $h$ compute $U(m_0 + h dm)$ by generating data for $m_0 + h dm$ and compute $(J[m_0 + h dm]^T\delta |d) $
6 - Plot in Loglog the two lines of equation above



In [ ]:

    
H=[1,0.1,0.01,.001,0.0001,0.00001,0.000001]
(D1,u0)=Forward(nt,nx,ny,m0,rho)
dub=Born(nt,nx,ny,m0,rho,dm)
error1=np.zeros((7))
error2=np.zeros((7))
for i in range(0,7):
    mloc=m0+H[i]*dm
    (d,u)=Forward(nt,nx,ny,mloc,rho)
    error1[i]  = LA.norm(d - D1,ord=1)
    error2[i]  = LA.norm(d - D1 - H[i]*dub,ord=1)



In [ ]:

    
hh=np.zeros((7))
for i in range(0,7):
    hh[i]=H[i]*H[i]
shotrec = plt.loglog(H,error1,H,H)   # this is how you'd plot a single line...
plt.show()
shotrec = plt.loglog(H,error2,H,hh)   # this is howyou'd plot a single line...
plt.show()

Gradient test

The last part is to check that the operators are consistent with the problem. There is then two properties to be satisfied

$ \Phi(m + hdm) = \Phi(m) + \mathcal{O} (h) \\ \Phi(m + h dm) = \Phi(m) + h (J[m]^T\delta |d)dm + \mathcal{O} (h^2) $

which are the linearization conditions for the objective. This is a bit slow to run here but here is the way to test it.

1 - Genrate data for the true model m
2 - Define a smooth initial model $m_0$ and comput the data $d_0$ for this model
3 - You now have $\Phi(m_0)$
4 - Define $ dm = m-m_0$ and $ h = {1,.1,.01,.001,...}$
5 - For each $h$ compute $\Phi(m_0 + h dm)$ by generating data for $m_0 + h dm$ and compute $(J[m_0 + h dm]^T\delta |d) $
6 - Plot in Loglog the two lines of equation above



In [ ]:

    
(DT,uT)=Forward(nt,nx,ny,mt,rho)
(D1,u0)=Forward(nt,nx,ny,m0,rho)
F0=.5*LA.norm(D1-DT)**2

g=Gradient(nt,nx,ny,m0,rho,D1-DT,u0)
G=np.dot(g.reshape(-1),dm.reshape(-1));
error21=np.zeros((7))
error22=np.zeros((7))
for i in range(0,7):
    mloc=m0+H[i]*dm
    (D,u)=Forward(nt,nx,ny,mloc,rho)
    error21[i]  = .5*LA.norm(D-DT)**2 -F0
    error22[i]  = .5*LA.norm(D-DT)**2 -F0 - H[i]*G



In [ ]:

    
shotrec = plt.loglog(H,error21,H,H)   # this is how you'd plot a single line...
plt.show()
shotrec = plt.loglog(H,abs(error22),H,hh) # this is how you'd plot a single line...
plt.show()



In [ ]:

    
error22



In [ ]:



In [ ]: