24小时热门版块排行榜

返回列表

chen17307

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 56.2
散金: 22
帖子: 12
在线: 9.8小时
虫号: 2932786
注册: 2014-01-15
性别: GG
专业: 计算数学与科学工程计算

[求助] 求助各路大神一个偏微分方程的数值计算问题（FPK方程）

求助各位虫友大神，用MATLAB写一个求解偏微分方程的程序，方程和离散化的推导过程如图1和图2，程序贴在下面了，能跑出结果但是结果不对，与精确解的结果不相符合啊。不知道哪里错了。跪求各路大神帮忙看看，小弟不胜感激啊。。。

%隐式差分格式求解FPK方程
%二阶差分格式
clear;clc;
%% 参数设置
Bata=0.1;D=0.01;delta=4;omega0=1;%参数设置

min_x=-1.5;
max_x=1.5;
min_y=-1.5;
max_y=1.5; %坐标轴范围

Nx=300;
Ny=300; %迭代步数

sx=(max_x-min_x)/Nx;
sy=(max_y-min_y)/Ny; %迭代步长

x=min_x:sx:max_x;
y=min_y:sy:max_y; %网格划分

for i=1:Nx+1
for j=1:Nx+1
   h(i,j)=Bata*y(j)-omega0*x(i)+delta*(x(i))^3; %化简式中的h(x,y)
end
end
%% 求精确解
%其精确解为p(x,y)=c*exp(-Bata/D*(y^2/2-omega0^2/2*x^2+delta/4*x^4))

%---精确解---%
for k=1:length(x)
for l=1:length(y)
      p(k,l)=exp(-Bata/D*(y(l)^2/2-omega0/2*x(k)^2+delta/4*x(k)^4));
end
end
%精确解的图像
figure(1)
mesh(x,y,p)
title('精确解')
xlabel('x','FontSize',12)
ylabel('y','FontSize',12)
zlabel('p(x,y)','FontSize',12)
axis([-1.5,1.5,-1.5,1.5,0,2])
%% 求数值解

%---赋初值和边值---%
%---初值和边界值来自求得的精确解的结果
p1=zeros(Nx+1);
p1(:,1)=p(:,1);%左边值
p1(:,end)=p(:,end);%右边值
p1(1,

=p(1,

;%初值

for ii=1:Nx
%---线性方程组AX=B的系数矩阵A---%
%A为三对角阵，主对角：Diag；上对角：Updiag
%下对角：Lowdiag
Diag=zeros(1,Nx-1);
%Up_diag=D.*ones(1,Nx-2);
%Low_diag=D.*ones(1,Nx-2);
Up_diag=D/sy^2.*ones(1,Nx-2);
Low_diag=D/sy^2.*ones(1,Nx-2);
for jj=1:Ny-1
      Diag(jj)=-2*D/sy^2-y(jj+1)/sx;
end

%---B---%
B=zeros(Nx-1,1);
for kk=1:Ny-1
   % B(kk,1)=(y(kk+1)*sx-h(ii,kk+1)*sy+Bata*sy^2)*p1(ii,kk+1)+h(ii,kk+1)*sy*p1(ii,kk+2);
   B(kk,1)=-((y(kk+1)/sx-h(ii,kk+1)/sy+Bata)*p1(ii,kk+1)+h(ii,kk+1)/sy*p1(ii,kk+2));
end
%B(1,1)= B(1,1)-D*p1(ii+1,1);
% B(end,1)= B(end,1)-D*p1(ii+1,end);
B(1,1)= B(1,1)-D/sy^2*p1(ii+1,1);
B(end,1)= B(end,1)-D/sy^2*p1(ii+1,end);
p1(ii+1,2:Nx)=zhuiganfa(Low_diag,Diag,Up_diag,B); %追赶法解线性方程组
% p1(ii+1,

=( p1(ii+1,

-min(p1(ii+1,

))./(max(p1(ii+1,

)-min(p1(ii+1,

));
%p1(p1<0)=0;
end
figure(2)
mesh(x,y,p1)
xlabel('x','FontSize',12)
ylabel('y','FontSize',12)
zlabel('p(x,y)','FontSize',12)