|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
以下程序每次运行,a的第二个值总是结果不同,我已经再循环时清零,我不知道为什么,请高手指点,谢谢!
A=[1 0.3 0.5 1];B=[0.6 0.5 0.4 0.7];C=[];
for i=1:length(A)
C=[C,abs(B(i).^2/(1-B(i).^2))];
l=ceil(max(C));
n=length(A)*(l+1)-1;
end; C,l,n
syms b;
y2=[], j=1;
for i=1:length(A)
product=1;y1=[];
y1=[y1,A(i).*((B(i)-b).^l)];
for j=1:4
if j~=i
product=product.*((1-b.*B(j)).^(l+1));
end
end
y2=[y2,y1.*product];
end
F=[];
for m=0:n
sum=0;
for i=1:length(A),
sum=sum+diff(y2(i),b,m);
end
if (m==0),
F=sum
elseif (m>0),
F=[F,1/prod(1:m).*sum];
end
f=subs(F,b,0);
end
syms b1;
sum1=0;
for r=1:length(f)
sum1=sum1+f(r)*b1.^(r-1);
end
b1=solve(sum1);
for i=1:length(b1)
if imag(b1(i))==0
b2=b1(i);
end
end
vpa(b2,3);
syms z;
H=0;
for j=1:length(A)
H=H+A(j)*z/(z-B(j));
end
L=[];L1=[];
for k=0:l
a=0;
L=[L,sqrt(1-b2^2)*(z-b2).^k./(1-b2*z).^(k+1)*z^(-1)];
L1=[L1,sqrt(1-b2^2)*(z^(-1)-b2).^k./(1-b2*z^(-1)).^(k+1)];
f=H.*L;
r1=limit(f*(z-0.6),z,0.6);
r2=limit(f*(z-0.5),z,0.5);
r3=limit(f*(z-0.4),z,0.4);
r4=limit(f*(z-0.7),z,0.7);
r5=limit(f*z,z,0);
a=(r1+r2+r3+r4+r5);
end;a
我觉得问题应该就在最后的FOR循环里面,可是检查好像也检查不出问题所在,先谢高手指点
[ 本帖最后由 lxq 于 2007-5-30 09:33 编辑 ] |
|