Topics

Leap-frog scheme for the wave equation – Matlab

14 July 2017

The function you will find at the end of this post solves the wave equation with homogeneous boundary conditions on a square domain with some initial Gaussian droplet chosen by the user. Discretisation is done using the leap-frog finite difference scheme.

Wave equation

Let $D = (a,b)^2 \in \mathbb{R}^2$, $a,b \in \mathbb{R}$, be a square domain. Consider the following hyperbolic problem \begin{aligned} %\begin{cases} \partial_{tt} u(\mathbf{x},t) - c^2 \Delta u(\mathbf{x},t) &= 0 && \mathbf{x} \in D, \; t \in (0,T], \\ u(\mathbf{x},0) &= u_0(\mathbf{x}), && \mathbf{x} \in D, \\ \partial_t u(\mathbf{x},0) &= 0 && \mathbf{x} \in D, \\ u(\mathbf{x},t) &= 0 && \mathbf{x} \in \partial\Omega, \; t \in (0,T], \end{aligned} %\end{cases} for a given fixed time $T>0$.

Discretisation

Let $N >0$. We first introduce a constant time step $\Delta t = t_{k+1} - t_k$, $k = 0,\dots,N$, so that $t_k = k \Delta t$ with $t_0 =0$ and $t_{N+1} = T$. Secondly, the spatial mesh is created introducing a set of equally spaced $M+1$ points along both x and y directions \begin{aligned} & a = x_0 < x_1 < \dots < x_M = b, \\ & a = y_0 < y_1 < \dots < y_M = b, \end{aligned} with $\Delta x = x_{i+1} - x_i$ and $\Delta y = y_{j+1} - y_j$, for all $i,j = 0,\dots,M$. The discrete solution at a given point $(x_i,y_j)$ at time $t_k$ is denoted by $u^{k}_{i,j}$. The equation is discretised using central difference approximations: $\frac{u^{k+1}_{i,j} - 2u^{k}_{i,j} + u^{k-1}_{i,j}}{\Delta t^2} - c^2 \left( \frac{u^k_{i+1,j} - 2u^k_{i,j} + u^k_{i-1,j}}{\Delta x^2} + \frac{u^k_{i,j+1} - 2u^k_{i,j} + u^k_{i,j-1}}{\Delta y^2} \right) = 0.$ Since we assumed $\Delta x = \Delta y$, we can rewrite the scheme in this way: $u^{k+1}_{i,j} = 2u^{k}_{i,j} - u^{k-1}_{i,j} + c^2\frac{\Delta t^2}{\Delta x^2} \left( u^k_{i+1,j} + u^k_{i,j+1} - 4u^k_{i,j} + u^k_{i-1,j} + u^k_{i,j-1} \right).$ for all $k = 1,\dots,N$.

Initial values

The above scheme requires the initial discrete values for $u^{1}_{i,j}$ and $u^{0}_{i,j}$. They are recovered using the initial conditions in the problem. In fact, $u^{0}_{i,j}$ is simply given by the known data $u_0(x_i,y_j)$. Regards $u^{1}_{i,j}$, we may use a first-order approximation in time: $\frac{u^{1}_{i,j} - u^0_{i,j}}{\Delta t} \approx \partial_t u((x_i,y_j),0) = 0,$ from which we get $u^{1}_{i,j} = u^0_{i,j}$.

Stability

Lastly, in order for the method to be stable, we need the following condition to be satisfied $\left(c \frac{\Delta t}{\Delta x}\right)^2 + \left(c \frac{\Delta t}{\Delta y}\right)^2 \leq 1,$ which becomes $\Delta t \leq \Delta x / c\sqrt{2}$ in our case.

Below you find a matlab function implementing the scheme presented above.


