Documentation

## 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

$\frac{{\partial }^{2}u}{\partial {t}^{2}}-\nabla \cdot \nabla u=0.$

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

$m\frac{{\partial }^{2}u}{\partial {t}^{2}}-\nabla \cdot \left(c\nabla u\right)+au=f.$

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);
geometryFromEdges(model,@squareg);
pdegplot(model,'EdgeLabels','on');
ylim([-1.1 1.1]);
axis equal
title 'Geometry With Edge Labels Displayed';
xlabel x
ylabel y

Specify PDE coefficients.

specifyCoefficients(model,'m',m,'d',0,'c',c,'a',a,'f',f);

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,'dirichlet','Edge',[2,4],'u',0);
applyBoundaryCondition(model,'neumann','Edge',([1 3]),'g',0);

Create and view a finite element mesh for the problem.

generateMesh(model);
figure
pdemesh(model);
ylim([-1.1 1.1]);
axis equal
xlabel x
ylabel y

Set the following initial conditions:

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

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

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

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);
457 successful steps
38 failed attempts
992 function evaluations
1 partial derivatives
112 LU decompositions
991 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.

figure
umax = max(max(u));
umin = min(min(u));
for i = 1:n
pdeplot(model,'XYData',u(:,i),'ZData',u(:,i),'ZStyle','continuous',...
'Mesh','off','XYGrid','on','ColorBar','off');
axis([-1 1 -1 1 umin umax]);
caxis([umin umax]);
xlabel x
ylabel y
zlabel u
M(i) = getframe;
end

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

Get trial now