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

10 次查看（过去 30 天）
Sagar Karmakar
Sagar Karmakar 2023-4-30
编辑： Sagar Karmakar 2023-5-2
Ran in:
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;

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

回答（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