import numpy as np
import pandas as pd
from numpy import linalg as LA
import math
from scipy.sparse import csr_matrix, diags
from scipy.sparse.linalg import spsolve, spsolve_triangular
from timeit import default_timer as timer

# Visualisation libraries

## plotly
from plotly.offline import init_notebook_mode, iplot 
import plotly.graph_objs as go
import plotly.offline as py
from plotly.subplots import make_subplots
import plotly.express as px

## seaborn
import seaborn as sns

## matplotlib
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse, Polygon
from matplotlib.font_manager import FontProperties
import matplotlib.colors as mcolors
from matplotlib import cm
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
plt.style.use('seaborn-whitegrid')
plt.rcParams['axes.labelsize'] = 14
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12
plt.rcParams['text.color'] = 'k'
%matplotlib inline

import warnings
warnings.filterwarnings("ignore")

Home / Finite Difference Methods

The Wave Equation in Time Domain

The wave equation in time domain can be stated as follows, \begin{equation} \frac{\partial^2 u}{\partial t^2}=c^2 \nabla^2 u+f,\qquad x\in \Omega \subset \mathbb{R}^d,~t\in(T_1,T_2] \end{equation} where $\Delta$ is the Laplacian, $f$ is a forcing function (for example our source) and $c$ is the wave velocity at which the time and spatially varying wave $u$ propagates.

• One Dimensional Wave Equation in Time Domain
  • Second-Order Finite-Difference Methods
  • A Fourth-Order Finite-Difference Method
• Two Dimensional Wave Equation in Time Domain
  • Second-Order Finite-Difference Methods
  • A Fourth-Order Finite-Difference Method
• Three Dimensional Wave Equation in Time Domain
  • Second-Order Finite-Difference Methods
  • A Fourth-Order Finite-Difference Method

A Fourth-Order Finite-Difference Method for
the One Dimensional Wave Equation in Time Domain

A one-dimensional form of the wave equation presented can be found as follows, \begin{align} \begin{cases} \dfrac{\partial^2 u}{\partial t^2} = c^2(x) \dfrac{\partial^2 u}{\partial x^2}+f(x,t),&(x,t)\in (a,b)\times(T_1,T_2],\\ u(x,T_1)=g(x),&x\in [a,b],\\ \dfrac{\partial }{\partial t}u(x,T_1) = h(x),&x\in [a,b],\\ u(a,t) = f_a(t),&t \in[T_1,T_2],\\ u(b,t) = f_b(t),&t \in[T_1,T_2]. \end{cases} \end{align}

The following notations are introduced to generate an algorithm for the problem,

• $h=\Delta x=\frac{b-a}{N_x}$ and $\tau=\Delta t=\frac{T_2-T_1}{N_t}$ where $M$ and $N$ are positive integers,
• $u_{i}^{n}=u(ih,n\tau)$, $f_{i}^{n}=f(ih,n\tau)$ for $0\leq i \leq N_x$ and $0\leq n \leq N_t$,
• $c_{i}=c(ih)$ for $0\leq i \leq N_x$ and $\lambda=\frac{\tau}{h}$.

