$$ \newcommand{\tp}{\thinspace .} $$

 

 

 

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$$