有会matlab迭代运算的吗

[dlong411]

有MATLAB高手吗，下面是我的一程序，不过那个迭代循环有问题，请高手们指点下啊，
程序片段：
Fi=[x1-F1*x1*(x1^2+x2^2)^((q-1)/2)+F2*x2*(x1^2+x2^2)^((p-1)/2);-x2 +F1*x2*(x1^2+x2^2)^((q-1)/2)+F2*x1*(x1^2+x2^2)^((p-1)/2)]-C;
Fxi=[diff(Fi,x1,1) diff(Fi,x2,1)];
x1=0;
x2=0;
% 迭代公式：
for k=1:4
Fi
F=subs(Fi);
Fx=subs(Fxi);
Fxx=inv(Fx);
X=subs(X)-Fxx*F
x1=X(1);
x2=X(2);
end

请会的高手赐教，谢了！




程序全部如下：
clear all
% 基本参数
miu=0.4;
m=1500.0;
K=2.79e5;
C=4.10e3;
p=1;
q=1;
roub=7800;
wide=0.2;
high=0.01;
L=3;
Eb=2.1e11;
Ab=wide*high;
mb=roub*Ab*L;
Ib=wide*high^3/12;
xc=miu*L;
x0=(1/K)^(1/q);
omiga0=sqrt(K*x0^(q-1)/m); % ω0
ksi=C*x0^(p-1)*omiga0^(p-2)/(2*m); % ξ
% xcx=xc/x0;
% roubx=x0*roub/m;
i=0;
for f=10;
i=i+1;
omiga=2*pi*f; % kf=sqrt(omiga*sqrt(roum*Am/(Emx*Im)));
syms x
n=1:1:2;
lambdan=pi*(n+0.5)/L;
betan=(cos(lambdan*L)-cosh(lambdan*L))./(sin(lambdan*L)-sinh(lambdan*L));
Wn=cosh(lambdan*x)-cos(lambdan*x)+betan.*(sin(lambdan*x)-sinh(lambdan*x)); %获得一个数组最大值取为8
Wn2=diff(Wn,x,2);
Kn=int(Eb*Ib*Wn2.^2,0,L);
Mn=int(roub*Wn.^2,0,L);
Kn=subs(Kn,xc,'x');
Mn=subs(Mn,xc,'x');
Wn=subs(Wn,xc,'x');
omigan=sqrt(Kn./Mn);
Omiga=omiga/omiga0; % Ω Omiga
Omigan=omigan./omiga0;
gamman=m./Mn; % γ gamma
% 求迭代公式：
syms s x1 x2
% int(exp(-s)*s^(q/2),0,inf);
lambdaq=double(2*int(exp(-s)*s^(q/2),0,inf)/(sqrt(pi)*int(exp(-s)*s^(q/2+0.5),0,inf))); % lambda lambda
lambdap=double(2*int(exp(-s)*s^(p/2),0,inf)/(sqrt(pi)*int(exp(-s)*s^(p/2+0.5),0,inf)));
delta1=sqrt(x1^2+x2^2); % δ delta
an=Omigan.^2-Omiga.^2;
an1=gamman.*Wn.*lambdaq;
an2=2*ksi*Omiga^p*gamman.*Wn.*lambdap;
An11=-an1*delta1^(q-1);
An12=an2*delta1^(p-1);
An21=an2*delta1^(p-1);
An22=an1*delta1^(q-1);
An=[An11 An12;An21 An22];
% AAn11=An11.*Wn./an;
% AAn12=An12.*Wn./an;
% AAn21=An21.*Wn./an;
% AAn22=An22.*Wn./an;
f1=sum(an1.*Wn./an);
f2=sum(an2.*Wn./an);
% AAn=[sum(An11) sum(An12);sum(An21) sum(An22)];
bn1=lambdaq*Omiga.^(-2);
bn2=2*ksi*Omiga.^(p-2)*lambdap;
% Bn=[1-bn1*delta1^(q-1) bn2*delta1^(p-1);bn2*delta1^(p-1) -1+bn1*delta1^(q-1)];
F1=bn1+f1;
F2=bn2+f2;
C=[-Omiga.^(-2);0];
X=[x1;x2];
Fi=[x1-F1*x1*(x1^2+x2^2)^((q-1)/2)+F2*x2*(x1^2+x2^2)^((p-1)/2);-x2+F1*x2*(x1^2+x2^2)^((q-1)/2)+F2*x1*(x1^2+x2^2)^((p-1)/2)]-C; %Fi=(Bn-AAn)*X-C
Fxi=[diff(Fi,x1,1) diff(Fi,x2,1)];
x1=0;
x2=0;
% 迭代公式：
for k=1:4
Fi
F=subs(Fi);
Fx=subs(Fxi);
% F=double(F);
% Fx=double(Fx);
Fxx=inv(Fx);
X=subs(X)-Fxx*F
x1=X(1);
x2=X(2); % F=subs(F,[1,3],'[x1,x2]')
end
An11=subs(An11);
An12=subs(An12);
An21=subs(An21);
An22=subs(An22);
qn1=-An11.*x1/an-An12.*x2/an;
qn2=-An21.*x1/an-An22.*x2/an;
% Qn=-[An11.*x1/an+An12.*x2/an;An21.*x1/an+An22.*x2/an];
% 功率流公式：
delta1=sqrt(X(1)^2+X(2)^2);
R0=0;
R1=lambdaq*delta1^q;
R2=-2*ksi*Omiga^p*lambdap*delta1^p;
R=sqrt(R1^2+R2^2);
thetafc=atan(R2/R1);
thetaw=atan(sum(Wn*qn2)./sum(Wn*qn1)); % θ theta
fai=atan(X(2)/X(1)); % φ fai
thetaz=fai-thetafc;
% deltaz=R/Omiga^2;
deltaw=sqrt((sum(Wn*qn2)).^2+(sum(Wn*qn1)).^2);
% XX(i)=delta1
Ps(i)=abs(-0.5*Omiga*(delta1*sin(fai)-deltaw*sin(thetaw)));
Pb(i)=abs(-0.5*Omiga*R*deltaw*sin(thetaz+thetaw));
end
f=logspace(0,3,30);
PPs=10*log(Ps/1e-12);
PPb=10*log(Pb/1e-12);
semilogx(f,PPs,'k-',f,PPb,'b-');
% legend('输入机器功率','输入基础梁的功率','f=[fy fz Mx]=fz*[1;1;1]');
title('Power Flow Through Inclined Isolated System');
xlabel('FREQUENCE(rad/s)');ylabel('Q(dB)');
hold on;
% grid
end