A higher-order scheme can be introduced as follows (the readers are encouraged to see [3], \begin{align} \label{eq1A5.01}\left(1+\frac{1^2}{12}\delta_t^2\right)^{-1}\frac{\delta_t^2}{\tau^2} u_{i}^{n}= c^2_i \left(1+\frac{1}{12}\delta_x^2\right)^{-1}\frac{\delta_x^2}{h^2} u_{i}^{n}+f_{i}^{n}, \end{align} where it will be shown that this scheme is 4th-order.

In this section, the scheme is analyzed and also an numerical algorithm is generated. This equation can be also written in the following form, \begin{align} \left(1+\frac{1}{12}\delta_x^2\right)\frac{\delta_t^2}{\tau^2} \frac{u_{i}^{n}}{c^2_i}= \left(1+\frac{1}{12}\delta_t^2\right)\frac{\delta_x^2}{h^2} u_{i}^{n} +\left(1+\frac{1}{12}\delta_t^2\right)\left(1+\frac{1}{12}\delta_x^2\right)\frac{f_{i}^{n}}{c^2_i}. \end{align} The left hand side of the equation (13) can be expanded as follows, \begin{align} \left(1+\frac{1}{12}\delta_x^2\right)\frac{\delta_t^2}{\tau^2} \frac{u_{i}^{n}}{c^2_i}= \frac{1}{12\tau^2}\left[ \left(\frac{u_{i+1}^{n+1}}{c^2_{i+1}}+10 \frac{u_{i}^{n+1}}{c^2_{i}} + \frac{u_{i-1}^{n+1}}{c^2_{i-1}}\right) -2\left(\frac{u_{i+1}^{n}}{c^2_{i+1}}+10 \frac{u_{i}^{n}}{c^2_{i}} + \frac{u_{i-1}^{n}}{c^2_{i-1}}\right) +\left(\frac{u_{i+1}^{n-1}}{c^2_{i+1}}+10 \frac{u_{i}^{n-1}}{c^2_{i}} + \frac{u_{i-1}^{n-1}}{c^2_{i-1}}\right)\right]. \end{align}

It follows that, \begin{align} \left[\frac{u_{i+1}^{n+1}}{c^2_{i+1}}+10 \frac{u_{i}^{n+1}}{c^2_{i}} + \frac{u_{i-1}^{n+1}}{c^2_{i-1}}\right] &-2 \left[\frac{u_{i+1}^{n}}{c^2_{i+1}}+10 \frac{u_{i}^{n}}{c^2_{i}} + \frac{u_{i-1}^{n}}{c^2_{i-1}}\right] +\left[\frac{u_{i+1}^{n-1}}{c^2_{i+1}}+10 \frac{u_{i}^{n-1}}{c^2_{i}} + \frac{u_{i-1}^{n-1}}{c^2_{i-1}}\right]= \lambda^2\left[u_{i+1}^{n+1}-2 u_{i}^{n+1} + u_{i-1}^{n+1}\right] \notag\\ & +10\lambda^2\left[u_{i+1}^{n}-2u_{i}^{n} + u_{i-1}^{n}\right] +\lambda^2 \left[u_{i+1}^{n-1}-2u_{i}^{n-1} + u_{i-1}^{n-1}\right] +\tau^2\mathbf{F}_{i}^{n}. \end{align}

where $$\mathbf{F}_{i}^{n}= \frac{1}{144}\left[\frac{1}{c^2_{i+1}}\left(f_{i+1}^{n+1}+10 f_{i+1}^{n}+f_{i+1}^{n-1}\right) +\frac{10}{c^2_{i}}(f_{i}^{n+1}+10f_{i}^{n}+f_{i}^{n-1})+\frac{1}{c^2_{i-1}}(f_{i-1}^{n+1}+10f_{i-1}^{n}+f_{i-1}^{n-1})\right] $$ The equation (14) can be simplified as follows,

\begin{align} \left(\frac{1}{c^2_{i+1}}-\lambda^2\right)u_{i+1}^{n+1}&+\left(\frac{10}{c^2_{i}}+2\lambda^2\right)u_{i}^{n+1}+\left(\frac{1}{c^2_{i-1}}-\lambda^2\right)u_{i-1}^{n+1}= \left(\frac{2}{c^2_{i+1}}+10\lambda^2\right)u_{i+1}^{n}+\left(\frac{20}{c^2_{i}}-20\lambda^2\right)u_{i}^{n} +\left(\frac{2}{c^2_{i-1}}+10\lambda^2\right)u_{i-1}^{n} \notag \\ & -\left(\frac{1}{c^2_{i+1}}-\lambda^2\right)u_{i+1}^{n-1} -\left(\frac{10}{c^2_{i}}+2\lambda^2\right)u_{i}^{n-1}-\left(\frac{1}{c^2_{i-1}}-\lambda^2\right)u_{i-1}^{n-1} +\tau^2\mathbf{F}_{i}^{n}. \end{align}

Consider the following matrices, \begin{align} A=\begin{bmatrix} \frac{10}{c^2_{1}}+2\lambda^2 & \frac{1}{c^2_{2}}-\lambda^2 & 0 & 0 &\dots & \dots &0 \\ \frac{1}{c^2_{1}}-\lambda^2 & \frac{10}{c^2_{2}}+2\lambda^2 & \frac{1}{c^2_{3}}-\lambda^2 & 0 & 0 &\dots &\vdots \\ 0 & \frac{1}{c^2_{2}}-\lambda^2 & \frac{10}{c^2_{3}}+2\lambda^2 & \frac{1}{c^2_{4}}-\lambda^2 & 0 & \dots &\vdots \\ 0 & 0& \ddots & \ddots & \ddots & 0 & \vdots \\ \vdots & \dots &0 & \frac{1}{c^2_{M-4}}-\lambda^2 & \frac{10}{c^2_{M-3}}+2\lambda^2 & \frac{1}{c^2_{N_x-2}}-\lambda^2 & 0\\ \vdots & \dots & 0 &0 & \frac{1}{c^2_{M-3}}-\lambda^2 & \frac{10}{c^2_{N_x-2}}+2\lambda^2 & \frac{1}{c^2_{N_x-1}}-\lambda^2 \\ 0 & \dots & \dots & 0 &0 & \frac{1}{c^2_{N_x-2}}-\lambda^2 & \frac{10}{c^2_{N_x-1}}+2\lambda^2\\ \end{bmatrix}, \end{align} \begin{align} B=\begin{bmatrix} \frac{20}{c^2_{1}}-20\lambda^2 & \frac{2}{c^2_{2}}+10\lambda^2 & 0 & 0 &\dots & \dots &0 \\ \frac{2}{c^2_{1}}+10\lambda^2 & \frac{20}{c^2_{2}}-20\lambda^2 & \frac{2}{c^2_{3}}+10\lambda^2 & 0 & 0 &\dots &\vdots \\ 0 & \frac{2}{c^2_{2}}+10\lambda^2 & \frac{20}{c^2_{3}}-20\lambda^2 & \frac{2}{c^2_{4}}+10\lambda^2 & 0 & 0 &\vdots \\ 0 & 0& \ddots & \ddots & \ddots & 0 & \vdots \\ \vdots & \dots &0 & \frac{2}{c^2_{M-4}}+10\lambda^2 & \frac{20}{c^2_{M-3}}-20\lambda^2 & \frac{2}{c^2_{N_x-2}}+10\lambda^2 & 0\\ \vdots & \dots & 0 &0 & \frac{2}{c^2_{M-3}}+10\lambda^2 & \frac{20}{c^2_{N_x-2}}-20\lambda^2 & \frac{2}{c^2_{N_x-1}}+10\lambda^2 \\ 0 & \dots & \dots & 0 &0 & \frac{2}{c^2_{N_x-2}}+10\lambda^2 & \frac{20}{c^2_{N_x-1}}-20\lambda^2\\ \end{bmatrix}, \end{align} \begin{align} C=-A \end{align} and also \begin{align} \textbf{u}^n=\begin{bmatrix} u_{1}^{n} \\ u_{2}^{n}\\ \vdots\\ u_{N_x-2}^{n}\\ u_{N_x-1}^{n} \end{bmatrix},~ \textbf{b}^n= \begin{bmatrix} -\left(\frac{1}{c^2_{0}}-\lambda^2\right)f_{a}^{n+1}+\left(\frac{2}{c^2_{0}}+10\lambda^2\right)f_{a}^{n} -\left(\frac{1}{c^2_{0}}-\lambda^2\right)f_{a}^{n-1}\\ 0\\ \vdots\\ 0\\ -\left(\frac{1}{c^2_{N_x}}-\lambda^2\right)f_{b}^{n+1}+\left(\frac{2}{c^2_{N_x}}+10\lambda^2\right)f_{b}^{n} -\left(\frac{1}{c^2_{0}}-\lambda^2\right)f_{b}^{n-1} \\ \end{bmatrix}, \textbf{F}^n=\tau^2 \begin{bmatrix} \mathbf{F}_{1}^{n}\\ \mathbf{F}_{2}^{n}\\ \vdots\\ \mathbf{F}_{N_x-2}^{n}\\ \mathbf{F}_{N_x-1}^{n} \end{bmatrix}. \end{align} Therefore, the presented scheme can be expressed in the following matrix form, \begin{align} \textbf{u}^{n+1}=A^{-1}\left(B \textbf{u}^{n}+C \textbf{u}^{n-1}+\textbf{b}^n+ \textbf{F}^n\right). \end{align}

For Stability and convergence analyses, see this file.

Example

Consider the following problem, \begin{align*} \frac{\partial^2 u}{\partial t^2} &= \mathrm{e}^{-x} \frac{\partial^2 u}{\partial x^2}+{\mathrm{e}}^{t+x}-{\mathrm{e}}^t ,&(x,t)\in (0,1)\times(0,1],\\ u(x,0)&=\mathrm{e}^{x},&x\in [0,1],\\ \frac{\partial }{\partial t}u(x,0) &= \mathrm{e}^{x},&x\in [0,1],\\ u(0,t) &= \mathrm{e}^{t} &t \in[0,1],\\ u(1,t) &= \mathrm{e}^{t+1},&t \in[0,1]. \end{align*} The exact solution corresponding to the problem can be found as follows, $$ u(x,t)=\mathrm{e}^{x+t} $$

The Fourth-order Scheme

def Figs(U_exact, U_Comp, Scheme, data, y_range = None):
    fig = go.Figure()
    # [a,b] and [T1,T2]
    a = 0.0
    b = 1.0
    T1 = 0.0
    T2 = 1.0
    x = np.linspace(a, b, U_exact.shape[0])
    fig.add_trace(go.Scatter(x = x, y = U_Comp, mode='lines', name='Approximated'))
    fig.add_trace(go.Scatter(x = x, y = U_exact, mode='lines', name='Exact'))
    fig['layout']['xaxis'].update(range=[x[0], x[-1]])
    if y_range != None :
        fig['layout']['yaxis'].update(range=y_range)
    fig.update_xaxes(showline=True, linewidth=1, linecolor='Lightgray', mirror=True,
                    showgrid=True, gridwidth=1, gridcolor='Lightgray',
                    title_text = '$x$')
    fig.update_yaxes(showline=True, linewidth=1, linecolor='Lightgray', mirror=True,
                    showgrid=True, gridwidth=1, gridcolor='Lightgray',
                    title_text = '$u(x,1)$')
    fig.update_xaxes(zeroline=True, zerolinewidth=1, zerolinecolor='Lightgray')
    fig.update_yaxes(zeroline=True, zerolinewidth=1, zerolinecolor='Lightgray')
    fig.update_layout(legend_orientation='v', plot_bgcolor= 'white', height= 600, width= 670)
    fig.update_traces(marker_line_color= 'Black', marker_line_width=0.5, opacity=1)
    fig.update_layout(title={'text':  r'%s Approximation with Nx = %i and Nt = %i' % 
                             (Scheme,data.Nx.values[-1], data.Nt.values[-1]),
                             'x':0.46, 'y':0.88, 'xanchor': 'center', 'yanchor': 'top'},
                      font=dict(size=10))
    fig.show()
    
def Example01_4thSolver1D(Nx,Nt):
    # Nx = the number of mesh points for space
    # Nt = the number of mesh points for time
    # [a,b] and [T1,T2]
    a = 0.0
    b = 1.0
    T1 = 0.0
    T2 = 1.0
    # Functions
    C2 =lambda x: np.exp(-x)
    g=lambda x: np.exp(x)
    h=lambda x: np.exp(x)
    fa=lambda t: np.exp(t)
    fb=lambda t: np.exp(t+1)
    f=lambda x,t: np.exp(t + x) - np.exp(t)
    Ue=lambda x,t: np.exp(x+t)
    # Additional functions
    gxx=lambda x: np.exp(x);
    gxxxx=lambda x: np.exp(x);
    hxx=lambda x: np.exp(x);
    ft=lambda x,t: np.exp(t + x) - np.exp(t);
    ftt=lambda x,t: np.exp(t + x) - np.exp(t);
    # dx, dt and lambda^2
    dx =(b-a)/Nx
    dt =(T2-T1)/Nt
    l2=(dt/dx)**2
    # discretizing [a,b] and [T1,T2]
    xx=np.linspace(a, b, Nx+1)
    tt=np.linspace(T1, T2, Nt+1)
    # dicretizing c^2(x)
    CC=C2(xx)
    # u initial time value
    u0 = g(xx);
    u=np.zeros(Nx+1, dtype=float)
    u1=np.zeros(Nx+1, dtype=float)
    # Computing the solution at the first time step
    GG= g(xx)+h(xx)*dt+0.5*(C2(xx)*gxx(xx)+f(xx,T1))*(dt**2)+\
    (1/6)*(C2(xx)*hxx(xx)+ft(xx,T1))*(dt**3)+(1/24)*(C2(xx)*gxxxx(xx)+ftt(xx,T1))*(dt**4)
    u1[1:Nx]=GG[1:Nx]
    u1[0]=fa(tt[1])
    u1[Nx]=fb(tt[1])
    # indeces
    mid=list(range(1,Nx))
    sub=list(range(1,Nx-1))
    sup=list(range(2,Nx))
    midp=[i+1 for i in mid]
    midm=[i-1 for i in mid]
    # Matrix A
    A=diags([(1/CC[sub])-l2,(10/CC[mid])+2*l2,(1/CC[sup])-l2], [-1,0,1])
    # Matrix B
    B=diags([(2/CC[sub])+10*l2,(20/CC[mid])-20*l2,(2/CC[sup])+10*l2], [-1,0,1])
    # Matrix C
    C=-A
    # Fa and Fb
    FA= fa(tt)
    FB= fb(tt)
    for n in range(1,Nt):
        FF=((dt**2)/12)*((1/CC[midp])*(f(xx[midp],tt[n+1])+10*f(xx[midp],tt[n])+f(xx[midp],tt[n-1]))+\
        (10/CC[mid])*(f(xx[mid],tt[n+1])+10*f(xx[mid],tt[n])+f(xx[mid],tt[n-1]))+\
        (1/CC[midm])*(f(xx[midm],tt[n+1])+10*f(xx[midm],tt[n])+f(xx[midm],tt[n-1])))
        u=np.zeros(Nx+1, dtype=float)
        bb=np.zeros(len(mid), dtype=float)
        bb[0]=-(1/CC[0]-l2)*FA[n+1]+(2/CC[0]+10*l2)*FA[n]-(1/CC[0]-l2)*FA[n-1]
        bb[Nx-2]=-(1/CC[Nx]-l2)*FB[n+1]+(2/CC[Nx]+10*l2)*FB[n]-(1/CC[Nx]-l2)*FB[n-1]
        u[mid]=spsolve(A,B*u1[mid]+C*u0[mid]+bb+FF)
        u[0]=fa(tt[n+1])
        u[Nx]=fb(tt[n+1])
        u0=u1
        u1=u
        del FF
    # end of the loop
    # Norm
    U_Comp = u
    U_exact = Ue(xx,tt[-1])
    Norm=LA.norm(U_Comp - U_exact, np.inf)
    return Norm, U_exact, U_Comp

it=4
Norm=np.zeros(it, dtype=float)
Nx0=10;
Nt0=200;
Nx=np.asarray([Nx0*2**n for n in range(0,it)])
Nt=np.asarray([Nt0 for n in range(0,it)])
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)
CPU_Time=np.zeros(it, dtype=float)
# iteration
for n in range(it):
    if n != (it-1):
        start = timer()
        Norm[n],_,_ =Example01_4thSolver1D(Nx[n],Nt[n])
        CPU_Time[n] = timer() - start
    else:
        start = timer()
        Norm[n], U_exact, U_Comp =Example01_4thSolver1D(Nx[n],Nt[n])
        CPU_Time[n] = timer() - start
    if (n>0):
        Ratio[n]=Norm[n-1]/Norm[n]
        LOG[n]=math.log(Ratio[n],2)
        
