本帖最后由 spritecoca 于 2010-8-24 21:09 编辑
利用MATLAB 函数pdepe求解一维偏微分方程中遇到问题:弄了好几天,还是没结果。求高手指出程序中错误或者给出该方程的差分解法,本人不胜感激。
方程形式如下:
Du/dt=-dQ/dx+s (x,t) x=[0,1500] t=[0,30] ;
其中 : Q=au^m, 当t
当t>tp时,s(x,t)=p-St-0.5-A ;
式中: a,m tp,S,A均为常数;
初始条件: u(x,0)= 0;
边界条件 : 左边界: u(0,t)=0 ; 右边界: u(xr,t)=[xr×( p-St-0.5-A)/a]m,m为常数;
本人编写程序如下,但程序有错误,在右边界点方程解出现波动,调整x在右边界附近的值(加密)任没有作用,还有就是系数函数m文件中s项和边界函数m文件中的qr项似乎都有问题,改变那些项目的值(如把qr=ur和pr=ur-(1500*(0.1*cos(15*pi/180)-S*t(i).^(-0.5)-A)/alfa).^0.6)两者对结果没有任何影响,不知道是哪里出了错,盼望高手解答,不人不胜感激!
clc
clear
%——---------------------------常数部分
M=0;
p=0.1; % rainfall tensity cm/min
Ks=0.007;% saturate hydraulic conductivity cm/min;
Sf=60;% wet front matric
ts=0.36;t0=0.06; % saturate moisture content and restuit
ti=0.10; %current water content
Fp=Sf*(ts-ti)./(p/Ks-1);
tp=Fp/p;
S=sqrt(2*Ks*(ts-ti)); % philip formula parameter soropity
A=1/2*Ks; % Philip formula parameter
J=sin(15*pi/180)^1/2;
n=0.025;
alfa=(1/n)*J;
m=5/3;
%---------------------------------
x=linspace(0,1500,100);
t=linspace(0.1,30.1,30); % rainfal time 30min
save timestep tp t S A alfa
sol=pdepe(M,@KenePara,@KeneInt,@KentBc,x,t);
u=sol(:,:,1);
surf(x,t,u);
xlabel('Ƴ¤/m');
ylabel('ʱ¼ä/min');
zlabel('ÆÂÃæË®Éî/h');
%--------------------------------------------------- ---------------系数c,f, s
function [c,f,s]=KenePara(x,t,u,dudx);
load timestep % load tp alfa S A
if nargin==2;
u=0;
dudx=0;
end
c=1; % ----------------------------paremeter c
f=-alfa*u.^(5/3); % --------------------parameter f
N=size(t,2);
for i=1:N;
if t(i)
s(i,:)=0;
else
s(i,:)=0.1*cos(15*pi/180)-S*t(i).^(-0.5)-A;
end
end
%-------------------------------------------边界条件
function [pl,ql,pr,qr]=KentBc(xl,ul,xr,ur,t)
load timestep % load tp t S A alfa
pl=ul;ql=0;
N=size(t,2);
for i=1:N
if t(i)
pr=ur;
qr=0;
else
pr(i,:)=ur-(1500*(0.1*cos(15*pi/180)-S*t(i).^(-0.5)-A)/alfa).^0.6;
qr=0;
end
end
%_-----------------------------------------------初始条件
Function h=KeneInt(x)
h=0;
问题已解决..................,利用差分方程求解。
untitled.jpg
(45.46 KB, 下载次数: 10)
2010-8-20 10:52 上传