matlab saturate,利用MATLAB求解一微偏微分方程

本帖最后由 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 上传