Appendix: Numerical solution method
Approximating the wave equation
We introduce a uniform mesh in space and time with spacings
\( \Delta x \) and \( \Delta t \), respectively. At each mesh
point
$$ (x_i,t_n),\quad x_i=i\Delta x,\ i=0,\ldots,N_x,\quad t_n=n\Delta t,
\ n=0,\ldots,N_t,$$
the wave equation is approximated by second-order finite differences,
$$
\begin{align*}
\frac{\partial^2}{\partial t^2}u(x_i,t_n) &\approx
\frac{u_i^{n+1} - 2u_i^n + u^{n-1}_i}{\Delta t^2},\\
\frac{\partial^2}{\partial x^2}u(x_i,t_n) &\approx
\frac{u_{i+1}^{n} - 2u_i^n + u^{n}_{i-1}}{\Delta x^2},
\end{align*}
$$
where \( u^n_i \) is the numerical approximation to the exact solution
at \( (x_i,t_n) \).
These approximations give rise to an explicit scheme where a new
value \( u^{n+1}_i \) is given by the formula
$$
\begin{equation}
u^{n+1}_i = -u^{n-1}_i + 2u^n_i + C^2
\left(u^{n}_{i+1}-2u^{n}_{i} + u^{n}_{i-1}\right),
\tag{6}
\end{equation}
$$
where the parameter \( C \),
$$
\begin{equation}
C = c\frac{\Delta t}{\Delta x},
\end{equation}
$$
is known as the dimensionless
Courant number. A stable time-stepping
is ensured only if \( C\leq 1 \). The value \( C \) governs also the accuracy
of the method: \( C=1 \) actually reproduces the exact solution (!), while
the accuracy is reduced when decreasing \( C \).
Approximating the initial conditions
The initial condition \( u(x,0)=I(x) \), where \( I(x) \) is a prescribed
mathematical function, is implemented in the numerical method by
$$ u(x_i,0) = I(x_i),\quad i=0,\ldots,N_x\tp$$
The other initial condition,
$$ \frac{\partial}{\partial t}(x,0) = 0,\hbox{ or }
\frac{\partial}{\partial t}(x,0) = V(x),$$
is approximated by a centered difference over an interval \( 2\Delta t \).
When this difference is combined with
(6),
we get a special formula for \( u^1_i \). Thereafter, for \( n\geq 1 \),
one can apply
(6).
Approximation of boundary conditions
The finite difference scheme
(6) is
applied at all
inner points in the spatial mesh, \( i=1,\ldots,N_x-1 \).
For \( i=1 \) or \( i=N_x-1 \),
(6) involves
the boundary points \( u^n_0 \) or \( u^n_{N_x} \) (respectively) at the
previous time step, but these are condidered known.
Dirichlet conditions
In case of \( u=0 \) (Dirichlet) conditions, we just set \( u^{n+1}_0=0 \)
and \( u^{n+1}_{N_x}=0 \). Feeding a wave in from the left is just
a matter of setting \( u^{n+1}_0 \) equal to the known boundary function of
time: \( u^{n+1}_0=U_0((n+1)\Delta t) \).
Neumann conditions
No-flux or Neumann conditions are
discretized by a centered finite difference and combined with
(6), yielding a modified form of
(6) at the boundary:
$$
\begin{align}
u^{n+1}_i &= -u^{n-1}_i + 2u^n_i + 2C^2
\left(u^{n}_{i+1}-u^{n}_{i}\right),\quad i=0\\
u^{n+1}_i &= -u^{n-1}_i + 2u^n_i + 2C^2
\left(u^{n}_{i-1}-u^{n}_{i}\right),\quad i=N_x\tp
\end{align}
$$
Open boundary conditions
Radiation, artificial, or open boundary conditions of the type
(3)-
(4) are
discretized by upstream first-order, forward in time differences,
resulting in explicit updating formulas:
$$
\begin{align}
u^{n+1}_i &= u^n_i + C(u^n_{i+1} - u^n_i),\quad i=0,\\
u^{n+1}_i &= u^n_i - C(u^n_{i} - u^n_{i-1}),\quad i=N_x\tp
\end{align}
$$
Even though the underlying finite differences are of first order only,
these conditions let the waves out of the domain perfectly and do not
lower the accuracy of the finite difference scheme used in the interior
of the domain or at no-flux (Neumann) boundaries.
Periodic conditions
The periodic condition \( u(0,t)=u(L,t) \) is trivial to incorporate in
the numerical method:
$$ u^{n+1}_0 = u^{n+1}_{N_x}\tp$$