12  Wave Equation

The same method of lines can be extended to obtain a numerical solution to the wave equation, with some caveats.

Consider a string of \(L\), the 1-dimensional wave equation that describes the movement of waves in the string is given by \[\frac{\partial^{2} u}{\partial t^{2}}=c^2 \frac{\partial^{2} u}{\partial x^{2}} \quad \text{with} \quad 0<x<L \quad \text{and} \quad t>0\] where \(u=u(x,t)\) is the vertical displacement of the string at location \(x\) at time \(t\) and \(c\) is the wave speed.

This partial differential equation has four derivatives in total, two derivatives in \(x\) and two derivative in \(t\), this means that four conditions are needed, two on \(x\) and two on \(t\):

12.1 The Method of Lines for the Wave Equation

The outline of the method of lines for the wave equation is as follows:

  1. Divide the spatial interval \([0,L]\) into \(N_x\) equally sized sections and label the points as \(x_0, x_1, x_2, \dots, x_{N_x}\) where \(x_n=nh_x\) and the spatial interval width is \(h_x=\frac{L}{N_x}\).

  2. Left Hand Side: For each point \(x_n\), define the approximation \(U_n(t) \approx u(x_n,t)\). Therefore the left hand side of the wave equation can be written as \[\frac{\partial^{2} u}{\partial t^{2}}(x_n,t) \approx \frac{\mathrm{d}^{2} U_n}{\mathrm{d} t^{2}}(t)\] and this holds for \(n=1, 2, \dots, N_x-1\) since \(U_0(t) \approx u(0,t)=u_l(t)\) and \(U_{N_x}(t) \approx u(L,t)=u_r(t)\) are already known from the boundary conditions. Notice that the derivative of \(U_n\) is an ordinary derivative since \(U_n\) is a function of \(t\) only.

  3. Right Hand Side: Use the finite difference approximation to approximate the spatial derivative in the differential equation. Here, the centred difference approximation for the second derivative will be used, namely \[\frac{\partial^{2} u}{\partial x^{2}}(x_n,t) \approx \frac{U_{n+1}(t)-2U_{n}(t)+U_{n-1}(t)}{h_x^2}.\] Therefore the right hand side of the wave equation will become \[c^2 \frac{\partial^{2} u}{\partial x^{2}}(x_n,t) \approx \frac{\alpha}{h_x^2}\left[ U_{n-1}(t)-2U_n(t)+U_{n+1}(t) \right].\] This holds for \(n=1, 2, \dots, N_x-1\) bearing in mind, once again, that \(U_0(t) \approx u(0,t)=u_l(t)\) and \(U_{N_x}(t) \approx u(L,t)=u_r(t)\) are known beforehand.

  4. These can be combined to give the discretised form of the heat equation \[\frac{\mathrm{d}^{2} U_n}{\mathrm{d} t^{2}}=\frac{c^2}{h_x^2}\left[ U_{n-1}-2U_n+U_{n+1} \right]\] for all \(n=1,2,\dots,N_x-1\) where \(U_n=U_n(t)\). This means that the partial differential equation has been split into \(N_x-1\) ordinary differential equations.

  5. This system has a second order derivative in \(t\). To reduce it to a single derivative in \(t\), let \[\frac{\mathrm{d} U_n}{\mathrm{d} t}=V_n,\] this means that \[\frac{\mathrm{d} V_n}{\mathrm{d} t}=\hat{c}\left[ U_{n-1}-2U_n+U_{n+1} \right] \quad \text{where} \quad \hat{c}=\frac{c^2}{h_x^2}.\] This substitution can convert the the wave equation into \(2(N_x-1)\) equations that can now be written in matrix form as \(\frac{\mathrm{d} \boldsymbol{U}}{\mathrm{d} t}=A\boldsymbol{U}+\boldsymbol{b}\) where \[\boldsymbol{U}=\begin{pmatrix} V_1 \\ V_2 \\ \vdots \\ V_{N_x-2} \\ V_{N_x-1} \\ U_1 \\ U_2 \\ \vdots \\ U_{N_x-2} \\ U_{N_x-1} \end{pmatrix} \quad \text{and} \quad \boldsymbol{b}(t)=\begin{pmatrix} \hat{c}u_l(t) \\ 0 \\ \vdots \\ 0 \\ \hat{c} u_r(t) \\ 0 \\ 0 \\ \vdots \\ 0 \\ 0 \end{pmatrix}\] and \(A\) is a block matrix of the form \[A=\begin{bmatrix} \mathcal{O} & \hat{c}A_1 \\ \mathcal{I} & \mathcal{O} \end{bmatrix}\] where \(\mathcal{O}\) is the zero matrix of size \(N_x-1\), \(\mathcal{I}\) is the identity matrix of size \(N_x-1\) and \(A_1\) is the same matrix as in the heat equation of size \(N_x-1\) (namely the matrix with \(-2\) on the diagonal and \(1\) on the super- and sub diagonals). This differential equation can be solved subject to the initial condition \[ \boldsymbol{U}_0=\begin{pmatrix} V_1(0) \\ V_2(0) \\ \vdots \\ V_{N_x-2}(0) \\ V_{N_x-1}(0) \\ U_1(0) \\ U_2(0) \\ \vdots \\ U_{N_x-2}(0) \\ U_{N_x-1}(0) \end{pmatrix}\approx\begin{pmatrix} \frac{\partial u}{\partial t}(x_1,0) \\ \frac{\partial u}{\partial t}(x_2,0) \\ \vdots \\ \frac{\partial u}{\partial t}(x_{N_x-2},0) \\ \frac{\partial u}{\partial t}(x_{N_x-1},0) \\ u(x_1,0) \\ u(x_2,0) \\ \vdots \\ u(x_{N_x-2},0) \\ u(x_{N_x-1},0) \end{pmatrix} =\begin{pmatrix} v_{init}(x_1) \\ v_{init}(x_2) \\ \vdots \\ v_{init}(x_{N_x-2}) \\ v_{init}(x_{N_x-1}) \\ u_{init}(x_1) \\ u_{init}(x_2) \\ \vdots \\ u_{init}(x_{N_x-2}) \\ u_{init}(x_{N_x-1}) \\ \end{pmatrix}. \] This system can now be solved using any of the IVP solvers with a temporal stepsize \(h_t\).

