Main Content

Wave Equation on Square Domain

This example shows how to solve the wave equation using the solvepde function.

The standard second-order wave equation is


To express this in toolbox form, note that the solvepde function solves problems of the form


So the standard wave equation has coefficients m=1, c=1, a=0, and f=0.

c = 1;
a = 0;
f = 0;
m = 1;

Solve the problem on a square domain. The squareg function describes this geometry. Create a model object and include the geometry. Plot the geometry and view the edge labels.

numberOfPDE = 1;
model = createpde(numberOfPDE);
ylim([-1.1 1.1]);
axis equal
title("Geometry With Edge Labels Displayed")

Figure contains an axes object. The axes object with title Geometry With Edge Labels Displayed, xlabel x, ylabel y contains 5 objects of type line, text.

Specify PDE coefficients.


Set zero Dirichlet boundary conditions on the left (edge 4) and right (edge 2) and zero Neumann boundary conditions on the top (edge 1) and bottom (edge 3).

applyBoundaryCondition(model,"neumann","Edge",([1 3]),"g",0);

Create and view a finite element mesh for the problem.

ylim([-1.1 1.1]);
axis equal
xlabel x
ylabel y

Figure contains an axes object. The axes object with xlabel x, ylabel y contains 2 objects of type line.

Set the following initial conditions:

  • u(x,0)=arctan(cos(πx2)).

  • ut|t=0=3sin(πx)exp(sin(πy2)).

u0 = @(location) atan(cos(pi/2*location.x));
ut0 = @(location) 3*sin(pi*location.x).*exp(sin(pi/2*location.y));

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

Specify the solution times as 31 equally-spaced points in time from 0 to 5.

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

Set the SolverOptions.ReportStatistics of model to 'on'.

model.SolverOptions.ReportStatistics ='on';
result = solvepde(model,tlist);
441 successful steps
34 failed attempts
952 function evaluations
1 partial derivatives
115 LU decompositions
951 solutions of linear systems
u = result.NodalSolution;

Create an animation to visualize the solution for all time steps. Keep a fixed vertical scale by first calculating the maximum and minimum values of u over all times, and scale all plots to use those z-axis limits.

umax = max(max(u));
umin = min(min(u));
for i = 1:n
    pdeplot(model,"XYData",u(:,i),"ZData",u(:,i), ...
    axis([-1 1 -1 1 umin umax]); 
    caxis([umin umax]);
    xlabel x
    ylabel y
    zlabel u
    M(i) = getframe;

Figure contains an axes object. The axes object with xlabel x, ylabel y contains an object of type patch.

To play the animation, use the movie(M) command.