Finite difference method for linear equation

We want to introduce the basic idea of the FDM for the linear equation, including the concept, schemes, global error and convergence.

Model problem

Consider the one-dimensional linear equation\[ \begin{align*} u_t+ au_x & = 0,\quad -\infty < x<\infty , t\geq 0\\ u(x,0) &= u_0(x). \end{align*} \tag{1}\]Here, \(a\neq 0\) is constant denoting the moving speed of the fluid. A fundamental result implies that the solution is \(u(x,t)=u_0(x-at),t>0\). As time evolves, the initial data moves to the right (\(a>0\)) or to the left (\(a<0\)) without any change.

We discrete the space-time domain by \((x_j,t_n)\) with\[ x_j = jh, \quad t_n=nk.\]The point-wise values of the exact solution will be denoted by \(u_j^n = u(x_j,t_n)\), and the discrete values are represented by \(U_j^n\).

Finite difference scheme

For brevity, we only consider two-level methods.\[ (\text{Forward Euler})\quad \frac{U_j^{n+1}-U_j^{n}}{k} + a\frac{U_{j+1}^n-U_{j-1}^n}{2h} = 0. \tag{2a}\]
\[ (\text{Backward Euler})\quad \frac{U_j^{n+1}-U_j^{n}}{k} + a\frac{U_{j+1}^{n+1}-U_{j-1}^{n+1}}{2h} = 0. \tag{2b}\]
A classical method is the Lax-Friedrichs scheme, which is obtained by replacing \(U^n_j\) by \((U_{j-1}^{n}+U_{j+1}^{n})/2\) and is stable provided \(k/h\) is small enough, as we will see later.\[ (\text{Lax-Friedrichs})\quad \frac{U_j^{n+1}-(U_{j-1}^{n}+U_{j+1}^{n})/2}{k} + a\frac{U_{j+1}^{n}-U_{j-1}^{n}}{2h} = 0. \tag{2c}\]
For brevity, we introduce the (linear) iteration operator for these methods:\[ U^{n+1} = \mathcal{H}_k(U^n). \tag{3}\]For example, the iteration matrix of the Lax-Friedrichs scheme is \(\mathcal{H}_k(U^n,j) = \frac{1}{2}(U_{j-1}^{n}+U_{j+1}^{n}) – \frac{ak}{2h}\left( U_{j+1}^{n+1}-U_{j-1}^{n+1} \right)\).

Error and norms

Define point-wise error\[ E_j^n = U_j^n-u_j^n.\]To measure the “distance” between the exact solution \(u\) and the numerical solution, we define the following norms\[ \begin{align*} \|U^n\|_{\infty} &= \sup_{j} |U_j^n|,\\ \|U^n\|_{1} &= \sum_{j} h|U_j^n|,\\ \|U^n\|_{2} &= \left( \sum_{j} h|U_j^n|^2 \right)^{1/2}. \end{align*}\]

In general, we approximate the continuous equation \(\mathcal{L}[u]=f\) by discrete equation\[ \mathcal{L}U_{j}^{n} =f_j^n. \tag{4}\]The error arises from the left, so it is natural to define the local truncation error (LTE) by substituting \(u\) into (4)\[ \tau_{j}^n = \mathcal{L}u_j^n-f_j^n. \tag{5a}\]For example, the forward Euler method has the LTE\[ \mathcal{L}u_j^n = \frac{1}{k}(k u_t(x_j,t_n)+O(k^2)) + \frac{a}{2h} (2hu_x(x_j,t_n)+O(h^3)) = O(k)+O(h^2)\]

