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!