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

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)

dx = x-x
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.

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.

• In the context of large sparse systems, Pade approximations are used in lieu of direct computation of the matrix exponential via eigendecomposition because the application of Pade approximants on a vector can be computed in a matrix-free manner. You are correct though that the deviation between the approximant and the exponential can causes stability issues though. May 14 at 1:20
• I thought that it's only an exponential integrator if you're really really actually no joking trying to do the matrix exponential -- an intermediate step of which might be to use Pade approximants, probably followed by some Krylov magic. This is expensive but has a big payoff. If you're just using regular old Pade approximants then it's no different from some kind of Runge-Kutta method right? I've never actually used exponential integrators, just going off of what I've read here. May 14 at 1:27
• It's been a while since I've read about this but I think that Pade + Krylov is actually one of the more stable ways of doing exponential integration for large problems because it requires you to evaluate the function $\phi(z) = e^{z}/z$ on the operator, which can be unstable. Even with Pade approximation, the stage coefficients of the ERK method still depend on the linear part of the evolution operator, whereas in standard RK they are fixed, so they aren't generally equivalent I'm pretty sure. May 14 at 1:40
• *should be $(e^z-1)/z$ May 14 at 2:12
• @DanielShapero thank you for the clear and understandable answer. Although I was wondering if the approximated exponentiation of the matrix in the question is not L stable is there other approximation which is obey this stability. A search lead me to an approximation which is $e^{z}=\frac{1+(1-c)z}{1-cz+(c-1/2)z^{2}}$ which goes to 0 as $z$ goes to $-\infty$. I was actually trying to implement exponential integrator for linear diffusion and found out the approximation mentioned in the question in stackexchange only. May 14 at 5:24