A Wave Guide

This example is an application from wave scattering theory. You send a plane wave into a wave guide that has a sharp bend the size of which is comparable to the wave-length. Due to the reflections on the walls of the wave-guide, the plane wave is scattered. Say you are interested to see how the shape of the bending affects the solution at the exit of the wave-guide. You compare 3 different shapes by computing the energy of the scattered waves at the exit.

You will see how working with the Graphical User Interface (GUI), which is specialized for PDE solving, is combined with the command-line mode, to make full use of MATLAB's general programming environment.

A wave guide

The Equation

The oscillations are modeled by the wave equation
U(x,y,t)' = div( grad(U) )
The time-dependence is eliminated if the wave is monochromatic and then
U(x,y,t) = exp(2*i*k*t*pi)*u(x,y). The wave equation is then reduced to the famous Helmholtz equation:
-4*pi^2*k^2*u = div( grad(u)),
which fits the framework of the PDE Toolbox.

The Boundary conditions

There are 3 types of boundary conditions: solid walls, entries and exits. Choose the wavelength 1/5 (and thus k = 5).

Note: the equation has complex solutions due to the boundary conditions.

The geometries

The simplest bend is an L-shaped domain. Drawn it with the GUI as a single geometrical object, a polygon. Then try a sloping side and finally a rounded corner.

Three geometries that bend a planar wave

The first two domains are polygons, and the third is defined by the set formula P1+(C1-SQ1).

Set the corresponding boundary conditions and export the information from the GUI to your MATLAB workspace. Call the geometry and boundary matrices g1, b1, g2, b2, g3, b3 denoting the first, second and third domain.

Working in command-line mode

The PDE Toolbox offers you the possibility to choose the environment that suits you best at different stages of your work. In this case you may want to split your work as follows:

Some operations are just as easy to execute both in the GUI and in the command-line mode. Here is an example of how to solve the Helmholtz equation:

You want to have 8 - 10 grid-points per wavelength If you put the triangle size of the initial mesh to 0.2 (a wavelength), then 3 uniform refinements produce the desired density of mesh-points.

[p,e,t] = initmesh(g3,'Hmax',0.2);
for j = 1:3
    [p,e,t] = refinemesh(g3,p,e,t);
The p, e, t matrices (points, edges, triangles) describe the triangular mesh. They are created by the function initmesh and by three successive calls to refinemesh The two rows of p contain the x and y coordinates of the mesh-points. The solution u is produced by calling assempde with the coefficients c = 1, a = -100*pi^2, f = 0:

u = assempde(b3,p,e,t,1,(2*pi*i*5)^2,0);
The pdeplot command produces the following plot:

Scattering in a bended wave-guide.

The vector u is a discrete approximation of the solution to the equation:
-div( grad(u) ) - 100*pi^2*u=0.
The entries of u are the values of the solution corresponding to the mesh-points given in p. Find the nodes at the exit, y == 0.8, and then the corresponding values in the solution:

index = find( p(2,:)==0.8 );              
x = p(1,index);
v = u(index);

x contains the grid-points at the exit, while v contains the complex values of the solution. Because the original mesh is unstructured, you must sort p (and v )

[x,index] = sort(x);
v = v(index);

Absolute values of the solutions at the exit.

At this stage, x and v are ordinary vectors containing grid-points and grid-values. The whole array of MATLAB facilities are available to process them. Compute the energy of the scattered waves at the exit and compare them with the energy of the ingoing wave (the ingoing energy is 0.4). The comparison is done by the following code, using the trapezoidal rule:

N = length(x);
w = abs(v).^2;
E = sum(0.5*(w(2:N)+w(1:N-1)).*diff(x))/0.4;
The L-shape bending recovers only 11.02% of the energy, compared to 80.42% and 89.53% in the plane mirror and curved mirror cases.