Simulation of Wave Motion


Ui,j Wave height at point i,j
i Grid index in x-direction
j Grid index in y-direction
n Iteration number
dt Timestep
mu Viscosity coefficient
il Number of grid points in x-direction
jl Number of grid points in y-direction



Wave behavior has always been a large area of study in science. Waves are a part of our everyday lives. You can't go through a day without experiencing a wave of one sort or another. There are waves in water, waves in the air, even waves in solids. There are different sizes and shapes of waves. Waves move at different speeds. There are different wave frequencies. You create waves every time you make a sound. Radio stations transmit signals in the form of radio waves. Depending on the surroundings, waves can behave in different ways. Sound waves are absorbed in a soundproof room yet they bounce of the sides of a metal drum creating a different sound. Over long distances, sound waves will eventually die out. Why does this happen? How can we study wave behavior using computers?

Fortunately, some simple forms of waves can be described mathematically. Waves in real life are all three-dimensional in nature. However, to simplify the problem, we will consider the behavior of a wave in two dimensions. A two dimensional wave can be likened somewhat to a surface wave on water. Most everyone has seen the ripple effect resulting from a stone being dropped in a puddle. There is even a name for the equation that describes wave motion in two dimensions. The "2D Wave Equation" is given below:

where "a" is the speed at which the wave travels (speed of sound if we are talking about sound waves).

This equation is what is known as a partial differential equation. The terms in the equation are called second partial derivatives. Derivatives represent how one variable changes with respect to a change in another variable. Velocity, for example, is the first derivative of position with respect to time. Likewise, acceleration is the first derivative of velocity with respect to time or the second derivative of position with respect to time.

To get into too much more detail about the theory behind this equation would be beyond the scope of this problem. Suffice it to say that the 2-D wave equation describes the motion of a wave by relating the time acceleration of the wave to its spatial gradients (derivatives in x and y).

The equation alone, however, is not enough to fully describe the behavior of a wave. What is the initial shape of the wave? What kind of boundaries will the wave run into? How will the initial shape and boundaries affect the wave? We will need to consider initial conditions and boundary conditions when we try to answer these questions.


How do we use an equation such as the 2-D wave equation on a computer? The first step is to "discretize" the equation. That is, approximate the equation at a single point in space so that it can be solved at that single point using values at surrounding points. To do this, we first need to set up a grid. In two dimensions this grid would look like this:

A discretized equation is calculated only at the intersections of the grid lines. Each one of the partial derivatives can be discretized using what are called "central differences." If we blow up the highlighted section of the grid we can see how the equation is discretized.

With the 2-D space discretized in this fashion, we can approximate the derivatives as follows:

If we substitute these approximations into the original equation we get



As was mentioned before, the behavior of a wave depends upon its initial shape. The initial shape can be likened to an initial disturbance. For this problem we will give an initial disturbance of a sine wave (actually a cosine wave) at the very center of the grid. The calculation will then track the propagation of the wave outward toward the boundaries. The wave should bounce or reflect off the boundaries and return toward the center.

Let's assume that the grid length is 1.0 in each direction, and that it is centered about the point (0, 0). This means that the left and right edges are at x=-0.5 and x=0.5 respectively. Likewise, the bottom and top edges are at y=-0.5 and y=0.5 respectively. If we want the wave to be symmetric about the origin, then we must use polar coordinates. The equation for our initial disturbance is:


This will give an initial distribution that looks like this

If we look again at the discretized equation, we see that to calculate the wave height at iteration 2 we need values at iterations 0 and 1. We know the values at n=1 as the initial conditions (see previous figure). What are the values at n=0? We apply a zero-gradient initial condition. Mathematically, this means that

In discrete form this means

If we substitute this into the original discretized equation we get

for the very first iteration. This is known as the "starting formula."


Another factor in the behavior of a wave is the presence of different types of boundaries. For this problem we will consider two types of boundary conditions. These two types of boundaries are fixed and free boundaries. To understand these two types of boundaries, picture a string that is attached to a wall at one end and free to move at the other end. If you pick up the string in the middle and shake it, the wave will move in both directions along the string. The wave will behave differently at the fixed end than it will at the free end. How do we describe these conditions mathematically?

Suppose the whole left edge of our 2-D grid is fixed. The discrete form of the boundary condition is

Likewise, suppose the whole right edge of the grid is free. The discrete form of the boundary condition is

Substituting this into the main discretized equation we get

We can apply similar conditions to the top and bottom edges of the grid.


Following is a general outline of the code to solve the wave equation:

  • Declare variables and set constants
  • Calculate x, y arrays.
  • Calculate initial values for U.
  • Use starting formula, boundary conditions for first iteration.
  • Apply general formula, boundary conditions until n = NEND, writing out a frame of data every NFREQ iterations
  • Done


Last iteration (somewhere around 300 - 500)


Number of iterations to wait between writing output (~10)


What good are numbers if we can't understand what they really mean? The old saying, "A picture is worth a thousand words," truly applies here. Our simulation will calculate the wave heights at 10000 grid points for several hundred time iterations. If we consider that each number describing a wave height is one word long, then one frame of the simulation (the solution frozen at an instant in time) contains 10000 words. To understand wave motion we need to animate all of the calculated frames. Assuming we have 180 frames, the total number of words is 1800000. So, in our case, a picture is worth 1.8 million words.

To visualize the wave, you will use two methods. The first and most illustrative method is to displace the 2-D grid in the 3rd dimension based on the value of the wave height. The second method is to color the grid points based on the value of the wave height.

You will be using AVS/Express to visualize the results. Since the grid is uniform, you only need to write out the values of wave height for each grid point.


  • Use different combinations of boundary conditions and initial conditions to further investigate wave behavior.
  • Use a non-uniform grid. This may require you to change how you do your visualizations.
  • Add the following term to the left side of the 2-D wave equation and vary mu to determine its effect on the wave:
  • Instead of the initial conditions shown above, initialize all values to zero and add a forcing function f(t) at some location in the grid. This can be either calculated (eg. f(t)=cos(2*pi*omega*t)) or read from a file (eg. data from a digital recording of a sound).






SI 1996

SI 1995

Troy Baer is the coordinator for the Wave Motion project. His office is in OSC, cubicle 420-6, phone 292-9701.

For assistance, write or call 614-292-0890.