data = pd.DataFrame({'Nx': Nx+1, 'Nt': Nt+1, 'Norm': Norm, 'Ratio': Ratio, '$\log_2$': LOG, 'CPU Time': CPU_Time})
del it, Norm, Nx0, Nt0, Nx, Nt, Ratio, LOG, CPU_Time, n, start
display(data.style.format({'Norm': "{:.4e}"}))
Figs(U_exact, U_Comp, Scheme = 'Explicit Scheme', data = data)

	Nx	Nt	Norm	Ratio	$\log_2$	CPU Time
0	11	201	2.4019e-07	0.000000	0.000000	0.060147
1	21	201	1.5106e-08	15.899709	3.990928	0.060584
2	41	201	9.3612e-10	16.137231	4.012321	0.063921
3	81	201	6.0715e-11	15.418350	3.946576	0.070082

it=4
Norm=np.zeros(it, dtype=float)
Nx0=10;
Nt0=20;
Nx=np.asarray([Nx0*2**n for n in range(0,it)])
Nt=np.asarray([Nt0*2**n for n in range(0,it)])
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)
CPU_Time=np.zeros(it, dtype=float)
# iteration
for n in range(it):
    if n != (it-1):
        start = timer()
        Norm[n],_,_ =Example01_4thSolver1D(Nx[n],Nt[n])
        CPU_Time[n] = timer() - start
    else:
        start = timer()
        Norm[n], U_exact, U_Comp =Example01_4thSolver1D(Nx[n],Nt[n])
        CPU_Time[n] = timer() - start
    if (n>0):
        Ratio[n]=Norm[n-1]/Norm[n]
        LOG[n]=math.log(Ratio[n],2)
        
