Stability#

Just like for the finite difference method we need to be very careful with the time step $\mathrm{\Delta }t$, in order to obtain a stable solution. A stable solution is one that remains bounded as $t⟶\mathrm{\infty }$. However, there is one major difference from the finite difference or collocation methods. The Galerkin method described above does not involve a mesh size $\mathrm{\Delta }x$! So the Courant number is not defined. How do we overcome this problem? In order to describe this problem it is helpful to first revisit the stability of the finite difference method.

import numpy as np
from scipy import sparse

N = 11
D2 = sparse.diags((1, 1, -2, 1, 1), (-N, -1, 0, 1, N), shape=(N+1, N+1))
Lambda = np.linalg.eig(D2.toarray())[0]
print(f"Minimum eigenvalue = {min(Lambda)}")

Minimum eigenvalue = -4.0

The diffusion equation#

Consider the diffusion equation with Dirichlet boundary conditions and a prescribed initial condition

In order to solve (75) we choose a global Galerkin method with mapped Legendre polynomials and basis functions ${\psi }_i(x)=P_i(X)-P_{i+2}(X)$ such that $V_N=\text{span}\{{\psi }_i{\}}_{i=0}^N$. We then use Eq. (73) and create a linear algebra problem by inserting for $u_N^{n+1},u_N^n$ and $v={\psi }_i$ and using integrating by parts on the diffusion term:

Since we are using basis functions that are all zero on the boundaries, the last term disappears. In order to incorporate the Dirichlet boundary condition we add a boundary term $B(x)$ and solve for ${\tilde{u}}_N^{n+1}=u_N^{n+1}-B(x)$, where ${\tilde{u}}_N^{n+1}(0,t)={\tilde{u}}_N^{n+1}(L,t)=0$. Note that we get exactly the same equation for $u_N^{n+1}$ as for ${\tilde{u}}_N^{n+1}$ since all terms with $B(x)$ will be eliminated when inserting $u_N^{n+1}={\tilde{u}}_N^{n+1}+B(x)$ and $u_N^n={\tilde{u}}_N^n+B(x)$. But we solve the equation for ${\tilde{u}}_N^{n+1}$, with the zero boundary conditions, and then to get $u_N^{n+1}$ we simply add $B(x)$.

We use Legendre basis functions with a mapping from $x\in [0,L]$ to $X\in [-1,1]$ such that

The stiffness matrix $s_{ij}=({\psi }_j^{\prime },{\psi }_i^{\prime })$ was given in Eq. (51), but it needs to be modified for the mapping

The mass matrix $a_{ij}=({\psi }_j,{\psi }_i)$ is similarly

The mass matrix is symmetric and we know that $(P_j,P_i)=\Vert P_i{\Vert }^2{\delta }_{ij}$, where $\Vert P_i{\Vert }^2=\frac{2}{2i+1}$ is the squared $L^2$ norm of $P_i$. Hence with $j=i$ we get

and with $j=i+2$ we get

The mass matrix is as such

On matrix form the equation to solve is

This equation is solved for $n=0,1,\dots ,N_t-1$, after first initializing ${\hat{\mathit{u}}}^0$. A solver is created in the class FE_Diffusion below, where we reuse the code from mandatory assignment 2. In the code we assemble the mass and stiffness matrices $A$ and $S$, and compute the eigenvalues $\lambda$ of $H=A^{-1}S$. We then use that $|1-\lambda \mathrm{\Delta }t|\le 1$ (where the minus is because we have used integration by parts) and the fact that here all eigenvalues are positive real numbers, to determine that

This means that for stability we need to choose a time step $\mathrm{\Delta }t\le 2/max(\lambda )$.

import numpy as np
import sympy as sp
from scipy import sparse
from galerkin import TrialFunction, TestFunction, DirichletLegendre, DirichletChebyshev, \
    project, inner

x, t, c, L = sp.symbols('x,t,c,L')

