📂Numerical Analysis

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)