Note. For a fixed spatial mesh \(h\) or fixed relation between \(k\) and \(h\) (f.e.,\(k/h=const\)), we only care the error behaviour with respect to time step \(k \to 0\). In this case, we can define the LTE of an explicit method as\[ \tau_{j}^n = \frac{1}{k}\left[ u(x_j,t_n+k)-\mathcal{H}_k(u(x_j,t_n)) \right]. \tag{5b}\]On the other hand, implicit method has a more complicated form\[ B_1 U^{n+1} = B_0 U^n+ k \tau^n \quad \Rightarrow \quad U^{n+1} = B_1^{-1}B_0 U^n+ k B_1^{-1} \tau^n,\]where \(B_1^{-1}B_0\) is the iteration operator \(\mathcal{H}_k\) in (3). Implicit methods, in general, lead to a more difficult convergence analysis. For generality, we consider the LTE equation for iteration process (3):\[ u^{n+1} = \mathcal{H}_k(u^n) + k B_1^{-1}\tau^n. \tag{3b}\]

Definition (Consistency). Given a refining path, the method is consistent if LTE satisfies\[ \|\tau^n\|\to 0 \quad \text{as}\quad k\to 0.\]

For example, (2a) is consistent w.r.t. \(\|\cdot \|_{\infty}\) with refining path \(k/ h^2=const\).

Stability

Now, we need to consider the key problem: convergence or error estimate. Intuitively, convergence means the difference between the exact and numerical solutions tends to zero w.r.t. some norm as \(k, h\to 0\). We also want to know at what speed this difference tends to zero. In many textbooks, it says that convergence is not easy to show, but we can instead show stability, as the Lax equivalence theorem states, to obtain convergence. In fact, we have

Lax Equivalence Theorem. For a consistent linear method, stability is necessary and sufficient for convergence.

Now, we first turn to stability and after that, we come back to convergence. Combing (3) and (3b), we obtain error equation at time level \(n\)\[ E^{n+1} = \mathcal{H}_k(E^n) – k B_1^{-1}\tau^n, \tag{6}\]which further implies by recursion\[ E^{n+1} = \mathcal{H}_k^n(E^0) – \sum_{j=0}^{n} k \mathcal{H}_k^j(B_1^{-1}\tau^{n-j}). \tag{6b}\]

Definition (Stability). A method is said to be stable if for each time \(T\) there is a constant \(C_s\) and a value \(k_0>0\) s.t.\[ \|\mathcal{H}_k^n\|\leq C_s \quad \forall nk\leq T,\ k < k_0.\]

In particular, the method is stable if \(\|\mathcal{H}_k\|\leq 1\). More generally, slight growth is allowed in the sense that \(\|\mathcal{H}_k\|\leq 1+\alpha k\) for all \(k< k_0\). In fact, \(\|\mathcal{H}_k^n\|\leq (1+\alpha k)^n \leq e^{\alpha nk}\leq e^{\alpha T} \).

Assume now the method is stable, that we have the following error estimate from (6b)\[\begin{align*} \|E^{n+1}\| &\leq \|\mathcal{H}_k^{n}\| \|E^0\| + k\|B_1^{-1}\| \sum_{j=0}^{n} \|\mathcal{H}_k^j\| \|\tau^{n-j}\| \\ &\leq C_s \|E^0\| + k\|B_1^{-1}\| \cdot C_s \cdot O (k^p+h^p) \cdot n \\ &\leq C_s \left( \|E^0\| + T\|B_1^{-1}\| \cdot O (k^p+h^p) \right) .\end{align*}\]Clearly, error \(E^{n+1}\) for all \(n\) tends to zero as \(k,h\to 0\). Hence, for a consistent (two-level) scheme, stability is sufficient for convergence.

Example 1. Prove Lax-Friedrichs shceme is convergent w.r.t. \(\|\|_{1}\) provided\[ ak/h\leq 1. \tag{7}\]Proof. It is clear that for (7)\[ \begin{align*} \|U^{n+1}\|_1 &= h \sum_{j}|U^{n+1}_j| \\ &= h \sum_{j}|\frac{1}{2}(U_{j-1}^{n}+U_{j+1}^{n}) – \frac{ak}{2h}\left( U_{j+1}^{n+1}-U_{j-1}^{n+1} \right)| \\ &\leq \frac{h}{2} \sum_{j}\left\{ \left( 1+\frac{ak}{h} \right)|U_{j-1}^{n}| + \left( 1-\frac{ak}{h} \right) |U_{j+1}^{n+1}| \right\} \\ &= \frac{1}{2} \left\{ \left( 1+\frac{ak}{h} \right) \|U^n\| + \left( 1-\frac{ak}{h} \right) \|U^n\| \right\} \\ &= \|U^n\|, \end{align*}\]which implies the iteration matrix \(\|\mathcal{H}_k\|_1\leq 1\) and thus the scheme is convergent. \(\heartsuit\)

