|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
function ex
%--------------------------------------------------------------------------
L = 50;
s1 = 0.5;
s2 = 0;
T = 5;
qr = 0.1699;
f = 0.4432;
a = 0.0526;
n = 1.3924;
ks = 31.6;
aa=0.0037;
bb=14.029;
x = linspace(0,L,50);
t = linspace(0,T,25);
options=odeset('RelTol',1e-4,'AbsTol',1e-4,'NormControl','off','InitialStep',1e-7)
sol = pdepe(0,@unsatpde,@unsatic,@unsatbc,x,t,options,s1,s2,qr,f,a,n,ks);
u1=real(sol(:,:,1))
figure
surf(x,t,u1)
% -------------------------------------------------------------------------
function [c,f,s] = unsatpde(x,t,u,DuDx,s1,s2,qr,f,a,n,ks)
[u,k,Dth] = sedprop(u,qr,f,a,n,ks);
c=[1];
f=[Dth.*DuDx-k];
s = [0];
% -------------------------------------------------------------------------
function u0 = unsatic(x,s1,s2,qr,f,a,n,ks)
u0 = [0.4432];
% -------------------------------------------------------------------------
function [pl,ql,pr,qr] = unsatbc(xl,ul,xr,ur,t,s1,s2,qr,f,a,n,ks)
pl = [-0.5];
ql = [-1];
pr = [ur(1)-s2];
qr = [0];
%-------------------------------------------------------------------------
function [u,k,Dth] = sedprop(u,qr,f,a,n,ks)
m = 1-1/n;
if u<qr
c=1e-20;
k=0;
u=qr;
Dth=aa*exp(bb*qr);
elseif u>f
c=1e-20;
k=ks;
u=f;
Dth=aa*exp(bb*f);
else
k=ks*((u-qr)/(f-qr))^0.5*(1-(1-((u-qr)/(f-qr))^(1/m))^m)^2;
Dth=aa*exp(bb*u);
end
我想通过子程序控制u的范围在[qr,f]之间,即[0.1699,0.4432],但是u的输出结果很多不在这个范围之间,是怎么回事?
如:
输出结果:
Columns 49 through 50
0.4432 0
0.8439 0
0.8176 0
0.7902 0
0.7748 0
0.7645 0
0.7491 0
0.7385 0
0.7302 0
0.7236 0
0.7168 0
0.7078 0
0.7037 0
0.7016 0
0.7017 0
0.7053 0
0.7054 0
0.7036 0
0.7064 0
0.7033 0
0.7060 0
0.7040 0
0.7058 0
0.7036 0
0.7058 0 |
|