# Discretize master equation

I'm trying to solve the following equation numerically using finite difference method.

$\frac{\partial P}{\partial t} = -\frac{\partial}{\partial x_1}[F_1(x_1,x_2)P] - \frac{\partial}{\partial x_2}[ F_2(x_1,x_2)P] + D(\frac{\partial^2P}{\partial x_1^2}+\frac{\partial^2P}{\partial x_2^2})$

where $P = P(x_1,x_2,t)$, $F_1 = \frac{\epsilon^2+x_1^2}{(1+x_1^2)(1+x_2)}-ax_1$, $F_2 = \frac{1}{\tau_0}(b-\frac{x_2}{1+cx_1^2})$.

The parameters $\epsilon, a,b,c, \tau_0$ are constants whose values can be 0.1, 0.1,0.1, 100, 5 respectively from Wang et al.. D can be a matrix, but for simplicity it is considered as constant, for example 0.001.

I tried the following approach but failed:

1. Discretize $F_1P$ and $F_2P$ using central difference at midpoint. For example, $\frac{\partial(F_1P)}{\partial x}|_k = \frac{1}{2}\bigg\{\frac{\partial(F_1P)}{\partial x}|_{k-1/2} + \frac{\partial(F_1P)}{\partial x}|_{k+1/2} \bigg\} =\frac{1}{2}\bigg\{\frac{-F_1(k-1,j)P(k-1,j)+F_1(k,j)P(k,j)}{dx} + \frac{-F_1(k,j)P(k,j)+F_1(k+1,j)P(k+1,j)}{dx} \bigg \}$
2. Discretize at midpoint $\frac{\partial^2P}{\partial x_1^2} = \frac{P(k+1,j)-2P(k,j)+P(k-1,j)}{dx^2}$

3. Neuman (relecting) boundary conditions, for example $\frac{\partial P}{\partial n}|_{y=0} = \frac{P(k,1)-P(k,-1)}{2dy} = 0 \Rightarrow P(k,1)= P(k,-1)$. This relationship is applied at boundaries along $x = y=0, x=y=L$

4. Initial condition $P_0$ is constructed via creating a random 2D matrix whose element $(k,j)$ in matrix corresponds to a probability $P_0(k,j)$. The sum of all elements is one, implying conservation of probability.

Can anyone help me how to discretize this problem? And how to check stability ? I would like to reproduce Fig 3 and 4 in Wang et al.

My Matlab code is shown for convenience:

for i = 1:Nt

Jx = F1.*P_previous;
Jy = F2.*P_previous;

for k = 2:Nx-1
for j = 2:Ny-1

%calculate partial of flux/ reaction

%Jx_partialx = (Jx(k+1,j) - Jx(k-1,j))/(2*dx);
%Jy_partialy = (Jy(k,j+1) - Jy(k,j-1))/(2*dy);
Jx_partialx = 1/2*((-F1(k-1,j)*P_previous(k-1,j)+F1(k,j)*P_previous(k,j))/dx + (-F1(k,j)*P_previous(k,j)+F1(k+1,j)*P_previous(k+1,j))/dx );
Jy_partialy = 1/2*((-F1(k,j-1)*P_previous(k,j-1)+F1(k,j)*P_previous(k,j))/dy + (-F1(k,j)*P_previous(k,j)+F1(k,j+1)*P_previous(k,j+1))/dy );

%calculate at each mesh point n for P_previous
%Fy*() --> partial derivative w.r.t y
%Fx*() --> partial derivative w.r.t x
P_partialx = P_previous(k+1,j) - 2*P_previous(k,j) + P_previous(k-1,j);
P_partialy = P_previous(k,j+1) - 2*P_previous(k,j) + P_previous(k,j-1);

P(k,j) = P_previous(k,j) - dt*(Jx_partialx + Jy_partialy) + Fy*P_partialy + Fx*P_partialx;    %Fx, Fy has dt
end
end

%P(1,:) = 0; P(Nx,:) = 0; P(:,1) = 0; P(:,Nx) = 0; %boundary condition
%reflecting boundary conditions
%Neumann boundary condition

%use partialP/partialy = (P(k,j+1) - P(k,j-1))/(2*dy)
%along y = 0
% --> P_previous(k,2) = P_previous(k,0)

