%% Finite Difference (FD) computation loop
% Important: Other parameters may not be mentioned here,
% and you might need to define them.
xm = 80; %km
dx = 1; %km
dt = 0.2; %s
B = 2; %km/s
sd = 4; %s
sp = 60; %km
tmax = 8; %s
U_0 = zeros(xm/dx,1); %previous i-1
U_1 = zeros(xm/dx,1); %current i
U_2 = zeros(xm/dx,1); %future i+1
nsp = sp/dx; % this is where the source will be positioned
for j=0:round(tmax/dt+1)
t = dt*j;
% Calculate U2 based on U1 and U0, see equation (7)
% Note that the t terms are not explicitly used in this part
for i=2:nx-1
U_2(i) = (B*dt/dx)^2*()
end
%% Boundary conditions
U_2(1) = 0.0; % Left side of grid
U_2(nx) = 0.0; % Right side of grid
%% Source function , see equation (8)
if t U_2(nsp) = 1.0 - cos(2*pi*t/sd);
end
%% Loop to set current displacement value to past displacement value
% and future displacement value to current displacement value. This is
% where the t aspect of equation 7 comes in.
for i=1:nx
U_0(i) = U_1(i);
U_1(i) = U_2(i);
end
%% Plots for each time step should be here
% ~~~Complete me!~~~
end
%% Plots to be saved should be here
% ~~~Complete me!~~~
%% Save the needed data to be used in the next program!
% ~~~Complete me!~~~Could you please fix my line
U_2(i) = (B*dt/dx)^2*()
this line is supposed to be the coded version of equation 7 above. That's all I need is equation 7 translated to code. Its based on the U_1, U_2, U_3 stuff above as well as all the parameters mentioned above to. So to clarify please put equation 7 in terms of code, U_1, U_2, U_3 are used as part of it for example the left side of equation 7 is U_2
We will now write a Matlab code to solve the 1D wave propagation equation 1 using the model space parameters and model propagation parameters listed below. A template code is provided later to which you will add several key elements. What is the point of these parameters? After we managed to re-write the partial derivatives (PD) of equation 1 using discrete values, we now want to apply it to a particular problem. In this case, imagine shear waves traveling along a rope, or a plane wave in the Earth. Make sure at this stage to use these exact values. We will see the impact of changing them later in Program 2b-1. According to the finite difference equation 7 , the position u at a future time dt can be obtained if we know the value of u at two previous time steps t and tdt. There are many ways of implementing the finite difference (FD) scheme in computer code. Our approach is to first define the motions of every spatial point on the 1D grid (i.e. along the rope) for an initial time, then calculate motions for those same spatial points at a subsequent time step, and so on. You will need to initialize three displacement arrays, defined as; U0:displacementofpointsonalineattimetdtU1:displacementofpointsonalineattimetU2:displacementofpointsonalineattimet+dt(previous)(current)(future) You can initialize arrays with the zeros command (use the help command. Note that Python will use np.zeros). The lengths of the vectors (U0, U1 and U2) are dependent on the model space dimension x and grid spacing dx. s(t)=1cos(T2t) where T is the source period. The value of the source time function should be added to the displacement field at the source position node at each time step. This function, however, is only added to the wave field from time zero through the source duration, after which propagation continues with no source addition. Boundary conditions Boundary conditions are used to specify what happens at the edges of your grid. In this program we will initially set both ends of our 1D grid to have fixed boundary conditions (i.e. reflecting energy) where: u(0kmposition,alltime)u(80kmposition,alltime)=0=0 This means that the displacement value of u will always be zero at the edges of the grid (imagine a rope with both ends attached to walls). In the second program we write, we will change the boundary conditions and see the effects. f(xh)=f(x)f(x)h+21f(x)h2 respectively. Subtracting (3) from (2) and solving for f(x) gives f(x)=2hf(x+h)f(xh). You can see clearly that this is simply the slope of a line between two points. On the other hand, adding (3) with (2) and solving for f(x) gives f(x)=h2f(x+h)2f(x)+f(xh) The form of equation (5) can be used to approximate the temporal as well as spatial derivatives in (1). Consider that as u is a function of x and t, we can equally have started with f(tn) in (2) and (3) and so derived a similar equation for t and n for (5). Now we can write equation 1 as: dt2u(x,t+dt)2u(x,t)+u(x,tdt)=2dx2u(x+dx,t)2u(x,t)+u(xdx,t), where u(x,t) means position of u at location x and time t. Note that x and t are spatial and temporal variables, dx and dt are temporal and spatial steps, and is the shear wave velocity. Based on equation 6 , we can get the value of u(x,t+dt) by rearanging (6). u(x,t+dt)=(dxdt)2{u(x+dx,t)2u(x,t)+u(xdx,t)}+2u(x,t)u(x,tdt). Take some time to appreciate this result. We started with a statement about position u(x,t), and its 2nd order partial derivatives in time and space. Then we approximated u with Taylor expansions, and re-wrote the partial derivatives in (1) into statements about u and points near u in time and space. Figure 1 illustrates the result