# Stability issue with MethodOfLines

**URL:** <https://discourse.julialang.org/t/stability-issue-with-methodoflines/84133>\
**Category:** Modelling & Simulations\
**Created:** [July 12, 2022, 11:07pm UTC](https://discourse.julialang.org/t/stability-issue-with-methodoflines/84133 "2022-07-12T23:07:10Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![marcofrancis](https://avatars.discourse-cdn.com/v4/letter/m/bcef8e/32.png) [@marcofrancis](https://discourse.julialang.org/u/marcofrancis)\
**Post date:** [July 12, 2022, 11:07pm UTC](https://discourse.julialang.org/t/stability-issue-with-methodoflines/84133/1 "2022-07-12T23:07:10Z")

</div>

Hi!  
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}  
x\in[-2,2]  
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](https://diffeq.sciml.ai/stable/solvers/ode_solve/), 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

 ![Screenshot 2022-07-13 003146](https://global.discourse-cdn.com/julialang/original/3X/9/f/9fbbe0564c5a878749e7b700fce6183df5eb7115.png)

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

```julia
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(τ)

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

plot(xs,g_lines[20,:])

```

---

<div class="post-metadata">

**Author:** ![JackDevine](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jackdevine/32/1048_2.png) [@JackDevine](https://discourse.julialang.org/u/JackDevine)\
**Post date:** [July 13, 2022, 12:57am UTC](https://discourse.julialang.org/t/stability-issue-with-methodoflines/84133/2 "2022-07-13T00:57:50Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![marcofrancis](https://avatars.discourse-cdn.com/v4/letter/m/bcef8e/32.png) [@marcofrancis](https://discourse.julialang.org/u/marcofrancis)\
**Post date:** [July 13, 2022, 7:26am UTC](https://discourse.julialang.org/t/stability-issue-with-methodoflines/84133/3 "2022-07-13T07:26:49Z")

</div>

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [July 13, 2022, 9:37am UTC](https://discourse.julialang.org/t/stability-issue-with-methodoflines/84133/4 "2022-07-13T09:37:16Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![skleinbo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skleinbo/32/36080_2.png) [@skleinbo](https://discourse.julialang.org/u/skleinbo)\
**Post date:** [July 13, 2022, 9:59am UTC](https://discourse.julialang.org/t/stability-issue-with-methodoflines/84133/5 "2022-07-13T09:59:51Z")

</div>

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

```julia
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)

![image](https://global.discourse-cdn.com/julialang/original/3X/5/d/5d7edfbdcd3bf8c204f58fac2850ad24af71d072.png)

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

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

```

Note that these BCs alone do not resolve the instability.

---

<div class="post-metadata">

**Author:** ![marcofrancis](https://avatars.discourse-cdn.com/v4/letter/m/bcef8e/32.png) [@marcofrancis](https://discourse.julialang.org/u/marcofrancis)\
**Post date:** [July 13, 2022, 10:14am UTC](https://discourse.julialang.org/t/stability-issue-with-methodoflines/84133/6 "2022-07-13T10:14:38Z")

</div>

OMG! You made my day 😃 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](https://discourse.julialang.org/t/dirac-delta-in-boundary-conditions-for-methodoflines-jl/80567/5) but somehow even putting the time derivative on the other side does not help

---

<div class="post-metadata">

**Author:** ![skleinbo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skleinbo/32/36080_2.png) [@skleinbo](https://discourse.julialang.org/u/skleinbo)\
**Post date:** [July 13, 2022, 12:18pm UTC](https://discourse.julialang.org/t/stability-issue-with-methodoflines/84133/7 "2022-07-13T12:18:04Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [July 13, 2022, 5:08pm UTC](https://discourse.julialang.org/t/stability-issue-with-methodoflines/84133/8 "2022-07-13T17:08:02Z")

</div>

Thank you for your interest in MethodOfLines.jl 🙂 .

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:

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

```

So this problem remaining would be very strange indeed.

---

<div class="post-metadata">

**Author:** ![marcofrancis](https://avatars.discourse-cdn.com/v4/letter/m/bcef8e/32.png) [@marcofrancis](https://discourse.julialang.org/u/marcofrancis)\
**Post date:** [July 13, 2022, 6:19pm UTC](https://discourse.julialang.org/t/stability-issue-with-methodoflines/84133/9 "2022-07-13T18:19:41Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [July 13, 2022, 6:33pm UTC](https://discourse.julialang.org/t/stability-issue-with-methodoflines/84133/10 "2022-07-13T18:33:31Z")

</div>

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