Problem in solving a Fokker-Planck equation

I have written a matlab code to solve a Fokker-Planck equation using finite difference method. The PDE is of the form

Where F1(x1,x2) = ((a*x1^n)/(S^n+x1^n) + (b*S^n)/(S^n+x2^n) - k1*x1)

F2(x1,x2) = ((a*x2^n)/(S^n+x2^n)+ (b*S^n)/(S^n+x1^n) - k1*x2)

and the finite difference scheme is given as

When the parameter a=1, there should be three peak for P(: , :, Nt) near three different spatial point and for a=0.2 there will be two peak (in surface plot)

But for my matlab code I am getting only one peak for any value of a.

I can't understand what going wrong....

clc;
clear all
Nx=100;Ny=100;Nt=1000;
xmin=0;xmax=3;
ymin=0;ymax=3;
dx=(xmax-xmin)/(Nx-1);dy=-(ymin-ymax)/(Ny-1);
tmin=0;tmax=1;
dt=-(tmin-tmax)/(Nt-1);
D=0.04;
%define grids
x=(0:Nx-1)*dx;
y=(0:Ny-1)*dy;
k=(0:Nt-1)*dt;
[X,Y]=meshgrid(x,y);
%parameter values
a=0.2;b=1;S=0.5;k1=1;k2=1;n=4;

% F1=@(u,v)((a*u^n)/(S^n+u^n)+(b*S^n)/(S^n+v^n)-k1*u);
% F2=@(u,v)((a*v^n)/(S^n+v^n)+(b*S^n)/(S^n+u^n)-k1*v);
F1 =((a*X.^n)./(S^n+X.^n)+(b*S^n)./(S^n+Y.^n)-k1*X);
F2 =((a*Y.^n)./(S^n+Y.^n)+(b*S^n)./(S^n+X.^n)-k1*Y);

P = zeros(Nx,Ny,Nt);

%Initial conditions
P(:,:,1)=exp(-((X-0.1).^2 / (2)) - ((Y).^2 / (2)));
%% Boundary conditions
P(1,:,:)=0;P(Nx,:,:)=0;
P(:,1,:)=0;P(:,Ny,:)=0;


for i=2:Nx-1
    for j=2:Ny-1
        for k=1:Nt-1
            P(i,j,k+1)=P(i,j,k)+dt*(-P(i,j,k)*((F1(i+1,j)-F1(i-1,j))/(2*dx))-P(i,j,k)*((F2(i,j+1)-F2(i,j-1))/(2*dy))-(F1(i,j))*...
                 ((P(i+1,j,k)-P(i-1,j,k))/(2*dx))-F2(i,j)*((P(i,j+1,k)-P(i,j-1,k))/(2*dx))+D*(((P(i+1,j,k)-2*P(i,j,k)+P(i-1,j,k))/(dx^2))+...
                 ((P(i,j+1,k)-2*P(i,j,k)+P(i,j-1,k))/(dy^2))));

        end
    end
end

 surf(X,Y,P(:,:,Nt-1));

Please help me to solve it correctly.

Sorry for my bad english.

[Torsten]

The boundary conditions and the initial conditions are not allowed to come into conflict for your equation. This seems to be the case for the example in Zhou et al. because it's written there: "Here, the Neumann or Dirichlet boundary conditions hav no influence upon the result." It is not the case for your initial function exp(-((X-0.1).^2 / (2)) - ((Y).^2 / (2))). Plot it, and you will see why.

What profile do you expect after t = 1 ?

[Sagar Karmakar]

At the final time expect to have two peak for the probability distribution when a =0.2 and three peak when a=1.

If the problem lies on the consideration of initial and boundary condiyions. Then what type of initial condition should take so that it allow my given BCs. As the P is a PDF , so is it good to have a delta function initially at t=0?

But, I don't know how to write a delta function in matlab.

Answer 1

[Torsten]

The order of your loops is wrong. Try

Nx=100;Ny=100;Nt=1000;
xmin=0;xmax=3;
ymin=0;ymax=3;
dx=(xmax-xmin)/(Nx-1);
dy=(ymax-ymin)/(Ny-1);
tmin=0;tmax=1;
dt=(tmax-tmin)/(Nt-1);
D=0.04;
%define grids
x=(0:Nx-1)*dx;
y=(0:Ny-1)*dy;

%parameter values
a=1;b=1;S=0.5;k1=1;k2=1;n=4;
F1 =@(x,y)(a*x.^n)./(S^n+x.^n)+(b*S^n)./(S^n+y.^n)-k1*x;
F2 =@(x,y)(a*y.^n)./(S^n+y.^n)+(b*S^n)./(S^n+x.^n)-k1*y;
P = zeros(Nx,Ny,Nt);
%Initial conditions
for i = 1:Nx
  for j = 1:Ny
     P(i,j,1)=exp(-((x(i)-0.1)^2 + y(j)^2)/2);
  end
end
%% Boundary conditions
P(1,:,:)=0;P(Nx,:,:)=0;
P(:,1,:)=0;P(:,Ny,:)=0;
for k=1:Nt-1
  for i=2:Nx-1
    for j=2:Ny-1
            P(i,j,k+1)=P(i,j,k)+dt*(-P(i,j,k)*((F1(x(i+1),y(j))-F1(x(i-1),y(j)))/(2*dx))...
                                    -P(i,j,k)*((F2(x(i),y(j+1))-F2(x(i),y(j-1)))/(2*dy))...
                                    -F1(x(i),y(j))*(P(i+1,j,k)-P(i-1,j,k))/(2*dx)...
                                    -F2(x(i),y(j))*(P(i,j+1,k)-P(i,j-1,k))/(2*dy)...
                                    +D*((P(i+1,j,k)-2*P(i,j,k)+P(i-1,j,k))/dx^2+...
                                       ((P(i,j+1,k)-2*P(i,j,k)+P(i,j-1,k))/dy^2)));
        end
    end
end
PNt = P(:,:,Nt);
m = max(PNt(:))

m = 1.2427

surf(x,y,P(:,:,Nt))

[Sagar Karmakar]

Hi! Thank you very much!

But one more question

If I want a delta function as initial condition, how can I express it in matlab code?

[Torsten]

If I want a delta function as initial condition, how can I express it in matlab code?

It doesn't work numerically. Try to use a Gaussian (as you did) with small variance and a fine mesh around it.