Problem in solving a Fokker-Planck equation

40 次查看（过去 30 天）
Sagar Karmakar
Sagar Karmakar 2023-7-21
编辑： Torsten 2023-7-21
Ran in:
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 个评论
Torsten
Torsten 2023-7-21
编辑：Torsten 2023-7-21
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 2023-7-21
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2023-7-21
编辑：Torsten 2023-7-21
Ran in:
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 个评论
Sagar Karmakar
Sagar Karmakar 2023-7-21
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 2023-7-21
编辑：Torsten 2023-7-21
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Boundary Conditions 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by