j = 1;
jy = j+1;
for k = 2:Nx-1
Jy_partialy = 1/2*((-F1(k,jy)*P_previous(k,jy)+F1(k,j)*P_previous(k,j))/dy + (-F1(k,j)*P_previous(k,j)+F1(k,j+1)*P_previous(k,j+1))/dy );
Jx_partialx = 1/2*((-F1(k-1,j)*P_previous(k-1,j)+F1(k,j)*P_previous(k,j))/dx + (-F1(k,j)*P_previous(k,j)+F1(k+1,j)*P_previous(k+1,j))/dx );
P_partialx = P_previous(k+1,j) - 2*P_previous(k,j) + P_previous(k-1,j);
P_partialy = P_previous(k,j+1) - 2*P_previous(k,j) + P_previous(k,jy);
P(k,j) = P_previous(k,j) - dt*(Jx_partialx + Jy_partialy) + Fy*P_partialy + Fx*P_partialx;
end

%along y = C
% --> P_previous(k,Ny+1) = P_previous(k,Ny-1)
% --> Jy_partialy = 0
j = Ny;
jy = j-1;
for k = 2:Nx-1
Jy_partialy = 1/2*((-F1(k,j-1)*P_previous(k,j-1)+F1(k,j)*P_previous(k,j))/dy + (-F1(k,j)*P_previous(k,j)+F1(k,jy)*P_previous(k,jy))/dy );
Jx_partialx = 1/2*((-F1(k-1,j)*P_previous(k-1,j)+F1(k,j)*P_previous(k,j))/dx + (-F1(k,j)*P_previous(k,j)+F1(k+1,j)*P_previous(k+1,j))/dx );
P_partialx = P_previous(k+1,j) - 2*P_previous(k,j) + P_previous(k-1,j);
P_partialy = P_previous(k,jy) - 2*P_previous(k,j) + P_previous(k,j-1);
P(k,j) = P_previous(k,j) - dt*(Jx_partialx + Jy_partialy) + Fy*P_partialy + Fx*P_partialx;
end

%along x = 0
% --> P_previous(2,j) = P_previous(0,j)
% --> Jx_partialx = 0
k = 1;
kx = k+1;

for j = 2:Ny-1
Jy_partialy = 1/2*((-F1(k,j-1)*P_previous(k,j-1)+F1(k,j)*P_previous(k,j))/dy + (-F1(k,j)*P_previous(k,j)+F1(k,j+1)*P_previous(k,j+1))/dy );
Jx_partialx = 1/2*((-F1(kx,j)*P_previous(kx,j)+F1(k,j)*P_previous(k,j))/dx + (-F1(k,j)*P_previous(k,j)+F1(k+1,j)*P_previous(k+1,j))/dx );
P_partialx = P_previous(k+1,j) - 2*P_previous(k,j) + P_previous(kx,j);
P_partialy = P_previous(k,j+1) - 2*P_previous(k,j) + P_previous(k,j-1);
P(k,j) = P_previous(k,j) - dt*(Jx_partialx + Jy_partialy) + Fy*P_partialy + Fx*P_partialx;
end

%along x = C
% --> P_previous(Nx+1,j) = P_previous(Nx-1,j)
% --> Jx_partialx = 0
k = Nx;
kx = k-1;

for j = 2:Ny-1
Jy_partialy = 1/2*((-F1(k,j-1)*P_previous(k,j-1)+F1(k,j)*P_previous(k,j))/dy + (-F1(k,j)*P_previous(k,j)+F1(k,j+1)*P_previous(k,j+1))/dy );
Jx_partialx = 1/2*((-F1(k-1,j)*P_previous(k-1,j)+F1(k,j)*P_previous(k,j))/dx + (-F1(k,j)*P_previous(k,j)+F1(kx,j)*P_previous(kx,j))/dx );
P_partialx = P_previous(kx,j) - 2*P_previous(k,j) + P_previous(k-1,j);
P_partialy = P_previous(k,j+1) - 2*P_previous(k,j) + P_previous(k,j-1);
P(k,j) = P_previous(k,j) - dt*(Jx_partialx + Jy_partialy) + Fy*P_partialy + Fx*P_partialx;
end

%4 corners , partialP/partialx and partialP/partialy = 0
k = 1; j = 1;
kx = k+1; jy = j+1;
Jy_partialy = 1/2*((-F1(k,jy)*P_previous(k,jy)+F1(k,j)*P_previous(k,j))/dy + (-F1(k,j)*P_previous(k,j)+F1(k,j+1)*P_previous(k,j+1))/dy );
Jx_partialx = 1/2*((-F1(kx,j)*P_previous(kx,j)+F1(k,j)*P_previous(k,j))/dx + (-F1(k,j)*P_previous(k,j)+F1(k+1,j)*P_previous(k+1,j))/dx );

