Clear Filters
Clear Filters

Problem in solving a Fokker-Planck equation

24 views (last 30 days)
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.
  2 Comments
Torsten
Torsten on 21 Jul 2023
Edited: Torsten on 21 Jul 2023
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
Sagar Karmakar on 21 Jul 2023
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.

Sign in to comment.

Accepted Answer

Torsten
Torsten on 21 Jul 2023
Edited: Torsten on 21 Jul 2023
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))
  2 Comments
Sagar Karmakar
Sagar Karmakar on 21 Jul 2023
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
Torsten on 21 Jul 2023
Edited: Torsten on 21 Jul 2023
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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!