Finite Element Analysis Toolbox
Interference and Diffraction

This classic double-slit experiment considers a planar and periodic oscillating wave, which hits and passes two narrow slits. Assuming the slits are narrow enough, the passing waves will bend and cause an interference pattern, while diffraction will attenuate the resulting off-axis amplitude.

In the example the Helmholtz equation is used to model the wave phenomena

\[ -( \frac{\partial^2 A}{\partial x^2} + \frac{\partial^2 A}{\partial y^2} ) - k^2 A = 0 \]

where A is the amplitude of the wave,and k the wave number (k = 2·π/λ).

Tutorial

This model is available as an automated tutorial by selecting Model Examples and Tutorials... > Classic PDE > Interference and Diffraction from the File menu, viewed as a video tutorial, or alternatively, follow the step-by-step instructions below.

  1. To start a new model click the New Model toolbar button, or select New Model... from the File menu.
  2. Select the Custom Equation physics mode from the Select Physics drop-down menu.
  3. Enter A into the Dependent Variable Names edit field. This will represent the amplitude of the wave.

  4. Press OK to finish the physics mode selection.

The geometry consists of two small rectangles, representing the slits, connected to a half circle domain. First create a circle centerd at the origin with radius 0.8.

  1. Select Circle from the Geometry menu.
  2. Enter 0.8 into the radius edit field.
  3. Press OK to finish and close the dialog box.

Then create and subtract a rectangle from the lower half of the circle.

  1. Select Rectangle from the Geometry menu.
  2. Enter -0.8 into the xmin edit field.
  3. Enter 0.8 into the xmax edit field.
  4. Enter -0.8 into the ymin edit field.
  5. Enter 0 into the ymax edit field.
  6. Press OK to finish and close the dialog box.
  7. Select C1 and R1 in the geometry object Selection list box.
  8. Press the - / Subtract geometry objects Toolbar button.

  9. Select Rectangle from the Geometry menu.
  10. Enter -0.08-0.02 into the xmin edit field.
  11. Enter -0.08+0.02 into the xmax edit field.
  12. Enter -0.2 into the ymin edit field.
  13. Enter 0 into the ymax edit field.
  14. Press OK to finish and close the dialog box.
  15. Select Rectangle from the Geometry menu.
  16. Enter 0.08-0.02 into the xmin edit field.
  17. Enter 0.08+0.02 into the xmax edit field.
  18. Enter -0.2 into the ymin edit field.
  19. Enter 0 into the ymax edit field.
  20. Press OK to finish and close the dialog box.
  21. Select CS1, R2, and R3 in the geometry object Selection list box.
  22. Press the + / Add geometry objects Toolbar button.

  23. Switch to Grid mode by clicking on the corresponding Mode Toolbar button.
  24. Enter 0.015 into the Grid Size edit field, and press the Generate button to call the grid generation algorithm.

  25. Switch to Equation mode by clicking on the corresponding Mode Toolbar button.
  26. Press the edit button.
  27. Enter -(Ax_x + Ay_y) - k^2*A_t = 0 into the Equation for A edit field.

  28. Press OK to finish and close the dialog box.
  29. Press OK to finish the equation and subdomain settings specification.
  30. Press the Constants Toolbar button, or select the corresponding entry from the Equation menu, and add the following modeling constants for the wave length wl, and wave number k in the Model Constants and Expressions dialog box.
Name Expression
wl 0.08
k pi*2/wl


  1. Switch to Boundary mode by clicking on the corresponding Mode Toolbar button.

First set homogeneous Neumann conditions for all the boundaries.

  1. Select all boundaries (1 - 11) in the Boundaries list box.
  2. Then select the Neumann, g_A radio button, and enter 0 into the Dirichlet/Neumann coefficient edit field.

An incoming planar wave is present at the inlet, with the complex boundary condition n·∇(A) + k·i·A = 2·k·i, which can be implemented as a Neumann flux type boundary condition.

  1. Select 4 and 8 in the Boundaries list box, and enter -k*i*A + 2*k*i into the Dirichlet/Neumann coefficient edit field.

The outlet is assumed non-reflective and therefore n·∇(A) + k·i·A = 0.

  1. Select 1 and 2 in the Boundaries list box.
  2. Enter -k*i*A into the Dirichlet/Neumann coefficient edit field.

  3. Press OK to finish the boundary condition specification.
  4. Switch to Solve mode by clicking on the corresponding Mode Toolbar button.
  5. Press the = Toolbar button to call the solver. After the problem has been solved FEATool will automatically switch to postprocessing mode and plot the computed solution.

After the problem has been solved FEATool will automatically switch to postprocessing mode, and display the computed wave amplitude A. The interference pattern can clearly be seen, as the four distinct lines where the waves have been completely canceled out.

The Point/Line Evaluation tool can be used to visualize the interference and diffraction pattern at the boundary.

  1. Select Point/Line Evaluation... from the Post menu.
  2. Select 10 and 11 in the Boundaries list box. This will automatically enter the corresponding boundary coordinates into the Evaluation Coordinates edit fields. Then press Apply or OK to finish and plot the amplitude curve.

Even though the values for boundary 11 are plotted in reverse, one can clearly see the four intersections with the zero amplitude line. From theory the minima and maxima are given as sin(θ) = (n+½)λ/L and nλ/L, respectively (for n=0, ±1, etc where L is the distance between the slits), which here should be 0 and 30 degrees for the maxima, consistent with the results of the simulation.

The interference and diffraction custom equation model has now been completed and can be saved as a binary (.fea) model file, or exported as a programmable MATLAB m-script text file, or GUI script (.fes) file.

Video