In the first part we considered a system composed of a chain of connected springs and masses, and came to the conclusion that the vertical position of the masses is given by the solution of a system of interdependent equations.

This system of equations is of the following form, where *a**i* is the acceleration and *U**i* the position of the mass at position *i*:

*m a*i* = k ( U**i-1** – 2 U**i** + U**i+1** )*

For each step of the computation, the new state of the system must be computed using its previous states. This is a basic mechanism for any kind of numeric simulation. The remaining difficulty is the correct representation of the acceleration term, which can be transformed using a first-order approximation, e.g. the derivative of *U* with respect to *t* can be given by:

*∂U(t)/∂t ≈ ( U(t+Δt) – U(t) ) / Δt*

We apply this formula two times to get the second order derivative (acceleration):

*∂*²*U(t)/∂t*²* ≈ ( U(t+2Δt) – U(t+Δt) – ( U(t+Δt) – U(t) ) ) / Δt / Δt*

*∂*²*U(t)/∂t*²* ≈ ( U(t+2Δt) – 2 U(t+Δt) + U(t) ) / Δt*²

Another way to write this formula:

*∂*²*U(t)/∂t*²* ≈ ( U(t+**Δt**) – 2 U(t) + U(t-Δt) ) / Δt*²

At this point, it is important to understand how to achieve the practical computation. The expected result is a set of values (the *U**i*) evolving over the time. The time between two consecutive steps will be exactly *Δt.* But the approximation of the acceleration shows that the computation of a new state of the *U**i* depends not only on the previous state, but also on the state before the previous state (i.e. 2 states in the past due to the *t-2Δt* value). Practically, this means that a total of 3 states of the *U**i* must be kept in memory: the current state being computed plus two steps in the past.

The usage of different steps of the computation shows that we now need to clearly identify the computed values using their temporal and spatial coordinates. An obvious simplification is to set:

*Δt = 1*

After that, let’s use the notation *U**i**(t)* to identify a computed value at position *i* and step *t*. We put the approximation in place of the acceleration. The resulting formula to compute the *U**i**(t+1)* is then:

*m ( U**i**(t+1) – 2 U**i**(t) + U**i**(t-1) ) = k ( U**i-1**(t)** **– 2 U**i**(t)** + U**i+1**(t)** )*

With this formula, the computation is straightforward because there is only one term at *t+1* and all other terms are given by the previous steps. It is astonishing to see the relative simplicity of the formula (a linear combination of the previous values, no special case to handle), keeping in mind that it mimics the behavior of a physical system.

This result shows a kind of symmetry between the parts on left and right: the same pattern *U(n+1)-2U(n)+U(n-1)* is present, with varying time index on the left hand side and varying spatial index on the right hand side. The part on the left is the approximation of a second order derivative with respect to *t*. So it appears natural to try to interpret the part on the right similarly.

This can be done and results in a well-known equation that can model any linear propagation phenomenon, i.e. where each propagated wave is independent of the others. This is the Wave Equation, given here in 1 dimension:

*∂**²U/**∂**t² = c² **∂**²U/**∂**x²*

To be continued…