Partial Differential Equation Toolbox

Wave Equation on a Square Domain

This example shows how to solve the wave equation using the hyperbolic function in the Partial Differential Equation Toolbox™.

We solve the standard second-order wave equation

$$ \frac{\partial^2 u}{\partial t^2} - \Delta u = 0$$

on a square domain with zero Dirichlet boundary conditions on left and right and zero Neumann boundary conditions on the top and bottom.

Problem Definition

The following variables will define our problem:

  • c, a, f, d: The coefficients of the PDE.

c = 1;
a = 0;
f = 0;
d = 1;

% Create a PDE Model with a single dependent variable
numberOfPDE = 1;
pdem = createpde(numberOfPDE);

% Create the geometry and append it to the PDE Model
% For more information, please see the documentation page for |squarereg|
% and |pdegeom|.
g = @squareg;

Boundary Conditions

Plot the geometry and display the edge labels for use in the boundary condition definition.

pdegplot(pdem, 'edgeLabels', 'on');
axis([-1.1 1.1 -1.1 1.1]);
title 'Geometry With Edge Labels Displayed';

% Edges 1 and 3 are free
bTopBot = applyBoundaryCondition(pdem,'Edge',([1 3]), 'g', 0);
% Solution is zero on edges 2 and 4
bLeftRight = applyBoundaryCondition(pdem,'Edge',([2 4]), 'u', 0);

Generate Mesh

msh = generateMesh(pdem);
axis equal

Generate Initial Conditions

The initial conditions:

  • $u(x,0) = \arctan\left(\cos\left(\frac{\pi x}{2}\right)\right)$.

  • $\left.\frac{\partial u}{\partial t}\right|_{t = 0} = 3\sin(\pi x) \exp\left(\sin\left(\frac{\pi y}{2}\right)\right)$.

This choice avoids putting energy into the higher vibration modes and permits a reasonable time step size.

[p,~,t] = meshToPet(msh);
x = p(1,:)';
y = p(2,:)';
u0 = atan(cos(pi/2*x));
ut0 = 3*sin(pi*x).*exp(sin(pi/2*y));

Define Time-Discretization

We want the solution at 31 points in time between 0 and 5.

n = 31;
tlist = linspace(0,5,n);

Find FEM Solution

uu = hyperbolic(u0,ut0,tlist,pdem,c,a,f,d);
428 successful steps
62 failed attempts
982 function evaluations
1 partial derivatives
142 LU decompositions
981 solutions of linear systems

Animate FEM Solution

To speed up the plotting, we interpolate to a rectangular grid.

delta = -1:0.1:1;
[uxy,tn,a2,a3] = tri2grid(p,t,uu(:,1),delta,delta);
gp = [tn;a2;a3];
umax = max(max(uu));
umin = min(min(uu));
for i = 1:n
    axis([-1 1 -1 1 umin umax]);
    caxis([umin umax]);
    M(i) = getframe;