Source: wave-equation.pdf

To describe the implementation of the wave equation in Octave we first have to derive an algorithm to numerically solve the differential equation. This is the general one-dimensional case of the wave equation:

$$ \frac{ \partial^2 u(x, t) }{ \partial x^2 } = \frac{1}{c^2} \frac{ \partial^2 u(x, t)}{\partial t^2}$$

Where $u(x,t)$ stands for the deflection at the point $x$ and the time $t$ and $c$ is the wave propagation speed in the medium. The second derivation can be approximated by a central difference:

$$ \frac{d^2 f(x)}{d x^2}\approx \frac{f(x+\Delta x) - 2f(x) + f(x-\Delta x)}{\Delta x^2}$$

which can be substituted in the first equation:

$$ \frac{u(x+\Delta x,t) - 2u(x,t) + u(x-\Delta x,t)}{\Delta x^{2}} = \frac{1}{c^{2}} \frac{u(x,t+\Delta t) - 2u(x,t) + u(x,t-\Delta t)}{\Delta t^{2}}$$

To maintain a clear view we will rename some functions as follows:

$$ c(x) = u(x,t) p(x) = u(x, t-\Delta t) f(x) = u(x, t+\Delta t)$$

So the equation will simplify to:

$$ \frac{c(x+\Delta x) - 2c(x) + c(x-\Delta x)}{\Delta x^{2}} = \frac{1}{c^{2}} \frac{f(x) - 2c(x) + p(x)}{\Delta t^{2}}$$

Now we have to transform the equation to express f(x). The reason for this step is to get a recursive form, which can be easily implemented into a program.

$$ step 1: \\ c^{2}\Delta t^{2} \frac{c(x-\Delta x) - 2c(x) + c(x+\Delta x)}{\Delta x^{2}} = f(x) - 2c(x) + p(x) \\ step 2: \\ f(x) = 2c(x) - p(x) + r^{2} \lbrack c(x-\Delta x) - 2c(x) + c(x+\Delta x) \rbrack$$

Where $r=\frac{c \Delta t}{\Delta x}$. This is the final equation which can now be used as an algorithm for a simulation.