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

Specific Domains Numerics
, ,
1

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:

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
If I include the final state of the solution it dwarfs everything else:
image

I solve the PDE using:

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)

2

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

2 Likes
3

I understand, but don’t think I’m getting all the way. I tried reading Upwind scheme - Wikipedia.

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

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

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

and

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

but the same thing happens with both of those.

4

You have a PDE:
image

and you discretize the advection term and the diffusion term.

Question: in your discretized model…
image

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

1 Like
5

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

1 Like
6

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.

1 Like
7

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!

3 Likes