How to solve a Fokker Planck equation in both space and time using finite difference method
9 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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
Antworten (0)
Siehe auch
Kategorien
Mehr zu Eigenvalue Problems finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!