class FE_Diffusion:
    """Class for solving the diffusion equation with

    Parameters
    ----------
    N : int
        Number of basis functions
    L0 : number
        The extent of the domain, which is [0, L0]
    bc : 2-tuple of numbers
        The Dirichlet boundary conditions at the two edges
    u0 : Sympy function of x, t, c and L
        Used for specifying initial condition
    """
    def __init__(self, N, L0=1, u0=sp.cos(sp.pi*x/L)+sp.cos(10*sp.pi*x/L)):
        self.N = N
        self.L = L0
        self.u0 = u0.subs({L: L0})
        self.bc = (self.u0.subs({x: 0, t: 0}).n(), self.u0.subs({x: L0, t: 0}).n())
        self.uh_np1 = np.zeros(N+1) # expansion coefficients \hat{u}^{n+1}
        self.uh_n = np.zeros(N+1)
        self.V = self.get_function_space()
        self.A, self.S = self.get_mats()
        self.lu = sparse.linalg.splu(self.A.tocsc())

    def get_function_space(self):
        return DirichletLegendre(self.N, domain=(0, self.L), bc=self.bc)

    def get_mats(self):
        """Return mass and stiffness matrices
        """
        N = self.N
        u = TrialFunction(self.V)
        v = TestFunction(self.V)
        a = inner(u, v)
        s = 2/self.L*sparse.diags([4*np.arange(N+1)+6], [0], shape=(N+1, N+1), format='csr')
        return a, s

    def get_eigenvalues(self):
        """Return eigenvalues of H=A^{-1}S"""
        return np.linalg.eig(np.linalg.inv(self.A.toarray()) @ self.S.toarray())[0]

    def __call__(self, Nt, dt=0.1, save_step=100):
        """Solve diffusion equation

        Parameters
        ----------
        Nt : int
            Number of time steps
        dt : number
            timestep
        save_step : int, optional
            Save solution every save_step time step

        Returns
        -------
        Dictionary with key, values as timestep, array of solution
        The number of items in the dictionary is Nt/save_step, and
        each value is an array of length N+1

        """
        # Initialize
        self.uh_n[:] = project(self.u0, self.V) # unm1 = u(x, 0)
        xj = self.V.mesh()
        plotdata = {0: self.V.eval(self.uh_n, xj)}

        # Solve
        f = np.zeros(self.N+1)
        for n in range(2, Nt+1):
            f[:] = self.A @ self.uh_n - dt*(self.S @ self.uh_n)
            self.uh_np1[:] = self.lu.solve(f)
            self.uh_n[:] = self.uh_np1
            if n % save_step == 0: # save every save_step timestep
                plotdata[n] = self.V.eval(self.uh_np1, xj)
        return plotdata

Lets solve the problem using initial condition $u(x,0)=\cos (\pi x/L)+\cos (10\pi x/L)$ and $dt=2/max(\lambda )$:

from utilities import plot_with_offset
sol = FE_Diffusion(40, L0=2)
dt = 2/max(sol.get_eigenvalues())
data = sol(1000, dt=dt, save_step=100)
plot_with_offset(data, sol.V.mesh(), figsize=(4, 3))

The solution is nice and smooth and the time step is

print(f'Maximum time step = {dt}')

Maximum time step = 2.1980578790345177e-05

Lets solve it using a slightly larger $dt$ to see if this really is the largest possible time step:

data = sol(1000, dt=1.01*dt, save_step=100)
plot_with_offset(data, sol.V.mesh(), figsize=(4, 3))

Obviously, the solver is now unstable and the solution diverges.

Changing from Legendre to Chebyshev makes the required time step even shorter.

class FE_Diffusion_Cheb(FE_Diffusion):

    def get_function_space(self):
        return DirichletChebyshev(self.N, domain=(0, self.L), bc=self.bc)

    def get_mats(self):
        """Return mass and stiffness matrices
        """
        N = self.N
        u = TrialFunction(self.V)
        v = TestFunction(self.V)
        a = inner(u, v)
        #s2 = -inner(u.diff(2), v) # slow
        k = np.arange(N+1, dtype=float)
        def diag(i):
            if i == 0:
                return 2*np.pi*(k+1)*(k+2)
            return 4*np.pi*(k[:-(i)]+1)
        diags = np.arange(0, N+1, 2)
        vals = [diag(i) for i in diags]
        s = 2/self.L * sparse.diags(vals, diags, shape=(N+1, N+1), format='csr')
        return a, s

sol = FE_Diffusion_Cheb(40, L0=2)
dt = 2/max(sol.get_eigenvalues())
data = sol(1000, dt=dt, save_step=100)
plot_with_offset(data, sol.V.mesh(), figsize=(4, 3))

print(f'Maximum time step = {dt}')

Maximum time step = 1.2332249161314778e-05

We can compare the absolute stability limits with the finite difference method, where $\mathrm{\Delta }t/\mathrm{\Delta }{x}^2\le 0.5$ for stability. Using a uniform mesh with $N=40$, we get $\mathrm{\Delta }x=2/N=0.05$ and $\mathrm{\Delta }t\le 0.5\cdot {0.05}^2=1.25\cdot {10}^{-3}$. Hence, the Galerkin method with Chebyshev basis functions requires approximately 100 times shorter time steps for stability!

