Numerical Solution of the Wave Equation: Finite Difference Method (FDM)
Method
Assume we are given the following one-dimensional wave equation.
$$ \dfrac{\partial^{2} u}{\partial t^{2}} = c^{2} \dfrac{\partial^{2} u}{\partial x^{2}}, \qquad 0 \le x \le 1, \quad t \ge 0 \tag{1} $$
Our goal is to approximate the above solution using a finite number of points. Let’s discretize the spacetime domain as follows.
$$ \left\{ (\ell \Delta x, n\Delta t) : \ell=0,1,\dots,d+1,\ n\ge 0 \right\}\quad \text{ where } \Delta x = \frac{1}{d+1} $$
We denote the approximation of $u(\ell \Delta x, n\Delta t)$ as $u_{\ell}^{n}$. Note that the superscript ${}^{n}$ is not a power but an index for time. The left-hand and right-hand sides of $(1)$ are approximated as follows. Reference.
$$ \dfrac{\partial^{2} u(x,t)}{\partial x^{2}} = \dfrac{1}{(\Delta x)^{2}} \big[ u(x-\Delta x, t) - 2u(x,t) + u(x+\Delta x, t) \big] + \mathcal{O}\left( (\Delta x)^{2} \right) $$
$$ \dfrac{\partial^{2} u(x,t)}{\partial t^{2}} = \dfrac{1}{(\Delta t)^{2}} \big[ u(x, t-\Delta t) - 2u(x,t) + u(x, t+\Delta t) \big] + \mathcal{O}\left( (\Delta t)^{2} \right) $$
Now, let us denote $\mu = \dfrac{\Delta t}{\Delta x}$. By substituting both expressions into $(1)$ and rearranging, we obtain the following.
$$ u(x, t-\Delta t) - 2u(x,t) + u(x, t+\Delta t) \approx \mu^{2} \big[ u(x-\Delta x, t) - 2u(x,t) + u(x+\Delta x, t) \big] $$
$$ \implies u(x, t+\Delta t) \approx 2u(x,t) - u(x, t-\Delta t) + \mu^{2} \big[ u(x-\Delta x, t) - 2u(x,t) + u(x+\Delta x, t) \big] $$
$$ \implies u(x, t+\Delta t) \approx 2(1 - \mu^{2})u(x,t) - u(x, t-\Delta t) + \mu^{2} \big[ u(x-\Delta x, t) + u(x+\Delta x, t) \big] $$
$$ \implies u_{\ell}^{n+1} = \mu^{2} u_{\ell+1}^{n} + 2(1 - \mu^{2})u_{\ell}^{n} + \mu^{2} u_{\ell-1}^{n} - u_{\ell}^{n-1},\quad \ell=1,\dots,d,\ n=0,1,\dots \tag{2} $$
Here, the ratio of the spatial interval to the temporal interval, $\mu = \dfrac{\Delta t}{\Delta x}$, is a crucial constant known as the Courant number, which determines the convergence of $(2)$.
Explanation
The above method converges when $\mu \le 1$.
Code
Setting the boundary conditions to $0$ and performing the calculations yields the following result.
using Plots
x = LinRange(-1, 1, 300)
Δx = x[2] - x[1]
t = LinRange(0, 4, 600)
Δt = t[2] - t[1]
μ = Δt / Δx
U = zeros(length(x), length(t))
U[:, 1] = exp.(-30 * x.^2)
U[:, 2] = exp.(-30 * x.^2)
for i ∈ 3:length(t)
U[2:end-1, i] = (μ^2)*U[3:end, i-1] + 2(1-μ^2)*U[2:end-1, i-1] + (μ^2)*U[1:end-2, i-1] - U[2:end-1, i-2]
end
heatmap(U)
anim = @animate for i ∈ 1:length(t)
plot(x, U[:, i], xlims=(-1, 1), ylims=(-2,2), legend=false)
end
gif(anim, fps=50)