Why do I get an oscillatory solution when applying the implicit trapezoidal method to the linear diffusion equation?

Question

I wish to solve the following equation, $$\frac{\partial f}{\partial t}=\frac{\partial}{\partial x}\left(D(x)\frac{\partial f}{\partial x}\right)$$ using an exponential integrator. I discretize this equation in the following way, $$\frac{\partial f}{\partial t}=\frac{D_{i+\frac{1}{2}}\frac{f^{n}_{i+1}-f^{n}_{i}}{x_{i+1}-x_{i}}-D_{i-\frac{1}{2}}\frac{f^{n}_{i}-f^{n}_{i-1}}{x_{i}-x_{i-1}}}{x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}}$$ $$\implies\frac{\partial f}{\partial t}={\mathcal A}f^{n}$$ where $\mathcal A$ is the tridiagonal matrix. The above equation can be solved in the following way, $$f^{n+1}=e^{{\mathcal A}dt}f^{n}$$ which upon approximation gives the following, $$(I-{\mathcal A}dt)f^{n+1}\approx(I+{\mathcal A}dt)f^{n}$$ I then apply Thomas algorithm to compute $f^{n+1}$.

Below I am attaching a python code to solve $$\frac{\partial f}{\partial t}=\frac{\partial}{\partial x}\left(\frac{\partial f}{\partial x}\right)$$ but with the above mentioned algorithm I am getting oscillations.

import numpy as np
import matplotlib.pyplot as plt
import sys

def tridiag(n,a,b,c,d):
    E = np.zeros(n)
    F = np.zeros(n)
    new = np.zeros(n)
    for i in range(n):
        if (i==0):
            E[i] = -c[i]/b[i]
            F[i] = d[i]/b[i]
        else:
            E[i] = -c[i]/(a[i]*E[i-1] + b[i])
            F[i] = (d[i]-a[i]*F[i-1])/(a[i]*E[i-1]+b[i])
    new[n-1] = F[n-1]
    for i in range(n-2,-1,-1):
        new[i] = E[i]*new[i+1] + F[i]
    return new

def initial(x):
    return np.exp(-(x-5.0)**2/(2*0.01))

def D(x,alpha):
    return x**alpha

def A(x,alpha):
    return x**alpha

res = 128
x = np.linspace(1,10,res)
f = np.zeros(res)
courant_diff = np.zeros(res)
courant_adv = np.zeros(res)

dx = x[2]-x[1]
dt_g = 1e-2
Tmin = 0
Tmax = 10
alpha = 1.0

for i in range(res):
    f[i] = initial(x[i])
plt.plot(x,f)

ap = np.zeros(res)   #upper diagonal
a0 = np.zeros(res)   #diagonal
am = np.zeros(res)   #lower diagonal
d = np.zeros(res)    #right hand side
for t in range(2):
    for i in range (res):
        if (i==0):
            ap[i] = 0.0
            a0[i] = (1+2*1e2*dt_g/dx**2)
            am[i] = -1e2*dt_g/dx**2
            d[i] = (1-2*1e2*dt_g/dx**2)*f[i]+(1e2*dt_g/dx**2)*f[i+1]
        elif (i==res-1):
            ap[i] = -1e2*dt_g/dx**2
            a0[i] = (1+2*1e2*dt_g/dx**2)
            am[i] = 0.0
            d[i] = (1-2*1e2*dt_g/dx**2)*f[i]+(1e2*dt_g/dx**2)*f[i-1]
        else:
            ap[i] = -1e2*dt_g/dx**2
            a0[i] = (1+2*1e2*dt_g/dx**2)
            am[i] = -1e2*dt_g/dx**2
            d[i] = (1-2*1e2*dt_g/dx**2)*f[i]+(1e2*dt_g/dx**2)*f[i+1]+(1e2*dt_g/dx**2)*f[i-1]
    g = tridiag(res,ap,a0,am,d)

    for k in range(res):
            f[k] = g[k]
    plt.plot(x,g)

Any help regarding the algorithm or the code will be highly appreciated.

Answer 1

I think you might have mixed up some terminology. An exponential integrator would use some type of eigensolver or related approach to exactly calculate

$$f^{n + 1} = e^{\delta t\cdot \mathcal A}f^n$$

What you've done (in effect) is to use the Pade approximation

$$e^z \approx \frac{1 + \frac{z}{2}}{1 - \frac{z}{2}},$$

which is not an exponential integrator. (Your formula is missing some factors of 1/2.) Remember that for the diffusion operator, the eigenvalues all lie along the negative real axis. To understand why your time discretization isn't giving the results you want, you have to look at how this Pade approximation behaves along the negative real axis. The function $e^z$ goes to 0 as $z \to -\infty$, while the Pade approximation goes to -1. Can you see why this might produce oscillations?

The fancy term for this is that you've used a method that isn't L-stable. The first-order backward scheme $(I - \delta t\cdot A)f^{n + 1} = f^n$ is the simplest L-stable method.