Stability issue with MethodOfLines

I’m trying to solve a (I would think) simple PDE with MethodOfLines.jl
The PDE in question is:

\partial_\tau g(τ,x) = \frac{1}{2}σ^2\partial^2_{x} g(τ,x)+ λ(\bar{x} - x)*\partial_x g(τ,x) - r*g(τ,x)
g(0,x) = e^{-kx}
where \sigma,\lambda,\bar{x},r,k are parameters.
If one wanted an analytical solution I’m pretty sure that by using Feynman-Kac and the formula for the expectation of a lognormal variable one can get it, but this is not why I’m asking your help.

In order to solve this problem numerically, I’m imposing some Neumann condition
\partial_x g(\tau,0)=\partial_x g(\tau,2)=0

If I solve this problem (I tried many solvers from ODE Solvers · DifferentialEquations.jl, like Tsit5(), Sundials.CVODE_BDF(), TRBDF2(), etc. ) on \tau\in [0,1]
I get some wigglines in the solution that then makes the whole solution unstable.
In the pic you can see an example where I plotted the cross section of the solution for \tau = 0.12

What am I doing wrong? Is it a problem of using an adaptive quadrature method?
Here a minimal code that reproduces the issue:

using MethodOfLines, OrdinaryDiffEq, Sundials, ModelingToolkit, Plots
import ModelingToolkit: Interval

τ_min= 0.
τ_max = 1
τ_span = (τ_min,τ_max)

x_min = -2.
x_max = 2.

dx = 0.01

σ = 0.03
λ = 0.1
xb = 1.
r = 0.01
k = 1.
h(x) = exp.(-k*x)

@parameters τ x
@variables g(..)

Dx = Differential(x)
Dxx = Differential(x)^2
Dτ = Differential(τ)

eq  = 0 ~ 0.5*σ^2*Dxx(g(τ,x)) + λ*(xb - x)*Dx(g(τ,x)) - Dτ(g(τ,x)) - r*g(τ,x) 

# Initial and boundary conditions
bcs = [g(τ_min,x) ~ h(x),
       Dx(g(τ,x_max)) ~ 0,
       Dx(g(τ,x_min)) ~ 0]

# Space and time domains
domains = [τ ∈ Interval(τ_min,τ_max),
           x ∈ Interval(x_min,x_max)]

@named pde_system = PDESystem(eq,bcs,domains,[τ,x],[g(τ, x)])

discretization = MOLFiniteDifference([x=>dx], τ)
prob = discretize(pde_system,discretization)
sol = MethodOfLines.solve(prob, Sundials.CVODE_BDF())
grid = get_discrete(pde_system, discretization)
xs = grid[x]
τs = sol.t
# Retrieve shaped solution
g_lines = [map(d -> sol[d][i], grid[g(τ, x)]) for i in 1:length(sol.t)]


That does look like quite a nasty instability.

It seems to me that your initial and boundary conditions don’t match unless you suppose that there is a discontinuity in the first derivative at the left hand side (where the instability is occuring). Have you tried the same problem with initial conditions that match the boundary conditions?


Mmm, something curious is going on here. I tried what @JackDevine suggested using for example:

h(x) = (x>= x_min+5*dx).*exp.(-k*x)+(1-(x>= x_min+5*dx)).*exp.(-k*(x_min+5*dx))

which has a kink, but satisfies the incriminated boundary condition on the left limit of the domain. Nothing really changed.
On the other hand, if I try to change the value of xb, i.e. where the coefficient of the first derivative changes sign to 0 from 1 (closer to the left boundary) then everything runs smoothly.

Can you try a smoother version ? MOL is basically finite differences so I expect big gradients if your initial condition is not smooth.

I believe MOL uses the wrong sign when discretizing the first derivate. It should generate the appropriate upwind operator, but gets confused by (xb-x) I reckon.

When I turn the sign around

eq  = Dτ(g(τ,x)) ~ 0.5*σ^2*Dxx(g(τ,x)) - λ*(x - xb)*Dx(g(τ,x)) - r*g(τ,x)

the instability is gone (I’ve extended the integration time to 5)


I have also used initial conditions that are commensurable with the boundary conditions.

h(x) = 1.0 + cos(π*(x-x_min)/(x_max-x_min))

Note that these BCs alone do not resolve the instability.


OMG! You made my day :smiley: Yeah I can confirm that even starting from a kosher initial condition (e.g. very peaked Gaussian far from the boundaries) will result in instability UNLESS the magic sign turning is performed.

I thought it could be related to Dirac delta in boundary conditions for MethodOfLines.jl - #5 by xtalax but somehow even putting the time derivative on the other side does not help

1 Like

Great! It might still be worth opening an issue with MOL. @xtalax

1 Like

Thank you for your interest in MethodOfLines.jl :slight_smile: .

This was a known issue in v0.3.0, but was fixed (I thought!) in v0.3.1. Which version of MOL are you using?

The first thing MOL does (in v0.3.1) is cast equations to a cardinal form:

pdeeqs = [eq.lhs - eq.rhs ~ 0 for eq in pdesys.eqs]

So this problem remaining would be very strange indeed.

1 Like

I have 0.3.1. As I said somehow even bringing everything on the lhs of the equation still makes it unstable.

Ah I see, there’s some weirdness going on here, Thank you for finding a valuable edge case!