Problem in solving a Fokker-Planck equation
24 views (last 30 days)
Show older comments
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
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 ?
Accepted Answer
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(:))
surf(x,y,P(:,:,Nt))
2 Comments
More Answers (0)
See Also
Categories
Find more on Boundary Conditions in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!