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

14 ビュー (過去 30 日間)
Sagar Karmakar
Sagar Karmakar 2023 年 4 月 30 日
編集済み: Sagar Karmakar 2023 年 5 月 2 日
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 件のコメント
John D'Errico
John D'Errico 2023 年 4 月 30 日
And what is your question? Why are you posting this?
Sagar Karmakar
Sagar Karmakar 2023 年 4 月 30 日
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?

サインインしてコメントする。

回答 (0 件)

カテゴリ

Help Center および File ExchangeBoundary Conditions についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by