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
end
function fokker_planck!(du, u, p, t)
D, x, dx = p
du .= gradient(u, dx) + x .* u
end
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:
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?