# Unstablities when solving initial value problem with Julia DifferentialEquations.jl

**URL:** <https://discourse.julialang.org/t/unstablities-when-solving-initial-value-problem-with-julia-differentialequations-jl/110910>\
**Category:** Modelling & Simulations\
**Tags:** differentialequation\
**Created:** [February 28, 2024, 4:37pm UTC](https://discourse.julialang.org/t/unstablities-when-solving-initial-value-problem-with-julia-differentialequations-jl/110910 "2024-02-28T16:37:04Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![paho](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paho/32/207349_2.png) [@paho](https://discourse.julialang.org/u/paho)\
**Post date:** [February 28, 2024, 4:37pm UTC](https://discourse.julialang.org/t/unstablities-when-solving-initial-value-problem-with-julia-differentialequations-jl/110910/1 "2024-02-28T16:37:04Z")

</div>

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:

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

[![Unstable solution](https://global.discourse-cdn.com/julialang/original/3X/e/b/eb2b9ae5ec2f504bbef54be00e8bdf7076c135e1.png)](https://i.stack.imgur.com/CzWFf.png)

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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 28, 2024, 7:56pm UTC](https://discourse.julialang.org/t/unstablities-when-solving-initial-value-problem-with-julia-differentialequations-jl/110910/2 "2024-02-28T19:56:46Z")

</div>

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

> <https://scicomp.stackexchange.com/questions/27737/advection-equation-with-finite-difference-importance-of-forward-backward-or-ce>

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

---

<div class="post-metadata">

**Author:** ![paho](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paho/32/207349_2.png) [@paho](https://discourse.julialang.org/u/paho)\
**Post date:** [February 28, 2024, 9:38pm UTC](https://discourse.julialang.org/t/unstablities-when-solving-initial-value-problem-with-julia-differentialequations-jl/110910/3 "2024-02-28T21:38:43Z")

</div>

@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](https://docs.sciml.ai/MethodOfLines/stable/) but the compile times grow too large for the system size I actually want to solve.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 18, 2024, 12:05am UTC](https://discourse.julialang.org/t/unstablities-when-solving-initial-value-problem-with-julia-differentialequations-jl/110910/4 "2024-03-18T00:05:23Z")

</div>

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