Partial Differential Equation Toolbox

Helmholtz's Equation on a Unit Disk with a Square Hole

This example shows how to solve a Helmholtz equation using the assempde function in the Partial Differential Toolbox™.

The Helmholtz equation, an elliptic equation that is the time-independent form of the wave equation, is

$-\Delta u-k^2u = 0$.

Solving this equation allows us to compute the waves reflected by a square object illuminated by incident waves that are coming from the right.

Problem Definition

The following variables will define our problem:

  • g: A specification function that is used by initmesh. For more information, please see the documentation pages for scatterg and pdegeom.

  • k, c, a, f: The coefficients and inhomogeneous term.

g = @scatterg;
k = 60;
c = 1;
a = -k^2;
f = 0;

Boundary Conditions

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

pdegplot(g, 'edgeLabels', 'on');
axis equal
title 'Geometry With Edge Labels Displayed';
% Create a pde entity for a PDE with a single dependent variable
numberOfPDE = 1;
pb = pde(numberOfPDE);
% Create a geometry entity
pg = pdeGeometryFromEdges(g);
bOuter = pdeBoundaryConditions(pg.Edges(5:8), 'g', 0, 'q', -60i);
innerBCFunc = @(thePde, loc, state) -exp(-1i*k*loc.x);
bInner = pdeBoundaryConditions(pg.Edges(1:4), 'u', innerBCFunc);
pb.BoundaryConditions = [bOuter bInner];

Create Mesh

We need a fine mesh to resolve the waves. To achieve this, we refine the initial mesh twice.

[p,e,t] = initmesh(g);
[p,e,t] = refinemesh(g,p,e,t);
[p,e,t] = refinemesh(g,p,e,t);
axis equal

Solve for Complex Amplitude

The real part of the vector u stores an approximation to a real-valued solution of the Helmholtz equation.

u = assempde(pb,p,e,t,c,a,f);

Plot FEM Solution


Animate Solution to Wave Equation

Using the solution to the Helmholtz equation, we construct an animation showing the corresponding solution to the time-dependent wave equation.

m = 10;
h = newplot;
hf = h.Parent;
axis tight
ax = gca;
ax.DataAspectRatio = [1 1 1];
axis off
maxu = max(abs(u));
for j = 1:m
    uu = real(exp(-j*2*pi/m*sqrt(-1))*u);
    caxis([-maxu maxu]);
    axis tight
    ax = gca;
    ax.DataAspectRatio = [1 1 1];
    axis off
    M(j) = getframe(hf);