Use a stencil matrix for the implementation#

It is rather messy to work with composite basis functions. A very smart trick is thus to use a stencil matrix $S$, such that (summation on repeated $j$ index)

This defines $S\in {\mathbb{R}}^{(N+1)\times (N+3)}$ for the Dirichlet problem as

With $\mathit{\psi }=\{{\psi }_i{\}}_{i=0}^N$ and $\mathit{P}=\{P_i{\}}_{i=0}^{N+2}$ (note the longer vector) we get all basis functions in matrix form:

But why is the stencil matrix so smart? Consider the mass matrix

You need to sort out the four (diagonal) matrices on the right. This is not very difficult, since they are all diagonal, but you need to watch out for the indices in order for the mass matrix to turn out correctly. Now use stencil matrix instead (with summation on $k$ and $l$):

Here you only need the one diagonal matrix $P=((P_l,P_k){)}_{k,l=0}^{N+2,N+2}$, where $(P_l,P_k)=\frac{2}{2l+1}{\delta }_{kl}$. The great advantage only becomes apparent in matrix form though

And this works for any stencil matrix, not just the Dirichlet one. That is the main advantage of the stencil matrix approach. For the Neumann basis (see lecture 11)

the stencil matrix is

where $\alpha =N(N+1)/((N+2)(N+3))$. The stiffness matrix can also be computed through the stencil matrix:

Again you need only the single matrix using orthogonal polynomials $q_{kl}=(P_l^{″},P_k)=-(P_l^{\prime },P_k^{\prime })$ and the matrix form is again much simpler:

Backward Euler#

The forward Euler method requires a very short time step for the diffusion equation. How about the backward Euler method

where the only difference from forward is that $(\mathcal{L}(u_N^{n+1}),v)$ is making use of the unknown $u_N^{n+1}$ and not the known $u_N^n$? The linear algebra problem is now

If we continue with the ansatz $u_N^{n+1}=g{u}_N^n$, and thus $u_N^{n+1}=g^{n+1}{u}_N^0$, we get

which on matrix form reads

If we multiply from left to right by $A^{-1}$ we get

Here we can once again use $H=A^{-1}S$ and $H{\hat{\mathit{u}}}^0=\lambda {\hat{\mathit{u}}}^0$ and thus we find

Since all time steps and eigenvalues are real numbers larger than 0, we find that $g$ is always less than 1 and the backward Euler method is as such unconditionally stable.

The Wave equation#

Lets consider the wave equation from lecture 5 in one spatial dimension

As you may recall, we solved the wave equation with Dirichlet, Neumann and open boundary conditions. The two time derivatives also required that the solution was initialized at two time-steps, or more mathematically precise, that both $u(x,0)$ and $\frac{\partial u}{\partial t}(x,0)$ were given.

We now use finite differences in time such that

A numerical residual ${\mathcal{R}}_N=\mathcal{R}(u_N^{n+1};u_N^n,u_N^{n-1})$ can now be defined as

where we use $\alpha =c\mathrm{\Delta }t$ and $u^{\prime }=\frac{\partial u}{\partial x}$ for simplicity, since there is only one spatial dimension. If we now initialize $u_N^0$ and $u_N^1$, the Galerkin method is to find $u_N^{n+1}\in V_N$ such that

for $n=1,\dots ,N_t-1$. We can also integrate the last term by parts

In order to solve the Galerkin problem we need to formulate the linear algebra problem. Hence we insert for $u_N^{n-1},u_N^n,u_N^{n+1}$ and $v={\psi }_i$ to obtain

We next rearrange such that the unknown ${\hat{u}}_j^{n+1}$ is on the left hand side and the rest of the known terms are on the right hand side

With matrix notation, $a_{ij}=({\psi }_j,{\psi }_i)$ and $s_{ij}=({\psi }_j^{\prime },{\psi }_i^{\prime })$, we get

where ${\mathit{b}}^n=((u_N^n{)}^{\prime }(L){\psi }_i(L)-(u_N^n{)}^{\prime }(0){\psi }_i(0){)}_{i=0}^N$. We can write this as

where the vector ${\mathit{f}}^n=A(2{\hat{\mathit{u}}}^n-{\hat{\mathit{u}}}^{n-1})-{\alpha }^2(S{\hat{\mathit{u}}}^n-{\mathit{b}}^n)$.