Unstablities when solving initial value problem with Julia DifferentialEquations.jl

I’m trying to numerically determine the stationary solution of Fokker-Plank equations \frac{\partial P}{\partial t}(x,t) = -\frac{\partial }{\partial x}\left[A(x,t) P(x,t)\right]+\frac{1}{2}\frac{\partial^2}{\partial^2 x}\left[B(x,t) P(x,t)\right] using DifferentialEquations.jl of Julia.

The idea is to evolve P with t until it converges/barely changes with t. For my question I have prepared a simpler model - finding the stationary solution of
\frac{\partial P}{\partial t}(x,t) = -\frac{\partial P}{\partial x}(x,t) + x P(x,t).
Analytically, the equation exhibits the stationary solution P_{st}=C \exp(-x^2/2).

To do this numerically I discretized the problem on an x-grid and employed periodic boundary conditions. Here is my code:

using DifferentialEquations
using Plots

function gradient(f, dx)
    ret = (circshift(f, -1) .- circshift(f, 1)) ./ (2 * dx)
    ret[1] = (f[end] - f[1]) / dx
    return ret

function fokker_planck!(du, u, p, t)
    D, x, dx = p
    du .= gradient(u, dx) + x .* u

N = 1000
x_min = -5.0
x_max = 5.0
dx = (x_max - x_min) / N
x = x_min:dx:x_max-dx
D = 1

# IC
initial_condition = exp.(-x .^ 2 ./ (2 * D))
# ODEProblem
prob = ODEProblem(fokker_planck!, initial_condition, (0.0, 4), [D, x, dx])
# solve
solution = solve(prob, Tsit5(), dt=1e-6);

I wanted to test this code by first initializing it with the stationary solution. So I expected that during the evolution P stays unchanged. However, I experienced that the solution is very unstable and when I initialize it differently it does not converge to the stationary solution. Here is an image of my solution:

Unstable solution

I have tested different solvers, stepsize and discretzations of x (larger intervals and finer grids) but so far I have not managed to obtain a stable solution. I suspected this is due to the error of the finite differences approximating the spatial derivative in the update but finer x-grid does not help. Is there an obvious mistake in my approach or are there other ways to approach this problem using julia?

The central difference discretization is not a stable discretization of the advection equation. This is a well-known phonomena:

Use a stable discretization like upwinding and you should be fine.

1 Like

@ChrisRackauckas Thank you! Is there a Julia library that implements an upwind scheme? Especially, if the coefficients are non-constant like for Fokker-Planck equations.

I have implemented a solution using MethodOfLines.jl but the compile times grow too large for the system size I actually want to solve.

MethodOfLines.jl does, and we’re working on the compile times with JuliaSimCompiler backends. FiniteVolumeMethod.jl might be a good solution for now.