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]
            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])

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]
            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]

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


1 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.

  • $\begingroup$ 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. $\endgroup$
    – whpowell96
    Commented May 14, 2023 at 1:20
  • 1
    $\begingroup$ 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. $\endgroup$ Commented May 14, 2023 at 1:27
  • $\begingroup$ 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. $\endgroup$
    – whpowell96
    Commented May 14, 2023 at 1:40
  • $\begingroup$ *should be $(e^z-1)/z$ $\endgroup$
    – whpowell96
    Commented May 14, 2023 at 2:12
  • $\begingroup$ @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. $\endgroup$
    – Sayan
    Commented May 14, 2023 at 5:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.