P_partialx = P_previous(k+1,j) - 2*P_previous(k,j) + P_previous(kx,j);
P_partialy = P_previous(k,j+1) - 2*P_previous(k,j) + P_previous(k,jy);
P(k,j) = P_previous(k,j) - dt*(Jx_partialx + Jy_partialy) + Fy*P_partialy + Fx*P_partialx;

k = 1; j = Ny;
kx = k+1; jy = j-1;
Jx_partialx = 1/2*((-F1(kx,j)*P_previous(kx,j)+F1(k,j)*P_previous(k,j))/dx + (-F1(k,j)*P_previous(k,j)+F1(k+1,j)*P_previous(k+1,j))/dx );
Jy_partialy = 1/2*((-F1(k,j-1)*P_previous(k,j-1)+F1(k,j)*P_previous(k,j))/dy + (-F1(k,j)*P_previous(k,j)+F1(k,jy)*P_previous(k,jy))/dy );
P_partialx = P_previous(k+1,j) - 2*P_previous(k,j) + P_previous(kx,j);
P_partialy = P_previous(k,jy) - 2*P_previous(k,j) + P_previous(k,j-1);
P(k,j) = P_previous(k,j) - dt*(Jx_partialx + Jy_partialy) + Fy*P_partialy + Fx*P_partialx;

k = Nx; j = 1;
kx = k-1; jy = j+1;
Jx_partialx = 1/2*((-F1(k-1,j)*P_previous(k-1,j)+F1(k,j)*P_previous(k,j))/dx + (-F1(k,j)*P_previous(k,j)+F1(kx,j)*P_previous(kx,j))/dx );
P_partialx = P_previous(kx,j) - 2*P_previous(k,j) + P_previous(k-1,j);
Jy_partialy = 1/2*((-F1(k,jy)*P_previous(k,jy)+F1(k,j)*P_previous(k,j))/dy + (-F1(k,j)*P_previous(k,j)+F1(k,j+1)*P_previous(k,j+1))/dy );
P_partialy = P_previous(k,j+1) - 2*P_previous(k,j) + P_previous(k,jy);
P(k,j) = P_previous(k,j) - dt*(Jx_partialx + Jy_partialy) + Fy*P_partialy + Fx*P_partialx;

k = Nx; j = Ny;
kx = k-1; jy = j-1;
Jx_partialx = 1/2*((-F1(k-1,j)*P_previous(k-1,j)+F1(k,j)*P_previous(k,j))/dx + (-F1(k,j)*P_previous(k,j)+F1(kx,j)*P_previous(kx,j))/dx );
Jy_partialy = 1/2*((-F1(k,j-1)*P_previous(k,j-1)+F1(k,j)*P_previous(k,j))/dy + (-F1(k,j)*P_previous(k,j)+F1(k,jy)*P_previous(k,jy))/dy );
P_partialx = P_previous(kx,j) - 2*P_previous(k,j) + P_previous(k-1,j);
P_partialy = P_previous(k,jy) - 2*P_previous(k,j) + P_previous(k,j-1);
P(k,j) = P_previous(k,j) - dt*(Jx_partialx + Jy_partialy)+ Fy*P_partialy + Fx*P_partialx;

%update P_previous
P_previous = P;
end


Thank you very much!

• It looks like you are using a forward Euler scheme? Compute the eigenvalues of the resulting matrix of your RHS and see if they are inside the region of absolute stability for forward Euler. – Kyle Mandli May 5 '18 at 15:36
• Yes, I'm using forward Euler scheme for P. The discretization method I used above produces negative values of probability distribution. I wonder whether this is related to stability issue. – canhochoi May 7 '18 at 5:21
• It could be, in general using a forward euler time stepping scheme with a diffusion term is ill-advised due to the stability restriction of the time step. There are multiple ways to proceed but I would try using a backward Euler method first to see if that works better for you. There are also fancier methods such as RK SSP methods that might help. – Kyle Mandli May 7 '18 at 15:46
• I've noticed ADI (Alternating direction implicit) can be used to solve 2D diffusion with high stability. With no advection/reaction terms, I've correctly produced time evolution of probability of random walk by using sparse matrix solving. However, when advection, reaction are included, it can't produce any result (I guess because some region of probability become negative). Do you know how to code the method you suggested? I'm familiar with matlab and a bit of Python. – canhochoi May 11 '18 at 8:46