请大家帮我看看错在哪里

markzhou

请大家帮我看看错在哪里。（我自己水平太臭，都搞了3天了，还是找不出错误！）
先谢谢了！
程序如下，出现的错误提示如下：
function F = myfun(x)%一个函数文件
g=9.81;
Ta=298.15;
HT=1000;
Ub2=6.5;Vw=5;Cps=3000;Cpf=1007;Hc=200;Ta=298.15;
Tg=Ta;
a=0.798;%a=(τα)2
o=0.0000000567;%σ
R=122;Acoll=pi*R^2;
Rc=10.16;Ac=pi*Rc^2;
h1=9.857956412;
h2=24.13074832;
Ub2=6.430030004;
Vs=0.1;%可变
Hs=0.1;%可变
ms=2*pi*R*Hs*Vs;
Ts0=273.15+80;
ep=0.94;%ep=εp
ec=0.9;%ec=εp
%Tf=273.15+25;%事先给定的
r=0.009753;
ra=0.0065;%ra=r∞
ff=0.76-0.118*Vw+0.0066*Vw^2;
hw=2.8+3*Vw;

y6=(0.2106+0.0026*(1.1614-0.00353*(2*x(4)-Ta-300))*x(11)*Ac/((1.1614-0.00353*(x(4)-300))*pi*R*4)*((1.1614-0.00353*(x(4)-300))*(x(3)+x(4))/2/((1.846+0.00472*(x(4)-300))*0.00001*g*(x(3)-x(4))))^(1/3))/((1.846+0.00472*(x(4)-300))*0.00001*x(4)/(g*(x(3)-x(4))*Cpf*(0.0263+0.000074*(x(4)-300))^2*(1.1614-0.00353*(x(4)-300))^2))^(1/3);%y6=h3;p=ρ,u=μ
y7=(0.2106+0.0026*Vw*((1.1614-0.00353*(x(4)-300))*(x(4)+x(5))/2/((1.846+0.00472*(x(4)-300))*0.00001*g*(x(4)-x(5))))^(1/3))/((1.846+0.00472*(x(4)-300))*0.00001*x(4)/(g*(x(4)-x(5))*Cpf*(0.0263+0.000074*(x(4)-300))^2*(1.1614-0.00353*(x(4)-300))^2))^(1/3);%y7=h4;p=ρ,u=μ
y8=o*(x(3)^2+x(1)^2)*(x(1)+x(1))/(1/ep+1/ec-1);%hrp1p2=y8
y9=o*(x(3)^2+x(5)^2)*(x(3)+x(5))/(1/ec+1/ec-1);%hrp2c=y9
y10=1/(1/1.2529*((x(3)-Ta)/(1+ff))^(-0.25)+1/hw)+o*(x(3)^2+Ta^2)*(x(3)+Ta)/(1/ep+(1+ff)/ep-1);%Ut2=y10,c=σ,ep=εp
y11=(2*(0.00353*g*((2*x(4)-Ta-Ta)*Hc-(r-ra)*Hc^2/2))/(3*(1.1614-0.00353*(2*x(4)-Ta-300))))^(1/2);%Vout=y11
y12=(1.1614-0.00353*(2*x(4)-Ta-300))*Ac*x(11);%mf=y12,pout=ρcoll,out
y1=(a*HT+h1*x(2)+x(8)*x(3)+Ub2*Tg)/(x(8)+Ub2+h1);%y1=Tp1,a=(τα)2
y2=(2*ms*Cps*Ts0/Acoll+h1*x(1)+h2*x(3))/(2*ms*Cps/Acoll+h1+h2);%y2=Ts
y3=(h2*x(2)+x(8)*x(1)+x(6)*x(4)+x(9)*x(5))/(h2+x(8)+x(6));%y3=Tp2
y4=(x(6)*x(3)+x(7)*x(5)+2*x(12)*Cpf*Ta/Acoll)/(2*x(12)*Cpf/Acoll+x(6)+x(7));%y4=Tf
y5=(x(7)*x(4)+x(9)*x(3)+x(10)*Ta)/(x(7)+x(9)+x(10));%y5=Tc

>> [x,y,h]=fsolve(@myfun,[330 340 330 315 305 10 10 6 6 6 6 10000])
??? Output argument "F" (and maybe others) not assigned during call to "E:\20061205\myfun.m (myfun)".

Error in ==> myfun at 2
g=9.81;

Error in ==> fsolve at 180
        fuser = feval(funfcn{3},x,varargin{:});


==========================eight=================================================
帮你修改了一下标题，类似“（我自己水平太臭，都搞了3天了，还是找不出错误！）”这样的字眼还是不要写在标题上为好
===============================================================================

[ 本帖最后由 eight 于 2007-1-12 19:19 编辑 ]

happy

myfun最后加上

F=[y1;y2;y3;y4;y5;y6;y7;y8;y9;y10;y11;y12];

coldspring

[quote]原帖由 happy 于 2006-12-6 11:13 发表
myfun最后加上

F=[y1;y2;y3;y4;y5;y6;y7;y8;y9;y10;y11;y12];

wf1254

好像确实有点臭啊，呵呵，这样的错误怎会看不出？

markzhou

谢谢！
我把它改了，因为它还要和复杂一些的程序结合，不过对改后的程序进行计算，算出来的结果既然都不是实数，这不知道为什么，用excel迭代可以算出这个方程组的答案阿!请大虾帮小弟想想办法。
Ta=298.15;Tg=298.15;
Vw=5;
cps=3000;
cpf=1007;
Hc=200;
Hm=4;%collector 中间距离地面高度
R=122;Acoll=pi*R^2;
Rc=10.16;Ac=pi*Rc^2;
Ub2=0.117;  %Ub=6.5; //1/(1/Ub+1/hg1) hg1 9.857956412/地面一层/注意:不能混淆
o=0.0000000567;%σ,Stefan-Boltzmann常数
g=9.81;
h1=24.13074832;
Vs=0.8;%可变
Hs=0.5;%可变
Ts0=353.15;%给定的,80
ep=0.94;%ep=εp
ec=0.9;%ec=εp
r=0.009753;
ra=0.0065;%ra=r∞  
ms=2*pi*R*Hs*Vs;  
ff=0.76-0.118*Vw+0.0066*power(Vw,2);
hw=2.8+3*Vw;

Tf=305.12;Ts=323;Tp=316.6;Tc=299.05;Vout=3.803;%事先给定的,25
for i=1:100
p=(1.1614-0.00353*(Tf-300));%p=ρ
u=(1.846+0.00472*(Tf-300))*0.00001;%u=μ
k=(0.0263+0.000074*(Tf-300));
Tout=2*Tf-Ta;
pout=(1.1614-0.00353*(Tout-300));%p=ρ
mf=pout*Ac*Vout;
Vm=mf/(pi*R*Hm*p);%事先给定:以R/2的速度来表示

h2=(0.2106+0.0026*Vm*power(p*(Tp+Tf)/2/(u*g*(Tp-Tf)),1/3))/power(u*(Tp+Tf)/2/(g*(Tp-Tf)*cpf*k*k*p*p),1/3);
h3=(0.2106+0.0026*Vm*power(p*(Tc+Tf)/2/(u*g*(Tf-Tc)),1/3))/power(u*(Tc+Tf)/2/(g*(Tf-Tc)*cpf*k*k*p*p),1/3);
h31=(0.2106+0.0026*Vm*power(p*(Tc+Tf)/2/(u*g*(Tf-Tc)),1/3));
h32=power(u*(Tc+Tf)/2/(g*(Tf-Tc)*cpf*k*k*p*p),1/3);
h321=u*(Tc+Tf)/2/(g*(Tf-Tc)*cpf*k*k*p*p);
hrpc=o*(power(Tp,2)+power(Tc,2))*(Tp+Tc)/(1/ep+1/ec-1);
Ut=1/(1/1.2529*power((Tc-Ta)/(1+ff),-0.25)+1/hw)+o*(power(Tc,2)+power(Ta,2))*(Tc+Ta)/(1/ec+(1+ff)/ec-1);
Tc2=((hrpc*(h2+h3+2*mf*cpf/Acoll)+h2*h3)*Tf+(h2*Ut-hrpc*2*mf*cpf/Acoll)*Ta)/((hrpc+h3+Ut)*h2+hrpc*h3);
Ts2=(h1*h2*Tf/(h1+h2+hrpc)+h1+hrpc*Tc/(h1+h2+hrpc)+0.117*Tg+2*ms*cps*Ts0/Acoll)/(h1+0.117+2*ms*cps/Acoll-h1*h1/(h1+h2+hrpc));
Tp2=(h1*Ts+h2*Tf+hrpc*Tc)/(h1+h2+hrpc);
Tf2=(h2*Tp+h3*Tc+2*mf*cpf/Acoll*Ta)/(h2+h3+2*mf*cpf/Acoll);
Vout2=power(2*(0.00353*g*((2*Tf2-2*Ta)*Hc-(r-ra)*Hc*Hc/2))/(3*p),0.5);

        if Tf2<=Tf+0.01 && Tf2>=Tf-0.01 && Ts2<=Ts+0.01 && Ts2>=Ts-0.01  && Tp2<=Tp+0.01 && Tp2>=Tp-0.01  && Tc2<=Tc+0.01 && Tc2>=Tc-0.01  && Vout2<=Vout+0.01 && Vout2>=Vout-0.01
                disp('yes')
            break;

        else
disp('no')
Tf=Tf2;Ts=Ts2;Tp=Tp2;Tc=Tc2;Vout=Vout2;
                continue;
end
end