请问MATLAB循环计算，结果发生飘移是为什么？
对t从9到16迭代计算，步长是1/6300，t算到9.05的时候，结果就不对了
下面的程序，中间换热计算的部分就是一些参数的定义，后面温度计算的部分算出来就不对了，感觉很简单的程序，不知道哪里有问题，求各位大佬帮帮忙
n=0;
x=0;
delta_t=1/6300;%4秒
for t=9:delta_t:16
%环境参数
C_in=900;
T_a_in=T_amb;
u_amb=1;
I=258.8*exp(-((t-4.654)/3.761)^2);
%换热计算
h_amb_g=5.7+3.8*u_amb;
q_amb_g=h_amb_g*A*(T_amb-T_g);
q_sky_g=sigma*epsilon_g*A*((T_sky+273)^4-(T_g+273)^4);
beta_f_g=1/((T_g+T_f)/2+273);
Ra_f_g=g*beta_f_g*(abs(T_g-T_f))*H^3/(nu_a*a);
if Ra_f_g<10^9
Nu_f_g=0.68+0.670*Ra_f_g^(1/4)/(1+(0.492/Pr_a)^(9/16))^(4/9);
else
Nu_f_g=(0.825+0.387*Ra_f_g^(1/4)/((1+(0.492/Pr_a)^(9/16))^(8/27)))^2;
end
h_f_g=Nu_f_g*lambda_a/H;
q_f_g=h_f_g*A*(T_f-T_g);
h_pv_g=sigma/(1/epsilon_pv+1/epsilon_g-1);
q_pv_g=h_pv_g*A*((T_pv+273)^4-(T_g+273)^4);
h_g_f=h_f_g;
q_g_f=h_g_f*A*(T_g-T_f);
beta_pv_f=1/((T_pv+T_f)/2+273);
Ra_pv_f=g*beta_pv_f*(abs(T_pv-T_f))*H^3/(nu_a*a);
if Ra_pv_f<10^9
Nu_pv_f=0.68+0.670*Ra_pv_f^(1/4)/(1+(0.492/Pr_a)^(9/16))^(4/9);
else
Nu_pv_f=(0.825+0.387*Ra_pv_f^(1/4)/((1+(0.492/Pr_a)^(9/16))^(8/27)))^2;
end
h_pv_f=Nu_pv_f*lambda_a/H;
q_pv_f=h_pv_f*A*(T_pv-T_f);
f_in=1.0;
f_out=1.5;
d=4*A_f/(2*delta_f+2*D);
Re_f=0.2*d/nu_a;
f=0.3164*Re_f^-0.25;
beta_f=1/(T_f+273);
u_f=sqrt(0.5*g*beta_f*(T_f_out-T_f_in)*H/(f_in*(A_f/A_in)^2+f_out*(A_f/A_out)^2+f*(H/d)))
q_f=rho_a*u_f*c_a*A_f*(T_f_out-T_f_in);
h_g_pv=h_pv_g;
q_g_pv=h_g_pv*A*((T_g+273)^4-(T_pv+273)^4);
h_f_pv=h_pv_f;
q_f_pv=h_f_pv*A*(T_f-T_pv);
E_pv=tau_g*eta_r*A*I*(1-0.0045*(T_pv+273.15-T_r));
beta_b_pv=1/((T_b+T_pv)/2+273);
Ra_b_pv=g*beta_b_pv*(abs(T_b-T_pv))*H^3/(nu_a*a);
if Ra_b_pv<10^9
Nu_b_pv=0.68+0.670*Ra_b_pv^(1/4)/(1+(0.492/Pr_a)^(9/16))^(4/9);
else
Nu_b_pv=(0.825+0.387*Ra_b_pv^(1/4)/((1+(0.492/Pr_a)^(9/16))^(8/27)))^2;
end
h_b_pv=Nu_b_pv*lambda_a/H;
q_b_pv=h_b_pv*A*(T_b-T_pv);
h_w_pv=sigma/(1/epsilon_w+1/epsilon_pv-1);
q_w_pv=h_w_pv*A*((T_w+273)^4-(T_pv+273)^4);
h_pv_b=h_b_pv;
q_pv_b=h_pv_b*A*(T_pv-T_b);
beta_w_b=1/((T_w+T_b)/2+273);
Ra_w_b=g*beta_w_b*(abs(T_w-T_b))*H^3/(nu_a*a);
if Ra_w_b<10^9
Nu_w_b=0.68+0.670*Ra_w_b^(1/4)/(1+(0.492/Pr_a)^(9/16))^(4/9);
else
Nu_w_b=(0.825+0.387*Ra_w_b^(1/4)/((1+(0.492/Pr_a)^(9/16))^(8/27)))^2;
end
h_w_b=Nu_w_b*lambda_a/H;
q_w_b=h_w_b*A*(T_w-T_b);
%自然对流
f_in=1.0;
f_out=1.5;
d=4*A_b/(2*delta_b+2*D);
Re_b=0.2*d/nu_a;%按u=0.2来算
f=0.3164*Re_b^-0.25;
beta_b=1/(T_b+273);
u_b=(1/100)*sqrt(0.5*g*beta_b*(T_b_out-T_b_in)*H/((f_in*(A_b/A_in)^2+f_out*(A_b/A_out)^2+f*(H/d))/10000))
q_b=rho_a*u_b*c_a*A_b*(T_b_out-T_b_in);
h_pv_w=h_w_pv;
q_pv_w=h_pv_w*A*((T_pv+273)^4-(T_w+273)^4);
h_b_w=h_w_b;
q_b_w=h_b_w*A*(T_b-T_w);
h_w_amb=29;%
R=delta_w/lambda_w+1/h_w_amb;
q_amb_w=(T_amb-T_w)*A/R;
%温度计算
T_1=(q_amb_g+q_sky_g+q_f_g+q_pv_g+alpha_g*A*I)*(delta_t*3600)/(m_g*c_g);
T_2=(q_g_f+q_pv_f-q_f)*(delta_t*3600)/(m_f*c_a);
T_3=(q_g_pv+q_f_pv+tau_g*alpha_pv*A*I-E_pv+q_b_pv+q_w_pv)*(delta_t*3600)/(m_pv*c_pv);
T_4=(q_pv_b+q_w_b-q_b)*(delta_t*3600)/(m_b*c_a);
T_5=(q_pv_w+q_b_w+q_amb_w)*(delta_t*3600)/(m_w*c_w);
T_g=T_g+T_1;
T_g=vpa(T_g,7)
T_f=T_f+T_2;
T_f=vpa(T_f,7)
T_pv=T_pv+T_3;
T_pv=vpa(T_pv,7)
T_b=T_b+T_4;
T_b=vpa(T_b,7)
T_w=T_w+T_5;
T_w=vpa(T_w,7)
T_f_out=2*T_f-T_f_in;
T_b_out=2*T_b-T_b_in;
if u_f==0||u_b==0
T_out=T_b_out;
else
T_out=(u_f*A_f*T_f_out+u_b*A_b*T_b_out)/(u_f*A_f+u_b*A_b);
end
plot(t,T_f,'*r');
hold on;
plot(t,T_b,'*y');
hold on;
plot(t,T_out,'*g');
hold on;
end