# Solving PDE Diffusion (using ODE solvers). Total mass increases towards infinity (should stay constant)

**URL:** <https://discourse.julialang.org/t/solving-pde-diffusion-using-ode-solvers-total-mass-increases-towards-infinity-should-stay-constant/67492>\
**Category:** Numerics\
**Tags:** diffeq, pde, ode\
**Created:** [September 1, 2021, 11:01am UTC](https://discourse.julialang.org/t/solving-pde-diffusion-using-ode-solvers-total-mass-increases-towards-infinity-should-stay-constant/67492 "2021-09-01T11:01:57Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [September 1, 2021, 11:01am UTC](https://discourse.julialang.org/t/solving-pde-diffusion-using-ode-solvers-total-mass-increases-towards-infinity-should-stay-constant/67492/1 "2021-09-01T11:01:57Z")

</div>

I am trying to solve a very simple case of the chemical Fokker-Planch equations. It is basically a diffusion equation of the probability density function of the state of a chemical reaction network. The process is very simple (a single molecule, X, is produced at constant rate λ and degraded linearly at rate β).

I am just trying to solve it once, and create a plot of the time development. The process is very simple, so performance is not an issue (as long as it doesn’t take 2 months to solve…). However, I do want a correct(ish) solution.

The diffusion is in 1d space, and the equation looks like this:  
∂P/∂t = - (∂/∂x)⋅[(λ-βx)⋅P] + (1/2)⋅(∂²/∂x²)⋅[(λ+βx)⋅P]  
I tried to make a discrete version:

```julia
grid = 0:0.05:200.0; Δx = grid[2]
A(i,u,λ,β,grid) = (0 < i <= length(grid)) ? u[i]*(λ-β*grid[i]) : 0.0
B(i,u,λ,β,grid) = (0 < i <= length(grid)) ? u[i]*(λ+β*grid[i]) : 0.0

function cfpe(du,u,p,t)
    λ,β = p
    for i = 1:length(grid)
        du[i] = u[i] - (0.5/Δx) * (A(i+1,u,λ,β,grid) - A(i-1,u,λ,β,grid)) + (0.5/Δx^2) * (B(i+1,u,λ,β,grid) -2*B(i,u,λ,β,grid) + B(i-1,u,λ,β,grid))
    end
end

```

using  
(∂/∂x)f ≈ (1/2) (1/Δx) [F(i+1)-F(i-1)]  
(∂/∂x²)f ≈ (1/2) (1/Δx²) [F(i+1)-2⋅F(i)+F(i-1)]

However, my solving doesn’t seem to work at all. For a starter, `sum(sol[i])`-\>Inf as `i` grows. This is plots of the solution for various times:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/f/d/fd5e9af6142e2ac3b2a82dd702e8110cb874c913.png)  
If I include the final state of the solution it dwarfs everything else:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/3/d36f193a09762e9d0debfc02024c831941ddb101.png)

I solve the PDE using:

```julia
using DifferentialEquations
u0 = zeros(length(grid)); u0[findfirst(grid.>=10.0)] = 1;
parameters = [10.,0.1]; tend = 10.0;
prob = ODEProblem(cfpe,u0,(0.,tend),parameters);
@time sol = solve(prob,tstops=0:0.1:tend);

```

(The functions are naturally =0 at 0, and I simply assume that they are 0 for very large x (x=200). The solution grows significantly before there’s a significant dnesity at these values, so I don’t think the boundary conditions are a problem)

(The initial condition is everything gathered in a single point. Initially, I thought the discontinuous initial condition might pose a problem, but the solution quickly becomes continuous, and the problem still persists)

---

<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:** [September 1, 2021, 12:16pm UTC](https://discourse.julialang.org/t/solving-pde-diffusion-using-ode-solvers-total-mass-increases-towards-infinity-should-stay-constant/67492/2 "2021-09-01T12:16:41Z")

</div>

> [@Torkel](#):
>
> ```julia
> du[i] = u[i] - (0.5/Δx) * (A(i+1,u,λ,β,grid) - A(i-1,u,λ,β,grid)) + (0.5/Δx^2) * (B(i+1,u,λ,β,grid) -2*B(i,u,λ,β,grid) + B(i-1,u,λ,β,grid))
> 
> ```

You should not use central difference for the advection term. It can give an unconditionally unstable discretization. Use upwinding.

---

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [September 1, 2021, 1:00pm UTC](https://discourse.julialang.org/t/solving-pde-diffusion-using-ode-solvers-total-mass-increases-towards-infinity-should-stay-constant/67492/3 "2021-09-01T13:00:32Z")

</div>

I understand, but don’t think I’m getting all the way. I tried reading [Upwind scheme - Wikipedia](https://en.wikipedia.org/wiki/Upwind_scheme).

Initially, I tried guessing their _a_ to be _(λ-βx)_, and so used:

```julia
advection = (1/Δx) * ( (λ-β*grid[i]) > 0 ? A(i+1,u,λ,β,grid) - A(i,u,λ,β,grid) : A(i,u,λ,β,grid) - A(i-1,u,λ,β,grid))
du[i] = u[i] - advection + (0.5/Δx^2) * (B(i+1,u,λ,β,grid) -2*B(i,u,λ,β,grid) + B(i-1,u,λ,β,grid))

```

but with the same phenomena happening. I also tried both

```julia
advection = (1/Δx) * (A(i+1,u,λ,β,grid) - A(i,u,λ,β,grid))

```

and

```julia
advection = (1/Δx) * (A(i,u,λ,β,grid) - A(i-1,u,λ,β,grid))

```

but the same thing happens with both of those.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [September 1, 2021, 2:34pm UTC](https://discourse.julialang.org/t/solving-pde-diffusion-using-ode-solvers-total-mass-increases-towards-infinity-should-stay-constant/67492/4 "2021-09-01T14:34:03Z")

</div>

You have a PDE:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/3/1359636bbdad955033c22f85526b42f6e8e764e6.png)

and you discretize the advection term and the diffusion term.

_Question_: in your discretized model…  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/6/7/67e7fc38ea1099664b3110d9425afa23736d8dba.png)

where does the first RHS term `u[i]` come from in the original PDE?

---

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [September 1, 2021, 3:00pm UTC](https://discourse.julialang.org/t/solving-pde-diffusion-using-ode-solvers-total-mass-increases-towards-infinity-should-stay-constant/67492/5 "2021-09-01T15:00:50Z")

</div>

You are right, I clearly messed that one up (I wrote to `du` the new value of `u` instead of `du`). Thanks 🙂

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [September 1, 2021, 3:32pm UTC](https://discourse.julialang.org/t/solving-pde-diffusion-using-ode-solvers-total-mass-increases-towards-infinity-should-stay-constant/67492/6 "2021-09-01T15:32:24Z")

</div>

Hope removing `u[i]` helped stabilize the simulation. This term gives an eigenvalue at +1, i.e., highly unstable, and removing it should definitely help.

---

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [September 1, 2021, 3:59pm UTC](https://discourse.julialang.org/t/solving-pde-diffusion-using-ode-solvers-total-mass-increases-towards-infinity-should-stay-constant/67492/7 "2021-09-01T15:59:45Z")

</div>

Yes, now it works well. I didn’t have that in the beginning, but then things didn’t work, and then I got this idea that it had to be there and that maybe was the problem, but I never realised that it was really stupid. Probably the initial error was solved by Chris’ suggestion, and then now when I removed the `u[i]` things works as intended!