The CFL condition

In the last section, for the Lax-Friedrichs method we have seen that condition (7) is sufficient for stability. Now, we wonder what is a necessary condition for stability. As indicated in [1,2], a necessary stability condition for any numerical method is that the domain of dependence of the finite difference method should include the domain of dependence of the PDE, at least in the limit as \(k,h\to 0\). This condition is known as the CFL condition after Courant, Friedrichs and Lewy.

In linear equation (1), the domain of dependence \(\mathcal{D}(\bar{x},\bar{t})\) consists of point \(\bar{x}-a \bar{t}\), because the wave moves from \(x-at\) at time \(0\) to \(x\) at \(t\). Now, consider Lax-Friedrichs method, we see \(\mathcal{D}_k(x_j,t_n)= \left\{ x:\ |x-x_j|\leq nh \right\} \), precisely,\[ \mathcal{D}_k(\bar{x},\bar{t})= \left\{ x:\ |x-\bar{x}|\leq (\bar{t}/k)h \right\}. \tag{8}\]Now, on what condition there is \(\bar{x}-a \bar{t} \in \mathcal{D}_k(\bar{x},\bar{t})\)? Substituting \(x=\bar{x}-a \bar{t}\) in (8) leads to\[ |a \bar{t}|\leq \frac{\bar{t}}{k}h \quad \Rightarrow \quad \frac{ak}{h}\leq 1.\]This is exactly (7), which is also necessary for stability.

References

[1] R. J. LeVeque. Numerical Methods for Conservation Laws. 1992.
[2] R. Courant, K. O. Friedrichs and H. Levy. On the partial difference equations of mathematical physics, IBM Journal, 11 (1967), pp. 215-234.

Appendix A. Python code

import numpy as np
import matplotlib.pyplot as plt
from numba import jit

# # pde example I=(0,L)
@jit
def u0(x):
    u=np.ones(x.size)
    for j in range(u.size):
        if x[j]>0:
            u[j]=0.
    return u


class Wave1d(object):
    def __init__(self):
        self.a = 1.
        self.L = 1.
        self.T = .5

    def u0(self,x):
        return u0(x)

    def u(self,x,t):
        return self.u0(x-self.a*t)
    
    def f(self,x,t):
        pass
    
    def plot(self):
        t = np.linspace(0, self.T, 50)
        x = np.linspace(0, self.L, 20)

        plt.figure(figsize=(8, 6), dpi=80)
        plt.ion()

        for i in range(t.shape[0]):
            plt.cla()
            plt.title("Wave")
            plt.grid(True)
            u = self.u(x,t[i])

            plt.plot(x,u, '-r', label="Exact solution: t="+str(int(100*t[i])/100))
            # plt.plot(v1*t[i], 0, 'or' , markersize=12)

            plt.xlim(x[0], x[-1])
            # plt.plot(v2*t[i], 0, 'ob' , markersize=12)
            plt.ylim(-1.2, 1.2)
            plt.legend(loc="upper left", shadow=True)
            plt.pause(.2)

        plt.ioff()
        plt.show()


# uniform spatial and time space
class Mesh():
    def __init__(self,x,t):
        self.x = x
        self.t = t


