## Solving nth order differential equation for a given boundary condition

We are going to solve the differential equation with the boundary conditions

$\psi_{xx} +[100 - \beta ]\psi = 0$

$\implies \psi_{xx} = [\beta - 100] \psi$

Let’s take the simpler boundary condition

$\psi(x=\pm 1) =0$.

Now, the first thing which we wanna do is to convert this equation in the first order form. After that, we can simply use the 4th order Runge-Kutta to obtain higher order solution.

To transform our 2nd order differential equation into first order, we define two variables. To transform nth order equations we need n such variables. This is the standard process to convert any differential equation into the system of first order differential equations.

Lets take

$y_1 = \psi$ (1)

$\implies y_1' = \psi_x$

$\implies y_1' = y_2$

$y_2= \psi_x$ (2) [derivative of $\psi$]

$\implies y_2' = \psi_{xx}$

$\implies y_2' = [\beta - 100] \psi$

So, we have transformed the 2nd order differential equation into two first order differential equations.

And our boundary condition will become $y_1(\pm 1) =0$

Now for the differential equations of the type

$\psi_{xx} = [\beta - 100] \psi$

$\implies \psi_{xx} = a \psi$

have the solution $\psi = e^{\pm \sqrt{a}x}$ (for positive $a$). Now this solution doesn’t satisfy our condition $\psi(x=\pm 1) =0$. Exponential functions are monotonic and hence cannot be 0 at two points. So, this tells us that $a$ is not positive, i.e., $\beta \ngtr 100$.

Now, if $\beta < 100$, then our eqn becomes

$\psi_{xx} = - a \psi$

The solution to this is of the form of cosines and sines, which has potential to satisfy the boundary condition. So, we have $\beta < 100$.

Defining Function in MATLAB:

function [ rhs ] = yfunc( x,ic,dummy,beta )
y1=ic(1);
y2=ic(2);

rhs = [y2; (beta -100)*y1];
end

Calculating Solutions:

clear; close all; clc
set(0,'defaultlinelinewidth',2);
set(0,'defaultaxeslinewidth',1.5);
set(0,'defaultpatchlinewidth',2);
set(0, 'DefaultAxesFontSize', 14)

tol=10^(-4); % define a tolerance level to be achieved
% by the shooting algorithm

xspan=[-1 1]; %boundary conditions % define the span of the computational domain

A=1; % define the initial slope at x=-1
ic=[0 A]; % initial conditions: x1(-1)=0, x1?(-1)=A

n0=100; % define the parameter n0
dbeta=1;
betasol=[];
beta_start=n0; % beginning value of beta
for jj=1:5 % begin mode loop
beta=beta_start; % initial value of eigenvalue beta
dbeta=n0/100; % default step size in beta
for j=1:1000 % begin convergence loop for beta
[t,y]=ode45('yfunc',xspan, ic,[],beta); % solve ODEs
if y(end,1)*((-1)^(jj+1)) > 0 %this checks if the beta needs to be higher or lower for convergence
beta = beta-dbeta;
else
beta = beta + dbeta/2; %uses bisection to converge to the solution
dbeta=dbeta/2;
end
if abs(y(end,1)) < tol % check for convergence
betasol=[betasol beta]; % write out eigenvalue
break % get out of convergence loop
end
end
beta_start=beta-0.1; % after finding eigenvalue, pick
% new starting value for next mode
norm=trapz(t,y(:,1).*y(:,1)); % calculate the normalization
plot(t,y(:,1)/sqrt(norm)); hold on
end
betasol
legend(sprintf('\\beta_1 = %.4f',betasol(1)),sprintf('\\beta_2 = %.4f',betasol(2)),sprintf('\\beta_3 = %.4f',betasol(3)),...
sprintf('\\beta_4 = %.4f',betasol(4)),sprintf('\\beta_5 = %.4f',betasol(5)));
xlabel('x'); ylabel('\psi (x)','FontWeight','bold')
saveas(gcf,'DiffEqnSol','jpg')

