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, 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(τ)
# 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,:])
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.
.