Computational Seismology
Numerical first derivative

This notebook is part of the supplementary material to Computational Seismology: A Practical Introduction, Oxford University Press, 2016.

### Excercise:¶

Initialise arbitrary analytical function and calculate numerical first derivative with 2-point operator, compare with analytical solution. Task: extend scheme to 4 point operator and demonstrate improvement (see Chapter 4)

This exercise covers the following aspects:

• Calculation of numerical first derivative
• Comparison with analytical solution

In [ ]:
# Import Libraries
%matplotlib notebook
import numpy as np
import math
import matplotlib.pyplot as plt

plt.style.use("ggplot")


We initialise a Gaussian function

$$f(t)=\dfrac{1}{\sqrt{2 \pi a}}e^{-\dfrac{t^2}{2a}}$$

In [ ]:
# Initial parameters
tmax=10.0                       # maximum time
nt=100                       # number of time sample
a=1                           # half-width
dt=tmax/nt                    # defining dt
t0 = tmax/2                   # defining t0

time=np.linspace(0,tmax,nt)   # defining time

# Initialization gaussian function with zeros
f=(1./np.sqrt(2*np.pi*a))*np.exp(-(((time-t0)**2)/(2*a)))  # eq.(4.32) p. 80

# Plotting of gaussian
plt.figure()
plt.plot(time, f)
plt.title('Gaussian function')
plt.xlabel('Time, s')
plt.ylabel('Amplitude')
plt.xlim((0, tmax))
plt.grid()
plt.show()


In the cell below we calculate numerical derivative using two points

$$f^{\prime}(t)=\dfrac{f(t+dt)-f(t-dt)}{2dt}$$

and analytical derivative $$f^{\prime}(t)=-\dfrac{t}{a}\dfrac{1}{\sqrt{2\pi a}}e^{-\dfrac{t^2}{2a}}$$

In [ ]:
# First derivative with two points

# Initiation of numerical and analytical derivatives
nder=np.zeros(nt)          # numerical derivative

# Numerical derivative of the given function
for it in range (1, nt-1):
nder[it]=(f[it+1]-f[it-1])/(2*dt)

# Analytical derivative of the given function

# Plot of the first derivative and analytical derivative
plt.figure()
plt.plot (time, nder,label="Numerical derivative", lw=2, color="yellow")
plt.plot (time, ader, label="Analytical derivative", ls="--",lw=2, color="red")
plt.title('First derivative')
plt.xlabel('Time, s')
plt.ylabel('Amplitude')
plt.legend()
plt.grid()
plt.show()


In the cell below calculation of the first derivative with four points is provided with the following weights:

$$f^{\prime}(t)=\dfrac{\dfrac{1}{12}f(t-2dt)-\dfrac{2}{3}f(t-dt)+\dfrac{2}{3}f(t+dt)-\dfrac{1}{12}f(t+2dt)}{dt}$$
In [ ]:
# First derivative with four points

# Initiation of derivative
ffder=np.zeros(nt)

for it in range (2, nt-2):
#ffder[it]=

# Plotting
plt.figure()
plt.plot (time, nder,label="Derivative, 2 points", lw=2, color="violet")
plt.plot (time, ffder, label="Derivative, 4 points", lw=2, ls="--")
plt.title('First derivative')
plt.xlabel('Time, s')
plt.ylabel('Amplitude')
plt.legend()
plt.grid()
plt.show()

In [ ]: