Labfans是一个针对大学生、工程师和科研工作者的技术社区。 论坛首页 | 联系我们(Contact Us)
MATLAB爱好者论坛-LabFans.com
返回   MATLAB爱好者论坛-LabFans.com > 工程数学软件 > MATLAB论坛
MATLAB论坛 一切MATLAB相关问题在此讨论。
 
 
主题工具 显示模式
旧 2007-08-16, 09:51   #1
水电校小张
初级会员
 
注册日期: 2007-08-16
帖子: 3
声望力: 0
水电校小张 正向着好的方向发展
默认 【求助】Matlab

求列位大侠多帮忙看看,错误很多,但公式正确。
我刚学MATLAB,定义函数,变量,数组和读入读出数据很迷惑。教材和书上的好像我都试过了,不好用。如果能指点一下基本技巧就更好了。
%smx.m
function f=smx(z,Q,s,z0,c1,c2)
syms A v P R Bs b Q Q0 Qp z z0 zp g ic n ns h t time s m
h=sym('h');z=sym('z');z0=sym('z0');s=sym('s');h=sym('h');
time=0;t=60;rr=1;n=11;g=9.8;Q=6*ones(1,n);ns=0.017;m=3;ic=0.0002;
zp=[];Qp=[];
while (rr>1.0e-3)
zp(1)=15.518;
time=time+t;
Qp(n)=30+0.1*time;
if Qp(n)>150.0
Qp(n)=150.0;
end
for i=1:n
h=z(i)-z0(i);
Bs=b+2*m*h;
A=(b+m*h)*h;
v=Q(i)/A;
P=b+(1+m^2)^0.5*h/2;
R=A/P;
c1=t-s(i)/(abs(v+g*h))^0.5;
c2=t-s(i)/(abs(v-g*h))^0.5;
if (c1<0|c2<0)
t=3/4*t;'warnning!'
end
N1=g*A-Bs*v^2;
N2=Bs*ic*v^2-g*(v*ns)^2*P^(4/3)/A^(1/3);
switch(i)
case 1
Qp(i)=Q(1)-2*v*t*(Q(i+1)-Q(i))/s(i+1)...
-N1*t*(z(i+1)-z(i))/s(i+1)+N2*t...
-t*(g/A/Bs)^0.5*((zp(i)-z(i))/t...
+(Q(i+1)-Q(i))/s(i+1));
case n
zp(i)=z(i)+t/(g/A*Bs)^0.5*((Qp(i)-Q(i))/t...
+2*v*(Q(i)-Q(i-1))/s(i)+...
N1*(z(i)-z(i-1))/s(i)-N2)-t*(Q(i)-Q(i-1))/s(i);
otherwise
z0=ar/2*(z(i+1)-z(i-1))+(1-ar)*z(i);
Q0=ar/2*(Q(i+1)-Q(i-1))+(1-ar)*Q(i);
zp(i)=z0-t/Bs*(Q(i+1)-Q(i-1))/(s(i)+s(i+1));
Qp(i)=Q0-2*v*t*(Q(i+1)-Q(i-1))/(s(i)+s(i+1))...
-N1*t*(z(i+1)-z(i-1))/(s(i)+s(i+1))+N2*t;
end
save smx zp(n),Qp(1)
for i=1:n
z(i)=zp(i);Q(i)=Qp(i);
end
end
rr=abs(Qp(1)-Q(1))
end
plot(zp,time,'-r')
hold on
plot(Qp,time,'。')



比如举一个比较全例子,用%加以说明各项是干什么的,就不胜感激了。
水电校小张 当前离线   回复时引用此帖
 


发帖规则
不可以发表新主题
不可以发表回复
不可以上传附件
不可以编辑自己的帖子

启用 BB 代码
论坛启用 表情符号
论坛启用 [IMG] 代码
论坛禁用 HTML 代码



所有时间均为北京时间。现在的时间是 11:07


Powered by vBulletin
版权所有 ©2000 - 2025,Jelsoft Enterprises Ltd.