# solver class
class Wave1DSolver():
    def __init__(self, pde, mesh):
        self.pde = pde
        self.mesh = mesh
        self.Nx = mesh.x.shape[0]-1
        self.Nt = mesh.t.shape[0]-1
        self.dx = mesh.x[1]-mesh.x[0]
        self.dt = mesh.t[1]-mesh.t[0]
        self.CFL = (pde.a*self.dt/self.dx)
        self.u_fdm = None
        
        self.u_exact = pde.u(mesh.x, pde.T)
        self.err_l2 = 0.
        self.err_l2_final = 0.

    def solve(self):
        dt, dx, CFL = self.dt, self.dx, self.CFL
        Nx = self.Nx
        uplus = np.zeros(self.Nx+1)
        u = np.zeros(self.Nx+1)
        err_l2 = 0.
        
        # set initial values
        x = self.mesh.x
        u = self.pde.u0(x)

        err_l2 += dt*dx * np.sum((u-self.pde.u(x,self.mesh.t[0]))**2)
        
        # gif plot - 1/3
        plt.figure(figsize=(8, 6), dpi=80)
        plt.ion()
        
        for j in range(1,self.Nt):
            # update_LaxFriedrichs(uplus, u, CFL)
            update_Upwind(uplus, u, CFL)
            # update_LaxWendroff(uplus, u, CFL)
            # update_BeamWarming(uplus, u, CFL)
                
            uplus[0]=1
            uplus[Nx]=0
            
            err_l2 += dt*dx * np.sum((uplus-self.pde.u(x,self.mesh.t[j+1]))**2)
            self.err_l2_final = np.sqrt( dx * np.sum((uplus-self.pde.u(x,self.mesh.t[j+1]))**2) )
            
            # gif plot - 2/3
            if j%np.ceil(t.shape[0]/20) == 1:
                plt.cla()
                plt.title("t="+str(int(1e4*self.mesh.t[j+1])/1e4))
                plt.grid(True)
                plt.plot(x,uplus, '-r', label="FDM solution")
                plt.plot(x,self.pde.u(x,self.mesh.t[j]), '--', label="Exact solution")
                plt.xlim(x[0], x[-1])
                plt.ylim(-0.5, 1.5)
                plt.legend(loc="upper left", shadow=True)
                plt.pause(.5)
            
            # Switch before next iteration
            u[:] = uplus
        # gif plot - 3/3
        plt.ioff()
        plt.show()
        
        self.u_fdm = uplus
        self.err_l2 = np.sqrt(err_l2)
        return 0
    
    def plot(self):
        plt.plot(self.mesh.x,self.u_fdm-self.u_exact, '-or', label='Error=FDM-Exact')
        # plt.plot(self.mesh.x, self.u_exact, '-db', label='Exact solution')
        plt.legend()
        plt.show()


@jit
def update_LaxFriedrichs(uplus, u, CFL):
    for i in range(1,u.size-1):
        uplus[i] = .5*(u[i-1]+u[i+1]) - .5*CFL*(u[i+1]-u[i-1])


@jit
def update_Upwind(uplus, u, CFL):
    for i in range(1,u.size-1):
        uplus[i] = CFL*u[i-1] + (1-CFL)*u[i]

@jit
def update_LaxWendroff(uplus, u, CFL):
    for i in range(1,u.size-1):
        uplus[i] = u[i]-.5*CFL*(u[i+1]-u[i-1]) + .5*CFL**2*(u[i-1]-2*u[i]+u[i+1])

@jit
def update_BeamWarming(uplus, u, CFL):
    i=1
    uplus[i] = u[i]-.5*CFL*(u[i+1]-u[i-1]) + .5*CFL**2*(u[i-1]-2*u[i]+u[i+1])
    for i in range(2,u.size-1):
        uplus[i] = u[i]-.5*CFL*(3*u[i]-4*u[i-1]+u[i-2]) + .5*CFL**2*(u[i]-2*u[i-1]+u[i-2])



if __name__ == "__main__":
    pde = Wave1d()
    # pde.plot()
    Nx = 40
    x = np.linspace(0, pde.L, num=Nx+1, endpoint=True)
    CFL = .1
    dt = CFL*(pde.L/Nx)/pde.a; print(dt)
    Nt = int(pde.T/dt+1); print(Nt)
    t = np.linspace(0, pde.T, num=Nt, endpoint=True)
    mesh = Mesh(x, t)

    solver = Wave1DSolver(pde, mesh)
    solver.solve()
    print(solver.err_l2, solver.err_l2_final)
    # solver.plot()