Solving the wave equation in 1D

From SEG Wiki
Jump to navigation Jump to search

This is the second lab of a series as listed in Seismology labs.
The wave equation describes the propagation of waves and it is a second order partial differential equation (PDE) that involves second derivatives in space and time.

One dimensional wave equation

Differential equation

In the one dimensional wave equation, there is only one independent variable in space. In this case we assume that x is the independent variable in space in the horizontal direction. Additionally, the wave equation also depends on time t. The displacement u=u(t,x) is the solution of the wave equation and it has a single component that depends on the position x and timet. In this case we assume that the motion (displacement) occurs along the vertical direction. We can visualize this solution as a string moving up and down.

The wave equation in one space dimension can be written as follows:

Where c is the speed at which the wave propagates.

For this case we assume that the length of the string is 1, meaning that: .

To find a unique solution, we need to impose initial conditions (IC) and boundary conditions (BC).

Initial conditions (IC)

To impose Initial conditions, we define the solution u at the initial time t=0 for every position x. Since the wave equation has 2 partial derivatives in time, we need to define not only the displacement but also its derivative respect to time.

In this case we assume that both displacement and its derivative respect to time are zero at t=0 for every position x

Boundary conditions (BC)

We also need to define boundary conditions meaning that the solution or its derivative respect to space has to be known at the boundaries of the domain for any time grater than zero. In this case the boundaries are the left x=0 and right x=1 ends of the string.

For this case we are defining displacement at both ends of the string. For the left end, x=0, we define the solution as a function that varies in time, and for the right end, x=1, we assume that this end is fixed and it does not move for any time greater than zero:

Discretization of the wave equation: finite difference (FD)

The wave equation as shown by (eq. 1) is a continuous analytical PDE, in which x can take infinite values between 0 and 1, similarly t can take infinite values greater than zero.

To solve the wave equation by numerical methods, in this case finite difference, we need to take discrete values of x and t: For instance we can take nx points for x and nt points for t, where nx and nt are positive integer numbers. Then the values that x and t take are:

x = {0, Δx, 2Δx,...,(nx-1)Δx};
t = {0, Δt, 2Δt,..., (nt-1)Δt};

where Δx = L/(nx-1) and Δt = T/(nt-1). For this case L=1, and T is the total time of simulation.

Then, when solving the wave equation, we are only solving for the defined points for x and t

Discrete differential equation

To approximate the wave equation (eq. 1), we will use Taylor series expansion. Taylor series is a way to approximate the value of a function at a given point by using the value it takes at a nearby point.

Space Discretization:

In this case we use Taylor series to approximate the displacement uat x +/- Δx knowing its value at x. Assuming that , the Taylor series for the forward and backward expansions for a given time i are:

Forward difference:

Backward difference:

Summing both expansions and finding  :

To simplify the nomenclature, we replace x by the counter j, then x+Δx becomes j+1, meaning the next position and x-Δx becomes j-1, meaning the previous position. Then the above formula can be expressed as:

Time discretization:

Similarly, we use the forward and backward Taylor series expansion to find the second derivative respect to time t for a given position j:

We replace t by the counter i, then t+Δt becomes i+1, meaning the next time and t-Δt becomes i-1, meaning the previous time. Then the above formula can be expressed as:

PDE discretization:

Replacing eq.4 and eq.5 into eq.1 we get:

Recursion formula:

In this case, we are using an explicit method to solve for u at a new time step i+1, knowing the solution for the previous time step i. Then, we can write eq. 6 as :


Discrete initial conditions

Only the IC corresponding to the derivative of u respect to time (eq. 2b) needs to be discretized. We use the forward and backward Taylor series as following:

Substracting the above terms, we get that:

According to eq.2 b when t=0 this IC is equal to zero, then:

Finally we get that:

And using the counter i+1 to replace t+Δt and i-1 to replace t-Δt we can write the above equation as

Discrete boundary conditions

There is no need for further treatment for the BC specified in eq. 3a and 3b.

