Last time we came to a formula that can be used to compute the simulation of a one-dimensional masses-springs system, and saw that it could be interpreted as approximation of the wave equation. This introduction was intended to show that a simple model, used primarily for a numeric simulation, leads to the same properties as a continuous system. We will use the wave equation as basis in the next posts and no more discrete system, because it yields just the same results.

The general form of the wave equation in 3 dimensions is:

*∂²U/∂t² = c² ∆U = c² ( ∂²U/∂x² + ∂²U/∂y² + ∂²U/∂z² )*

Where *c* is the speed of the wave (i.e. "celerity"), *∆ *the Laplacian operator and U a scalar function of *x*, *y*, *z* and *t*. In order to make the simulation more interesting while keeping it computable, let’s look at approximated formulas in two dimensions. We are looking for the computation of *U(t,x,y)*, and assume that we can approximate the second-order derivatives by discretization using *∆x = ∆y = ∆t = 1*:

*∂²U/∂t² = c² ( ∂²U/∂x² + ∂²U/∂y² )*

*U(t+1,x,y) – 2 U(t,x,y) + U(t-1,x,y) = c² ( U(t,x+1,y) – 2 U(t,x,y) + U(t,x-1,y) + U(t,x,y+1) – 2 U(t,x,y) + U(t,x,y-1) )*

*U(t+1,x,y) = c² ( U(t,x+1,y) + U(t,x-1,y) + U(t,x,y+1) + U(t,x,y-1) ) + (2-4c²) U(t,x,y) – U(t-1,x,y)*

An obvious simplification of the formula can be reached by assuming *c = 1*, but as a consequence the propagation speed is fixed hence a phenomenon involving an interface between different mediums (e.g. refraction) cannot be simulated. Furthermore, the propagation speed would be exactly the same as the discretization intervals (e.g. *∆t*). The result is that the waves move of exactly one table cell per step, what may lead to unwanted effects. With this assumption, the formula becomes:

*U(t+1,x,y) = U(t,x+1,y) + U(t,x-1,y) + U(t,x,y+1) + U(t,x,y-1) – 2 U(t,x,y) – U(t-1,x,y)*

Graphically, the computation of this formula can be represented as following:

Now, to define the computation completely, it is necessary to clarify the behavior at the borders and the initial conditions. For the borders, the simplest approach is to keep a border of constant cells (not processed in the computation) set to 0. The resulting behavior on the borders is a reflection of all incoming waves. Because this behavior may not be wished for a particular simulation activity, algorithms have been developed to "absorb" the incoming waves on the borders, but I won’t address them in this post. Another method, simpler but increasing the computation time, is just to extend the simulated area to a such a size that reflected waves won’t reach the interesting area before the end of the simulated time.

The initial condition is another interesting topic because 2 steps are necessary. Basically, the 2 initial steps can be set to 0 for a stable start (nothing moves unless a perturbation is introduced). Any perturbation can be inserted in the computation in either one or two consecutive steps. If only one step is changed, the perturbation has no speed and will propagate equally in all directions. If two computation steps are modified, it is possible to give a direction to the introduced perturbation, under the condition that the celerity of the medium is respected.

That’s all for now regarding wave propagation. There is not a lot of information in these 3 posts, but enough to give primary ideas and to start to develop a simple simulation. Many other methods exist for calculations based on finite differences. The goal here was to just look at one that can be understood with basic mathematics. I will come back to this topic sometimes in the (maybe far) future.