function waveLeapFrog2D(a,b,T,c,np)
%LEAPFROG2D solves the wave equation using the leapfrog method
% input:
%     a    first point of the 1d interval
%     b    second point of the 1d interval
%     T    final time
%     c    parameter of the pde
%    np    number of Gaussian initial pulses
%
% The function solves the wave equation on a square domain \Omega with
% homogeneous boundary conditions:
%  {  u(x,t)_{tt} - c^2*\nabla{u(x,t)} = 0   in x \in \Omega, t \in (0,T],
%  {  u(x,0)   = u0(x)                       in x \in Omega,
%  {  u_t(x,0) = 0                           in x \in Omega,
%  {  u(x,t)   = 0                           in x \in \partial\Omega, t \in (0,T]
% using the 2nd order finite difference leapfrog scheme; the initial data
% u0(x) is chosen to be np gaussian droplet(s).
%
% EXAMPLE 1: waveLeapFrog2D(-1,1,4,1,5);
% Solves (1) on D=(-1,1)x(-1,1) with c=1 and with u0(x) being 5 initial
% gaussian droplets on D
%
% EXAMPLE 2: waveLeapFrog2D(-1,1,4,1);
% As above, but u0(x) is 1 initial gaussian droplet on D (default)
%
% Copyright (c) 2014, L. Rocchi;

if nargin < 5
% using a single initial gaussian droplet
np = 1;
if nargin < 4
error('too less input data!');
end
end

% Discretisation
csq = c^2;             % parameter in the pde
N   = 101;             % number of nodes along one direction
dx  = (b-a)/(N-1);     % spatial step size
xx  = a:dx:b;          % nodes on x-axis
yy  = xx;              % nodes on y-axis

% In order to get a stable method we should have (cdt/dx)^2 + (cdt/dy)^2 \leq 1,
% which gives dt \leq dx/(c\sqrt(2)) (since dx=dy)
dt = dx/(c*sqrt(2));
t  = 0:dt:T;           % vector of times

% Initialisation
[u0,height] = initialu0(xx,yy,np);

% Leapfrog scheme
unew = zeros(N,N);
u1 = u0 + dt * 0;      % u1 = u0
for k = 2:length(t)
% Internal nodes
for i = 2:N-1
for j = 2:N-1
unew(i,j) = 2*u1(i,j) - u0(i,j) + csq*((dt/dx)^2)*( ...
u1(i+1,j) + u1(i-1,j) - 4*u1(i,j) + u1(i,j+1) + u1(i,j-1) );
end
end
% Plot solution
figure(1)
hi = surf(xx,yy,unew);
axis equal; lighting phong; light; grid on;
xlabel('x'); ylabel('y'); axis([xx(1) xx(end) yy(1) yy(end) -height height]);
title(['time = ',num2str(t(k))],'fontweight','bold','fontsize',16);
set(hi,'EdgeColor','none','FaceColor','c');
set(gca,'Fontsize',16);
pause(0.02)
% Update for ext iteration
u0 = u1;
u1 = unew;
end

end % end function

% Child function
function [u0,ampli] = initialu0(xx,yy,np)
%Initialise u(x,t=0) = u0(x)

[X,Y] = meshgrid(xx,yy);
u0 = zeros(size(X));

% np Gaussian droplets centered at [xcyc(1), xcyc(2)]
sigma = 0.05;  % width
ampli = 0.15;  % amplitude
for i = 1:np
xcyc = xx(1) + (xx(end)-xx(1))*rand(1,2);
u0 = u0 + ampli*exp( -((X-xcyc(1)).^2 + (Y-xcyc(2)).^2)/(sigma^2) );
end

% Plot initial data
figure(1)
hi = surf(xx,yy,u0);
axis equal; lighting phong; light; grid on;
xlabel('x'); ylabel('y'); axis([xx(1) xx(end) yy(1) yy(end) -ampli ampli]);
title('u(x,0) (press any key to start)','Fontsize',16);
set(hi,'EdgeColor','none','FaceColor','c');
set(gca,'Fontsize',16);
pause;

end % end child function


This is a simulation using $4$ initial droplets!