MATLAB implementation

The overall goal is to find the solution u . Vertical displacement for every discrete position of x at every discrete time t. The propose strategy is to define u as a matrix with as many columns as discrete positions and as many rows as discrete times, and find the corresponding solutions.

Example: Let's define dx as the step in position and dt as the time step. Let's also define nt as the total number of rows and nx as the total number of columns then:

row 1 will correspond to the solution at t=0, 
row 2 at t=dt, 
row 3 at t=2dt 
row i  at t=(i-1)dt
and row nt at t=(nt-1)dt;


column 1  will correspond to the solution at x=0, 
column 2 at x=dx, 
column 3 at x=2dx,
column j at x=(j-1)dx
and column nx at x=(nx-1)dx;

To calculate the solution for every discrete time and position we will use a modified version of (eq. 7)


where i corresponds to the i-th row and j to the j-th column

Input data

MATLAB code:

L=1;  % string length 
c=1;   % wave velocity
f=5;  % frequency of the input wave

time=4;   % total simulation time
dx=0.005;  % step in position
dt=0.005  % time step

Defining the solution matrix

MATLAB code:

nx=ceil(L/dx)+1;  % total number of columns
nt=ceil(time/dt)+1;  % total number of rows
u=zeros(nt,nx);   % initializing matrix u with zeros 

Boundary conditions

Boundary conditions are defined at both ends of the string:

  • Right end BC

According to (eq. 3b) , at x=1, for any time greater than zero t>0, the displacement is equal to zero u =0.

The solution for x=1 corresponds to the last column of our solution matrix u. To enforce this BC we need to set all the rows of this last column to zero, except the first row since the fist column corresponds to t=0. However this step was already done when we initialized the matrix u with zeros.

  • Left end BC

According to (eq. 3a), at x=0 for any time greater than zero t>0, the displacement is equal to a time dependent function u = 3*sin(2*pi*f*t)*exp(-t).

In our matrix u, we need to set the solution for the first column, since this column correspond to x=0. Furthermore, this solution varies with time, meaning that the value of u will be different for each row. We will not set any value for the first row since this row corresponds to t=0.

This BC can be set using a for loop that begins at row 2 till the last row of matrix u. For every row, we need to calculate the corrosponding discrete time t and then evaluate the solution using (eq. 3a).

Initial conditions

We have 2 IC. According to (eq. 2a), the solution is u=0 at t=0 for all x . This means that the first row of the matrix u should be set to zero. However this step has already been done when we defined the solution matrix and initialized it with zeros.

The second IC (eq. 2b) will be implemented directly in the first time step. The solution for first time step, t=dt, corresponds to the second row of matrix u,i=2. Then using (eq. 9) we get:

Notice that in the above equation we have the term u(0,j) at the right hand side. However the counter i in MATLAB is defined from 1 not 0. To fix this issue we will use (eq. 8). This equation is valid at t=0, corresponding to the first row of matrix u, i=1. Then (eq. 8) becomes:

Finally replacing (eq. 11) into (eq. 10) we get :

We implement (eq. 12) in the MATLAB code as follows:


for j=2:nx-1  %loop for every column
  u(2,j) =u(1,j)+1/2*r*(u(1,j+1) - 2*u(1,j)+ u(1,j-1)); % solution for the first time step 

Populating the solution matrix

To find the solution for the remaining time steps, i≥3, we will implement (eq. 9) using 2 for loops, one for the rows (time) and one for the columns (space). Remember that the first and last column correspond to the BC which are already set.

complete the MATLAB code accordingly

for i =  % complete this part. For loop for  time (rows)
    for j= % complete this part. For loop for space (columns)
     % use eq. 9 to complete this part

Notice that the loop for columns j does not include the first and last column because they were considered when setting the BC.


The plotSimulation function simulates the string displacement .

function plotSimulation(u)

close all;
hold on;
for i = 3 : nt
    axis([0 nx -5 5]);
    hold off;