%**********************************一单元测试**********************************
syms c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c4_ c5_ c6_ real;
syms alpha1 alpha2 alpha3 alpha4 alpha5 alpha6 real;
syms L T0 rh r0 E S I positive;
syms sga yta lda lda0
syms Q1 R1 K1 Q2 R2 K2 tol
lda=0.05011;
%L=10;S=100;I=1;E=20;
sga=100;yta=0.1;r0=0;
%**********************************代数运算************************************
alpha1=yta^2+lda^2;alpha2=alpha1*sga^2;alpha3=yta^2*sga^2;alpha4=lda^2;alpha5=0;alpha6=1/(2*lda*yta);
c1=alpha1-(alpha3*(r0+0.5)+alpha5);c2=alpha3*r0;c4=alpha3*r0;c3=0.5*alpha3;c5=alpha3;
c6=-alpha1*(alpha3*(r0+0.5)+alpha5)-alpha2;c7=alpha1*alpha3*r0;c9=alpha1*alpha3*r0;
c8=0.5*alpha1*alpha3;c10=alpha1*alpha3;c11=-(alpha1*alpha2-4*alpha3*alpha4);
c4_=3*c4;c5_=3*c5;
c6_=3*alpha3+c6;
a=sym(zeros(250,1));
b=sym(zeros(250,1));
delta=sym(zeros(6,1));
F=sym(zeros(6,1));
for i=1:250 %符号向量
a(i)=sym (['a',num2str(i)]);
end
for i=1:250 %符号向量
b(i)=sym (['b',num2str(i)]);
end
% a=sym('a%d',[200,1]);
% b=sym('b%d',[200,1]);
for ii=6:249
a(ii+1)=-(a(ii-1)*(ii-2)*(ii-3)*(ii-4)*(ii-5)*c1+a(ii-2)*(ii-3)*(ii-4)*(ii-5)*(c2*(ii-6)+c4)+a(ii-3)*(ii-4)*(ii-5)*(c3*(ii-6)*(ii-7)+c5*(ii-6)+c6)+a(ii-4)*(ii-5)*(c7*(ii-6)+c9)+a(ii-5)*(c8*(ii-6)*(ii-7)+c10*(ii-6)+c11))/(ii*(ii-1)*(ii-2)*(ii-3)*(ii-4)*(ii-5));
end