How to solve a Fokker Planck equation in both space and time using finite difference method

9 Ansichten (letzte 30 Tage)
I want the solution of a Fokker-Planck equation using Finite Difference Scheme to find the probability distribution P and steady state distribution P(ss)(which is the value of P for very large time t) .The numerical scheme is as given in the following image
Where F1 and F2 are given as
I derive the first term of the RHS of (5) by finding partial derivative of F1 and F2 and the matlab code become as shown in below but I did not get the probability distribution P(x,y,t) in x-y plane where colorbar will denote the density of P. I don't know my code is right or not as I am a beginner, So I need a help from someone who can verify this code and tell me where I am wrong
clc;
clear all;
Nt=1000;Nx=30;Ny=30;
Tmax=1;Lx=3;Ly=3;
dt=Tmax/Nt;
dx=Lx/Nx;
dy=Ly/Ny;
t=(0:Nt)*dt;
x=(0:Nx)*dx;
y=(0:Ny)*dy;
% P=zeros(Nx,Ny,Nt);
a=1;b1=1;b2=1;k1=1;k2=1;S=0.5;n=4;D=0.04; %a=a1=a2
% Boundary conditons
for k=1:Nt+1
for i=1:Nx+1
P(i,1,k)=0;
P(i,Ny+1,k)=0;
end
for j=1:Ny+1
P(1,j,k)=0;
P(Nx+1,j,k)=0;
end
end
% ICs
for i=1:Nx+1
for j=1:Ny+1
P(i,j,1)=40;
end
end
% P(:,:,1)=ones(Nx+1,1)*exp(-(x.^2+y.^2)/2);
%Compute P
for k=1:Nt
for j=2:Ny
for i=2:Nx
P(i,j,k+1)=P(i,j,k)+dt*...
(P(i,j,k)*(((n*a*S^n*x(i)^(n-1))./((S^n)+x(i)^n)^2)+((n*a*(S^n)*y(j)^(n-1))./((S^n)+y(j)^n)^2)-(k1+k2)))-(((a*x(i)^n)./(S^n+x(i)^n))+((b1*S^n)./(S^n+y(j)^n))-k1*x(i))*((P(i+1,j,k)-P(i-1,j,k))./(2*dx))...
-(((a*y(j)^n)./(S^n+y(j)^n))+((b2*S^n)./(S^n+x(i)^n))-k2*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
%plot(t,P(x,y,:));
imagesc(x,y,P(:,:,500))
colormap jet;
colorbar;
% pcolor(P); shading interp;
% colorbar;
% colormap jet;
  2 Kommentare
Sagar Karmakar
Sagar Karmakar am 30 Apr. 2023
I want to find an image like this
But it is not possible for me to find such image from the above code. I want to know if I have to change anything in the code?

Melden Sie sich an, um zu kommentieren.

Antworten (0)

Community Treasure Hunt

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

Start Hunting!

Translated by