How to solve a Fokker Planck equation in both space and time using finite difference method
9 views (last 30 days)
Show older comments
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 Comments
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!