12.2 Stability

All the diagonal entries of the matrix \(A\) are zero, this means that using the Gershgorin Circle Theorem will produce circles all centred at the origin which would not guarantee that all the eigenvalues have negative real parts. Another approach is needed.

All the eigenvalues \(\lambda\) of the matrix \(A\) must satisfy \[\mathrm{det}\left( A-\lambda \tilde{\mathcal{I}} \right)=0\] where \(\tilde{\mathcal{I}}\) is the identity matrix of size \(2(N_x-1) \times 2(N_x-1)\). The matrices \(A\) and \(\tilde{\mathcal{I}}\) can be written in block form as \[A=\begin{bmatrix} \mathcal{O} & \hat{c}A_1 \\ \mathcal{I} & \mathcal{O} \end{bmatrix} \quad \text{and} \quad \tilde{\mathcal{I}}=\begin{bmatrix} \mathcal{I} & \mathcal{O} \\ \mathcal{O} & \mathcal{I} \end{bmatrix}\] where \(\mathcal{O}\) is the zero matrix, \(\mathcal{I}\) is the identity matrix and \(A_1\) is the same as the heat matrix, all of which have size \((N_x-1) \times (N_x-1)\). This means that the eigenvalues must satisfy \[\mathrm{det}\left( A-\lambda \tilde{\mathcal{I}} \right)=0 \quad \implies \quad\mathrm{det} \begin{bmatrix} -\lambda \mathcal{I} & \hat{c}A_1 \\ \mathcal{I} & -\lambda\mathcal{I} \end{bmatrix}=0.\]

Since all the submatrices are commutative, then \[\mathrm{det} \begin{bmatrix} -\lambda \mathcal{I} & \hat{c}A_1 \\ \mathcal{I} & -\lambda\mathcal{I} \end{bmatrix}=\mathrm{det}\left( \lambda^2 \mathcal{I}-\hat{c}A \right).\]

Detereminant of Block Matrices

Consider the block matrix \(M\) given by \[M=\begin{bmatrix} A & B \\ C & D \end{bmatrix}\] where \(A\) and \(D\) are square submatrices. There are multiple formulas for obtaining the determinant of such a matrix:

  • If the matrix \(A\) is invertible, then \[\mathrm{det}\begin{bmatrix} A & B \\ C & D \end{bmatrix}=\mathrm{det}\left( A \right)\mathrm{det}\left( D-CA^{-1}B \right)\]

  • If the matrix \(D\) is invertible, then \[\mathrm{det}\begin{bmatrix} A & B \\ C & D \end{bmatrix}=\mathrm{det}\left( D \right)\mathrm{det}\left( A-BD^{-1}C \right)\]

  • If all the submatrices are square and of the same size:

    • if \(CD=DC\) \[\mathrm{det}\begin{bmatrix} A & B \\ C & D \end{bmatrix}=\mathrm{det}\left( AD-BC \right)\]
    • if \(CA=AC\) \[\mathrm{det}\begin{bmatrix} A & B \\ C & D \end{bmatrix}=\mathrm{det}\left( AD-CB \right)\]
    • if \(BD=DB\) \[\mathrm{det}\begin{bmatrix} A & B \\ C & D \end{bmatrix}=\mathrm{det}\left( DA-BC \right)\]
    • if \(AB=BA\) \[\mathrm{det}\begin{bmatrix} A & B \\ C & D \end{bmatrix}=\mathrm{det}\left( DA-CB \right)\]