GEOL0018 - Numerical methods
Practical week 6 – The Wave Equation
In this practical you will model the wave equation in a 1-D homogeneous spring using a
centered in time – centered in space (CTCS) – scheme, also known as leapfrog scheme.
The spring is L = 1 m long and is made of a material with wave velocity c = 1 m s-1. We will
start by consider homogeneous Dirichlet boundary conditions (i.e., the ends of the spring
are fixed) and will consider separately two sets of initial conditions:
𝜕𝑢
a) 𝑢(0, 𝑥) = sin(2𝜋𝑥) ; (0, 𝑥) = 0
𝜕𝑡
2 𝜕𝑢
b) 𝑢(0, 𝑥) = 𝑒 −400(𝑥−0.3) ; (0, 𝑥) = 0
𝜕𝑡
Recall that the wave equation is:
𝜕2𝑢 2
𝜕2𝑢
=𝑐
𝜕𝑡 2 𝜕𝑥 2
where u is displacement and 𝑥 and 𝑡 are distance and time. Using the short-hand notation
from lectures that:
𝑢𝑖𝑛 = 𝑢(𝑥𝑖 , 𝑡𝑛 )
where xi is the ith space node and 𝑡n is the nth timestep. The CTCS scheme can be written
as:
𝑢𝑖𝑛+1 − 2𝑢𝑖𝑛 + 𝑢𝑖𝑛−1 𝑢 𝑛 − 2𝑢𝑖𝑛 + 𝑢𝑖−1
2 𝑖+1
𝑛
= 𝑐
∆𝑡 2 ∆𝑥 2
1. Re-arrange the CTCS scheme (on paper) into an explicit formulation, with 𝑢𝑖𝑛+1 terms
on the left-hand side and the other terms on the right-hand side. Fully derive both the
scalar and matrix/vector formulations and show all your working by considering a simple
example with unknown points in space.
[See working below]
2. Write a MATLAB script to solve the wave equation using the CTCS scheme. The
MATLAB script geol2029_wave_ctcs_week6_starting.m on Moodle outlines how to
solve the problem for homogeneous Dirichlet BCs. Add the initial conditions a). Program
and calculate the first time step u1,j and plot it, superimposed on the initial condition u0,j .
[See MATLAB script below]
3. Add a new section to step forward in time. For the problem with initial conditions a, plot
the numerical solutions at t=0.01s, t=0.02s, t=0.03s, etc, and compare with the analytical
solution. Discuss the results.
[See MATLAB script below – Section (4) Exact solution of a sinusoidal oscillation:
𝜋 𝜋
𝑢𝑒 (𝑥, 𝑡) = 𝐴sin ( 𝑥) cos( 𝑐𝑡)]
𝐿 𝐿
Numerical solution lags behind analytical solution due to numerical dispersion. This
happens because numerical phase velocity is less than the real/physical phase velocity.
1
GEOL0018 - Numerical methods
4. What happens if you increase the space increment dx? Why does this happen?
Numerical instability – solution explodes / grows exponentially - von Neumann stability
analysis -Courant–Friedrichs–Levy (CFL) condition. We are constrained in our choices of
time and space step.
5. Repeat questions 2 and 3 above, but now for the problem with initial conditions b).
Discuss the physical meaning of the evolution of the displacement u with time.
Gaussian pulse – like plucking a string on a musical instrument.
6. Double the time step (Δt) and re-run the model, plotting the results. What happens?
Why does this happen?
Numerical instability – solution explodes / grows exponentially - von Neumann stability
analysis -Courant–Friedrichs–Levy (CFL) condition. The effects of numerical dispersion
also become greater with time and space step. Therefore, we are constrained in our
choices of time and space step.
7. Change your Matlab script so that you now consider Neumann boundary conditions, i.e.,
the ends of the spring are free, with:
𝜕𝑢 𝜕𝑢
(𝑡, 0) = (𝑡, 𝐿) = 0
𝜕𝑥 𝜕𝑥
along with initial conditions b). What happens? Discuss the results
[See MATLAB script below – Section (4)]
B(1,2) = 2*s2;
B(I,I-1) = 2*s2;
Wave changes polarity when it reflects at the boundaries.
The analytical solution predicts that the shape of the pulse shouldn’t change with time and
should keep bouncing back and forth. But we see some spurious oscillations / ‘ringing’.
This suggests that not all of the energy is propagating at the physical speed, and trails
behind the main wavefront. This is a result of numerical dispersion.
2
GEOL0018 - Numerical methods
Working for Question 1:
3
GEOL0018 - Numerical methods
4
GEOL0018 - Numerical methods
5
GEOL0018 - Numerical methods
MATLAB solution
clear;clc
%%%%%% Solving the wave equation in 1d using the CTCS method %%%%%
%%%%%% (CTCS = centered in time and space)
%%%% (1) DEFINE PARAMETER and CONDITIONS %%%%
v = 1; % Wave-speed (m/s)
L = 1; % Length of string (m)
% boundary conditions (@ = named function)
a = @(t) 0;
b = @(t) 0;
% Define initial conditions
%f = @(x) sin(2*pi*x); % Uncomment for boundary condition (a)
f = @(x) exp(-400*(x-0.3).^2); % Uncomment for boundary condition (b)
g = @(x) 0;
%%%% (2) DISCRETISATION %%%%
% discretisation in space
I = 90; % Number of unknowns (only the interior points)
x = linspace(0,L,I+2); %add boundary points
Dx = L/I;
xint = x(2:end-1); %without boundary points
% Discretistation in time (with CFL criterion)
Dtmax = Dx/v;
Dt = Dtmax/2;
N = 500;
Dt = 0.01;
t = Dt*(0:N-1); % Vector of times
% Matrix B
s2 = (v*Dt/Dx)^2; % sigma-squared
6
GEOL0018 - Numerical methods
B = s2*diag(ones(I-1,1),1) + 2*(1-s2)*diag(ones(I,1),0) + s2*diag(ones(I-1,1),-1); % B matrix with Dirichlet
BC
%B(1,2) = 2*s2; % Uncomment this for B matrix with Neumman BC
%B(I,I-1) = 2*s2; % Uncomment this for B matrix with Neumman BC
%%%% (3) INITIALISATION %%%%
% displacement at first timestep
u = zeros(I,N);
u(:,1) = f(xint);
% Boundary conditions for first time-step
c = zeros(I,1); % for all time-steps
c(1) = s2*a(t(1));
c(I) = s2*b(t(1));
% Displacement at second timestep
u(:,2) = 0.5*B*u(:,1) + Dt*g(xint) + c/2 ;
%% Question (2)
%figure;
%plot(xint, u(:,1), color="blue")
%hold on
%plot(xint, u(:,2), color="red")
%hold off
%%%% (4) COMPUTE REST OF DISPLACEMENTS - main loop (Question 3) %%%%
for n=2:N-1
% Account for varying boundaries
c(1) = s2*a(t(n));
c(I) = s2*b(t(n));
u(:,n+1) = B * u(:,n) - u(:,n-1) + c; % u(:,n+1)=unew (:,n)=u u(:,n-1)=uold
u(1, n+1) = u(2, n+1); % Uncomment for free ends (Neumann BCs) - spatial derivative of the
displacement at the boundaries - ends are free to move vertically but not horizontally
7
GEOL0018 - Numerical methods
u(I, n+1) = u(I-1, n+1); % Uncomment for free ends (Neumann BCs) - spatial derivative of the
displacement at the boundaries - ends are free to move vertically but not horizontally
end
% Plot snapshots
%k = 33;
%figure
%plot(xint, u(:,k), color="k");
% xlabel("Distance (m)")
% ylabel("Amplitude (m)")
%hold on
% Analytical solution for initial condition (a)
%plot(xint, f(xint)*cos(2*v*pi*t(k)), color="red", LineStyle="--")
% Analytical solution for initial condition (b)
%plot(xint, 0.5*f(xint + v*t(k)), color="red", LineStyle="--")
%plot(xint, 0.5*f(xint - v*t(k)), color="red", LineStyle="--")
%hold off
%axis([0 L -1 1])
%legend('Numerical sol.','Analytical sol.')
%title(['Time: ',num2str(t(k))])
%plot solution as animation
figure;
for k=1:1:N
plot(xint, u(:,k), 'b-');
hold on
xlabel("Distance (m)")
ylabel("Amplitude (m)")
% Analytical solution for initial condition (a)
% plot(xint, f(xint)*cos(2*v*pi*t(k)), color="red", LineStyle="--")
% Analytical solution for initial condition (b)
%plot(xint, 0.5*f(xint + v*t(k)), color="red", LineStyle="--")
%plot(xint, 0.5*f(xint - v*t(k)), color="red", LineStyle="--")
8
GEOL0018 - Numerical methods
hold off
axis([0 L -1 1])
legend('Numerical sol.','Analytical sol.')
%ylim([-1,1]);
title(['Time: ',num2str(t(k))])
drawnow;
pause(0.04)
%waitforbuttonpress
end
% Plot some individual snapshots
%figure;
%plot(xint, u(:,160), 'k.-');
%hold on
%plot(xint, 0.5*f(xint - v*t(160)), 'r')
%plot(xint, 0.5*f(xint + v*t(160)), 'r')
%xlabel("Distance (m)")
%ylabel("Amplitude (m)")