data = pd.DataFrame({'Nx': Nx+1, 'Nt': Nt+1, 'Norm': Norm, 'Ratio': Ratio, 'Log2': LOG, 'CPU Time': CPU_Time})
del it, Norm, Nx0, Nt0, Nx, Nt, Ratio, LOG, CPU_Time, n, start
display(data.style.format({'Norm': "{:.4e}"}))
Figs(U_exact, U_Comp, Scheme = 'Explicit Scheme', data = data)

	Nx	Nt	Norm	Ratio	Log2	CPU Time
0	11	21	2.4176e-07	0.000000	0.000000	0.006511
1	21	41	1.5199e-08	15.906608	3.991554	0.012664
2	41	81	9.5222e-10	15.961532	3.996527	0.026077
3	81	161	5.9392e-11	16.032825	4.002957	0.057179

Convergence Rate

For a fixed $N_t$, we have

it=4
Norm=np.zeros(it, dtype=float)
Nx0=10
Nt0=200
Nx=np.asarray([Nx0*2**n for n in range(0,it)])
Nt=np.asarray([Nt0 for n in range(0,it)])
#
Exact=np.zeros((it,Nx0+1), dtype=float)
Fourth_order_Scheme=np.zeros((it,Nx0+1), dtype=float)

#
for n in range(it):
    _,U_exact, U_Comp = Example01_4thSolver1D(Nx[n],Nt[n])
    Exact[n,:] = U_exact[0:(Nx[n]+1):(2**n)]
    Fourth_order_Scheme[n,:] = U_Comp[0:(Nx[n]+1):(2**n)]
    del U_exact, U_Comp

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(16, 6))
for n in range(it):
    tag= "$Nx = %d$" % (Nx[n])
    _ = ax.loglog(np.linspace(0, 1, Nx0+1), np.abs(Exact[n,:]-Fourth_order_Scheme[n,:]),label=tag)
    _ = ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0., fontsize = 14)
    _ = ax.set_title('Fourth-Order Scheme', fontsize = 14)

---

References

1. Das, Sambit, Wenyuan Liao, and Anirudh Gupta. "An efficient fourth-order low dispersive finite difference scheme for a 2-D acoustic wave equation." Journal of computational and Applied Mathematics 258 (2014): 151-167.

2. Liao, Wenyuan. "On the dispersion, stability and accuracy of a compact higher-order finite difference scheme for 3D acoustic wave equation." Journal of Computational and Applied Mathematics 270 (2014): 571-583.

3. Liao, Wenyuan. "An implicit fourth-order compact finite difference scheme for one-dimensional Burgers’ equation." Applied Mathematics and Computation 206.2 (2008): 755-764.

---