## Introduction to Genetics Algorithm (GA) (Part 2)

To find a basic introduction of GA, the first part can be found here.

III. Examples using Genetics Algorithm

In these examples, we will use Matlab and its function ga to apply GA for the optimization problem. For the manual of using this function, you can find it at https://www.mathworks.com/help/gads/ga.html or type in Matlab:

help ga

III.1. Optimization of a polynomial function

Problem: find the integer x  value that minimizes function:

y = 0.2x2 + 50/x

First, we define the function:

function y = simple_fitness(x)
y = 0.2*x.^2+(50./x);

And use Matlab to finish the job:

clear; close all; clc

x=linspace(0,40,100);
y=simple_fitness(x);
plot(x,y);

rng(1,'twister') % for reproducibility
lb=[1];
[x,fval] = ga(@simple_fitness,1,[],[],[],[],lb

xlabel('Size of cup');
ylabel('Objective Function Value');
title(sprintf('Minimum of the function is %.2f at x = %.2f',min(y), x))
grid on

And here is the result:

III.2. Find the best location of the earthquakes

Problem: we have a set of 5 seismic stations with coordinates: [-2 3 0; 1 3 0; -2 -1 0; 0 -3 0; 2 -2 0]

The travel time of seismic ray from EQ to stations is calculated from:

To simplify the problem, we choose v = 6 km

clear; close all; clc
stnx=[-2 1 -2 0 2];
stny=[3 3 -1 -3 -2];
stnz=[0 0 0 0 0];
EQ=[2 2 -2];

to=sqrt((EQ(1)-stnx).^2+(EQ(2)-stny).^2+(EQ(3)-stnz).^2)/6
to=to+0.05*randn(size(to))
rms = @(mp) sum(sqrt((1./(length(to)-4)).*((sqrt((mp(1)-stnx).^2+(mp(2)-stny).^2+(mp(3)-stnz).^2)/6)-to).^2));
% rms([2 2 -2])
ub=[10 10 0]; %upper bound
lb=[0 0 -10]; %lower bound
options = gaoptimset('Display','iter','PlotFcn',@gaplotbestf,'Generations',150)
% options = gaoptimset('Display','iter','PlotFcn',@gaplotbestf,'TolCon',1e-3);

mp = ga(rms,3,[],[],[],[],lb,ub,[],options)

figure; hold on; grid on;
scatter3(stnx, stny, stnz, 'b', 'filled', '^');
scatter3(mp(1), mp(2), mp(3), 'r', '*');
scatter3(EQ(1), EQ(2), EQ(3), 'g', '*');
view(3); xlabel('X axis'); ylabel('Y axis'); zlabel('Z axis');
legend('Stations','EQ location from GA','Real EQ location');



Utpal Kumar

IESAS

## Simple 1D velocity model inversion from P arrival time

Refer to Chapter 5, Introduction to Seismology, Shearer 2009.

Problem:

From the P-wave travel time data below (note that the reduction velocity of 8km/s), inverse for the 1D velocity model using T(X) curve fitting (fit the T(X) curve with lines, then invert for the ray parameter p and delay time τ(p), then solve for the velocity and depth).

Solution:

We construct the T(X) plot from the data: y = T  – X/8 => T = y + X/8 and then divide the T(X) curve into fitting lines (Fig. 2).

The ray parameter and delay is calculated for each line by using:

p = dT/dX

τ = T – pX

These parameters p and τ are slope and intercept of a line, respectively.

The τ(p) can be written as:

as in the matrix form of:

or, in short,

$\tau = Gh$

We can invert for the thickness of each layer using least square inversion:

$h = (G^{_{T}}G)^{-1}G^{T}\tau$

If we assume ray slowness u(i) = p(i) for each line, we can solve the h and v for the 1D velocity model.

Here is the result of the inversion:

The sharp increase in Vp from ~6.6 km/s to 7.7 km/s at ~ 30km depth suggest for the Moho discontinuity.

The codes for the problem is written in python and can be downloaded: source codeτ inversion code, data points.

Nguyen Cong Nghia, IESAS

## Ray tracing through a 1-D velocity model

Refer to Chapter 4 of Shearer, Introduction to Seismology.

For a ray piercing through Earth, the ray parameter (or horizontal slowness) p is defined by several expressions:

where u = 1/v is the slowness, θ is the ray incidence angle, T is the travel time, X is the horizontal range and utp is the slowness at the ray turning point.

The vertical slowness is defined as:

and integral expressions for the surface-to-surface travel time are:

With these equations we can calculate the travel time (T) and horizontal distance (X) for a given ray with ray parameter p and velocity model v(z).

Apply these to a problem (Exercise 4.8 Shearer):

(COMPUTER) Consider MARMOD, a velocity-versus-depth model, which is typical of much of the oceanic crust (Table 4.1). Linear velocity gradients are assumed to exist at intermediate depths in the model; for example, the P velocity at 3.75 km is 6.9 km/s. Write a computer program to trace rays through this model and produce a P-wave T(X) curve, using 100 values of the ray parameter p equally spaced between 0.1236 and 0.2217 s/km. You will find it helpful to use subroutine LAYERXT (provided in Fortran in Appendix D and in the supplemental web material as a Matlab script), which gives dx and dt as a function of p for layers with linear velocity gradients. Your program will involve an outer loop over ray parameter and an inner loop over depth in the model. For each ray, set x and t to zero and then, starting with the surface layer and proceeding downward, sum the contributions, dx and dt, from LAYERXT for each layer until the ray turns. This will give x and t for the ray from the surface to the turning point. Multiply by two to obtain the total surface-to-surface values of X(p) and T(p). Now produce plots of: (a) T(X) plotted with a reduction velocity of 8 km/s, (b) X(p), and (c) τ(p). On each plot, label the prograde and retrograde branches. Where might one anticipate that the largest amplitudes will occur?

Using the LAYERXT subroutine (Appendix D of Shearer), the FORTRAN and MATLAB can be downloaded from here. In this example, we will use python to calculate and plot the result.

The main program is as follow:

from ttcal import layerxt
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

p1 = 0.1236
p2 = 0.2217
n = 100

result = [['p','X','T']]
for p in np.linspace(p1, p2, n):
X = 0; T = 0;
for i in range(0, len(v_model)-1):
p = p
h = v_model.h[i+1] - v_model.h[i]
utop = 1/v_model.vp[i]
ubot = 1/v_model.vp[i+1]
dx, dt, irtr = layerxt(p,h,utop,ubot)
X += dx; T += dt
#If the ray reflected at layer the dx, dt calculation stops
if irtr == 2:
break
result.append([p,X,T])

result_df = pd.DataFrame(result, columns=result.pop(0))
result_df = result_df[::-1]

#Multiply by 2 to get total surface-to-surface value ot X(p) and T(p)
result_df['T'] = 2*result_df['T']
result_df['X'] = 2*result_df['X']

plt.figure(1)
result_df['RT'] = result_df['T']-result_df['X']/8
result_df.plot(x='X', y='RT', legend=False)
plt.title('Reduction time')
plt.ylabel('Reduce time T - X/8 (s)')
plt.xlabel('Distance (km)')
plt.savefig('ex4.8.1.png')

plt.figure(2)
result_df.plot(x='p', y='X', legend=False)
plt.title('X(p)')
plt.ylabel('Distance X (km)')
plt.xlabel('Ray parameter (km/s)')
plt.xlim([p1,p2])
plt.savefig('ex4.8.2.png')

plt.figure(3)
result_df['Taup'] = result_df['T'] - result_df['p']*result_df['X']
result_df.plot(x='p', y='Taup', legend=False)
plt.title('Tau(p)')
plt.ylabel('Tau(p) (s)')
plt.xlabel('Ray parameter (km/s)')
plt.xlim([p1,p2])
plt.savefig('ex4.8.3.png')

plt.figure(4)
v_model.h = v_model.h * (-1)
v_model.plot(x='vp', y='h', legend=False)
plt.ylabel('Depth (km)')
plt.xlabel('Velocity (km/s)')
plt.xlim([4,9])
plt.ylim([-10,0])
plt.savefig('ex4.8.4.png')

#plt.show()

LAYERXT subroutine written in python:

from math import sqrt, log

def layerxt(p,h,utop,ubot):
''' LAYERXT calculates dx and dt for ray in layer with linear velocity gradient

Inputs:   p     =  horizontal slowness
h     =  layer thickness
utop  =  slowness at top of layer
ubot  =  slowness at bottom of layer
Returns:  dx    =  range offset
dt    =  travel time
irtr  =  return code
= -1,  zero thickness layer
=  0,  ray turned above layer
=  1,  ray passed through layer
=  2,  ray turned within layer, 1 segment counted in dx,dt'''
# ray turned above layer
if p >= utop:
dx = 0
dt = 0
irtr = 0
return dx,dt,irtr

#Zero thickness layer
elif h == 0:
dx = 0
dt = 0
irtr = -1
return dx,dt,irtr

#Calculate some parameters
u1 = utop
u2 = ubot
v1 = 1.0/u1
v2 = 1.0/u2
b = (v2 - v1)/h
eta1 = sqrt(u1**2 - p**2)

#Constant velocity layer
if b == 0:
dx = h*p/eta1
dt = h*u1**2/eta1
irtr = 1
return dx,dt,irtr

x1 = eta1/(u1*b*p)
tau1=(log((u1+eta1)/p)-eta1/u1)/b

#Ray turned within layer, no contribution to integral from bottom point
if p >= ubot:
dx = x1
dtau = tau1
dt = dtau + p*dx
irtr = 2
return dx,dt,irtr

irtr=1
eta2=sqrt(u2**2-p**2)
x2=eta2/(u2*b*p)
tau2=(log((u2+eta2)/p)-eta2/u2)/b

dx=x1-x2
dtau=tau1-tau2

dt=dtau+p*dx
return dx,dt,irtr

def flat(zsph,vsph):
''' Calculate the flatten transformation from a spherical Earth.
zf, vf = flat(zsph,vsph)'''

a = 6371.0

rsph = a - zsph
zf = -a*log(rsph/float(a))
vf = (a/float(rsph))*vsph
return zf, vf

Input velocity model:

h,vp,vs,den
0.0,4.5,2.4,2
1.5,6.8,3.75,2.8
6.0,7.0,3.5,2.9
6.5,8.0,4.6,3.1
10.0,8.1,4.7,3.1

The codes can be downloaded here: main program, python LAYERXT subroutine, velocity model.

You need to install numpy, pandas and matplotlib package to run this. Normally it can be installed in the terminal by (presume python has been installed):

pip install numpy pandas matplotlib

Nguyen Cong Nghia – IESAS

## Modeling a wave on a string using finite differences

Based on problem 3.7 chapter 3 of Introduction to Seismology (Shearer)

(COMPUTER) In the case of plane-wave propagation in the x direction within a uniform medium, the homogeneous momentum equation (3.9) for shear waves can be expressed as

,where u is the displacement. Write a computer program that uses finite differences to solve this equation for a bar 100 km in length, assuming β = 4 km/s. Use dx = 1 km for the length spacing and dt = 0.1 s for the time spacing. Assume a source-time function  at u (50 km) of the form:

Apply a stress-free boundary condition at u(0 km) and a fixed boundary condition at u (100 km).  Approximate the second derivatives using the finite difference scheme:

Plot u(x) at 4 s intervals from 1 to 33 s. Verify that the pulses travel at velocities of 4 km/s. What happens to the reflected pulse at each endpoint? What happens when the pulses cross?

!**********************************************************************
!* *
!* 1-D Plane wave propagation with fixed and free boundaries *
!* *
!**********************************************************************

!----initialzie t, dx, dt, tlen, beta, and u1, u2, u3 arrays----c

implicit none
integer:: i, counter, nx, nx1
integer:: it, ntmax
real(8):: t, dt, dt2, dx, dx2, beta, beta2, tlen, tmax, rhs
real*8, allocatable:: u1(:), u2(:), u3(:)
character(10):: data
character(4):: nf
!----parameters ----c
! dt must be less than (dx/beta) for numerical stability

nx = 100 ! number of grid
dt = 0.10d0 ! time interval
dx = 1.0d0 ! grid interval
beta = 4.0d0 ! wave velocity
tlen = 5.0d0 ! time length of wave source
ntmax = 1000 ! maximum time step (finish at it = ntmax, or)
tmax = 33.0d0 ! maximum calculation time (finish at t< tmax)

!---- allocate dimension variables ----c

nx1 = nx + 1
allocate (u1(nx1), u2(nx1), u3(nx1))

!---- initialize t, counter, u1, u2, u3 ----c

t = 0.0d0
counter = 1

do i = 1, nx1
u1(i) = 0.0d0
u2(i) = 0.0d0
u3(i) = 0.0d0
end do

!---- calculate c^2, dt^2, dx^2 ----c
beta2 = beta**2
dt2 = dt**2
dx2 = dx**2

! ============= Time marching loop ===========

do it = 1, ntmax

t = t + dt

!---- calculate u3(i) ----c

do i= 2, nx
rhs=beta2*2*(u2(i+1)-2.0d0*u2(i)+u2(i-1))/dx2
u3(i)=dt2*rhs+2.0d0*u2(i)-u1(i)
end do

!---- free boundary (du/dx=0) ----c
u3(1)=u3(2)
!---- fixed boundary (u=0) ----c
u3(nx1)=0.0d0
!---- source time function c----c

if (t.le.tlen) then
u3(51) = sin(3.1415927*t/tlen)**2
end if

!---- change new and old variables ----c

do i=1,nx1
u1(i) = u2(i)
u2(i) = u3(i)
end do
!---- make data file ----c
write(nf,'(i3.3)') counter
data = "data"//nf

open(7,file=data,status='replace')
do i=1, nx1
write(7,*) i, u2(i)
end do
!---- output u2 at desired intervals, stop when t is big enough

write(6,'(a5,1x,i3,1x,f6.3)') 'loop', counter, t

if (t.gt.tmax) exit

counter = counter + 1

end do

stop
end

The following is used to plot the results:

#!/usr/bin/bash
gnuplot << eof
set term gif animate
set output "animate.gif"
n=9 #n frames
set xrange [0:102]
set yrange [-2:2]
i=0
list = system('ls data*')
do for [file in list] {
plot file w lines lt 1 lw 1.5 t sprintf("t=%i sec",i/10)
i = i + 1
}
set output
eof

Nguyen Cong Nghia – IESAS

## Calling SAC(Seismic Analysis Code) (in Perl)

For seismologists, using a SAC for sac data manipulation is essential (though there are few alternatives). Here, we see how can we call SAC from a perl script:

#!/usr/bin/perl
open(SAC, "| sac ") or die "Error opening sac\n";
print SAC qq[
echo on
*fg seismogram    #sample seismic signal in SAC's memory
fg sine 2 npts 2000 delta 0.01
*fg impulse npts 100 delta 0.01
bandpass bessel corner 0.1 0.3 npole 4
ppk
fft
plotsp am
save sine_fft.pdf
];
print SAC "quit\n";
close(SAC);

Example script to call SAC functions in perl

## Understanding Seismograms

Seismograms are basic information about earthquakes, chemical and nuclear explosions, mining induced earthquakes, rock bursts and other events generating seismic waves.

Seismograms reflect the combined influence of the seismic source, the propagation paths, the frequency response of the recording instrument and the ambient noise at the recording site.

u(t) = s(t)*g(t)*i(t)*n(t).         “*” means convolution.

where, u is the seismic record, s is the source effect, g is the propagation effect, i is the instrument effect and n is the noise.

So, by understanding seismogram we can understand the seismic source, earth structure, and the noise in the medium.

This is our short attempt to understand the seismograms.

Understanding Seismograms 1

For understanding the seismogram, there are two main things we should notice- record duration and the dispersion. Due to different nature of the propagation velocity of the seismic waves and the different propagation paths taken by them to the station, travel time differences between the main group usually grows with distance. Since the body wave groups do not disperses, so their individual duration remains more or less constant, only the time difference between them changes with distance.

Figure 1: Event-Station Distribution map for the event in Mindanao, Phillipines Island (2005-02-05, Mw-7). The source is 540 km deep.

Figure 2

Figure 3

Figure 4

The epicentral distance for QIZ, BJT and DGAR is ~ 19, 35, 52 degrees respectively. For the station BJT (Figure 3), the epicentral distance is 35.2 degrees and for the station DGAR (Figure 4), the epicentral distance is 52.4 degrees. The P-S arrival times for the station BJT is nearly 300 sec and for the station DGAR, the P-S arrival time is nearly 400 sec. For the station QIZ, it is only around 200 sec. So, with the increase in epicentral distance, the difference in individual phase arrival increases. Also for the station DGAR, the P arrival is later than the BJT.

The time difference between main body wave onsets is:

 Distance (degrees) less than Time-difference in body wave onsets (P-S) in sec 10 180 60 (Fig 1,2,3) 960 100 1800 180 2700

Figure 5 : Travel-Time curve for depth 540 km

For the deep events, we do not see the first arrival as P wave (Figure 6) for short distances (D

Figure 6

Figure 6, is the seismogram at 2 degrees epicentral distance from the deep source. Here, we can notice that the first arrival is not P but the core reflected phase PcP. It is also cevident from the travel time curve(Figure 5).

Figure 7: Event-Station Distribution Map of Near Coast of Peru (2001-06-23, Mw-8.3). The source is at 2 km depth.

When the source is near the surface then we observe the surface waves as well in the seismogram.The surface waves from earthquakes at intermediate (> 70 km) or great depth (> 300 km) may have amplitudes smaller than those of body waves or may not even be detected on seismic records.

In contrast to the body waves, the velocity of surface waves is frequency dependent and hence it is dispersed. The duration of Love and Rayleigh waves increases with distance. As for station NNA (Figure 8), the duration of the first group of Rayleigh wave is about 230 seconds and for station HKT (Figure 9), the duration is around 600 seconds.

We can also notice that for shallower earthquake event, P-wave is the first arrival for stations less than 10 degrees of the source (Figure 10).

Figure 8

In Figure 8, we can notice that for surface waves, the largers periods arrive first than the lowers periods. This is the general case for normal layering. The larger periods sample the higher depths and since the velocity for higher depths is larger so it arrives before the shorter period waves.

Figure 9

In figure 9, for the surface waves, we can observe the wave groups or the beating effect (analogous to acoustics). This is because of the interference of the two or more harmonic waves with closer frequencies or periods.

Figure 10: Travel time curve for depth 2.2km

For certain distance ranges, the travel time curves for different types of seismic phases are close to each other or can even overlap (Figure 5,11).

(a)

(b)

Figure 11: S, SKS and ScS are very close to each other.

The amplitude of P-wave depends not only on the epicentral distance but on the site effect (underground geology and crustal heterogeneity), azimuth etc.

The most obvious indication on a seismogram that a large earthquake has a deep focus is the small amplitude of the surface waves with respect to the body-wave amplitudes and the P and S waveforms often have impulsive onsets.

A more precise determination of the depth h of a seismic source, however, requires either the availability of a seismic network with at least one station being very near to the source, e.g., at an epicentral distance D < h (because only in the near range the travel time t(D, h) of the direct P wave varies strongly with source depth h), or the identification of seismic depthphases on the seismic record.

At distant seismograph stations, the depth phases pP or sP follow the direct P wave by a time interval that changes only slowly with distance but rapidly with depth.