cxwxy1
2008-04-19, 11:05
我是新学MATLAB,我想编一个用牛顿法解一个非线性方程组,14个迭代方程、14个未知数。下面是我写的程序,但无法运行,始终不明白为什么。求助各位高手,帮忙指点一下。见笑了,谢谢!
function y=G(maxit,q0)
f=0.505; m=1.75; b=4.75;
a=-0.0891; c=-0.0775; x=0.5575; k=0.3152;
L=0.8; D=14.87; r=1.35; ct=14.20; dz=0.008;
Q0=934.20;
h(15)=27.20;
maxit=10;
q=0
Q=Q0+q;
h=((q/(k*r^a*ct^c)))^(1/x);
Kp=(D^b/(f*L))*Q^(1-m);
Ke=(k*r^a*ct^c)^(1/x)*q^(1-1/x);
for K=1:maxit
y(1)=(Kp(1)+Ke(1))*h(1)-Kp(1)*h(2)-Kp(1)*dz;
y(2)=-Kp(1)*h(1)+(Kp(1)+Kp(2)+Ke(2))*h(2)-Kp(2)*h(3)-(Kp(2)*dz-Kp(1)*dz);
y(3)=-Kp(2)*h(2)+(Kp(2)+Kp(3)+Ke(3))*h(3)-Kp(3)*h(4)-(Kp(3)*dz-Kp(2)*dz);
y(4)=-Kp(3)*h(3)+(Kp(3)+Kp(4)+Ke(4))*h(4)-Kp(4)*h(5)-(Kp(4)*dz-Kp(3)*dz);
y(5)=-Kp(4)*h(4)+(Kp(4)+Kp(5)+Ke(5))*h(5)-Kp(5)*h(6)-(Kp(5)*dz-Kp(4)*dz);
y(6)=-Kp(5)*h(5)+(Kp(5)+Kp(6)+Ke(6))*h(6)-Kp(6)*h(7)-(Kp(6)*dz-Kp(5)*dz);
y(7)=-Kp(6)*h(6)+(Kp(6)+Kp(7)+Ke(7))*h(7)-Kp(7)*h(8)-(Kp(7)*dz-Kp(6)*dz);
y(8)=-Kp(7)*h(7)+(Kp(7)+Kp(8)+Ke(8))*h(8)-Kp(8)*h(9)-(Kp(8)*dz-Kp(7)*dz);
y(9)=-Kp(8)*h(8)+(Kp(8)+Kp(9)+Ke(9))*h(9)-Kp(9)*h(10)-(Kp(9)*dz-Kp(8)*dz);
y(10)=-Kp(9)*h(9)+(Kp(9)+Kp(10)+Ke(10))*h(10)-Kp(10)*h(11)-(Kp(10)*dz-Kp(9)*dz);
y(11)=-Kp(10)*h(10)+(Kp(10)+Kp(11)+Ke(11))*h(11)-Kp(11)*h(12)-(Kp(11)*dz-Kp(10)*dz);
y(12)=-Kp(11)*h(11)+(Kp(11)+Kp(12)+Ke(12))*h(12)-Kp(12)*h(13)-(Kp(12)*dz-Kp(11)*dz);
y(13)=-Kp(12)*h(12)+(Kp(12)+Kp(13)+Ke(13))*h(13)-Kp(13)*h(14)-(Kp(13)*dz-Kp(12)*dz);
y(14)=-Kp(13)*h(13)+(Kp(13)+Kp(14)+Ke(14))*h(14)-Kp(14)*h(15)-(Kp(14)*dz-Kp(13)*dz);
J=jacobian(y,q);
dq=-J\f;
q=q+dq;
end
function y=G(maxit,q0)
f=0.505; m=1.75; b=4.75;
a=-0.0891; c=-0.0775; x=0.5575; k=0.3152;
L=0.8; D=14.87; r=1.35; ct=14.20; dz=0.008;
Q0=934.20;
h(15)=27.20;
maxit=10;
q=0
Q=Q0+q;
h=((q/(k*r^a*ct^c)))^(1/x);
Kp=(D^b/(f*L))*Q^(1-m);
Ke=(k*r^a*ct^c)^(1/x)*q^(1-1/x);
for K=1:maxit
y(1)=(Kp(1)+Ke(1))*h(1)-Kp(1)*h(2)-Kp(1)*dz;
y(2)=-Kp(1)*h(1)+(Kp(1)+Kp(2)+Ke(2))*h(2)-Kp(2)*h(3)-(Kp(2)*dz-Kp(1)*dz);
y(3)=-Kp(2)*h(2)+(Kp(2)+Kp(3)+Ke(3))*h(3)-Kp(3)*h(4)-(Kp(3)*dz-Kp(2)*dz);
y(4)=-Kp(3)*h(3)+(Kp(3)+Kp(4)+Ke(4))*h(4)-Kp(4)*h(5)-(Kp(4)*dz-Kp(3)*dz);
y(5)=-Kp(4)*h(4)+(Kp(4)+Kp(5)+Ke(5))*h(5)-Kp(5)*h(6)-(Kp(5)*dz-Kp(4)*dz);
y(6)=-Kp(5)*h(5)+(Kp(5)+Kp(6)+Ke(6))*h(6)-Kp(6)*h(7)-(Kp(6)*dz-Kp(5)*dz);
y(7)=-Kp(6)*h(6)+(Kp(6)+Kp(7)+Ke(7))*h(7)-Kp(7)*h(8)-(Kp(7)*dz-Kp(6)*dz);
y(8)=-Kp(7)*h(7)+(Kp(7)+Kp(8)+Ke(8))*h(8)-Kp(8)*h(9)-(Kp(8)*dz-Kp(7)*dz);
y(9)=-Kp(8)*h(8)+(Kp(8)+Kp(9)+Ke(9))*h(9)-Kp(9)*h(10)-(Kp(9)*dz-Kp(8)*dz);
y(10)=-Kp(9)*h(9)+(Kp(9)+Kp(10)+Ke(10))*h(10)-Kp(10)*h(11)-(Kp(10)*dz-Kp(9)*dz);
y(11)=-Kp(10)*h(10)+(Kp(10)+Kp(11)+Ke(11))*h(11)-Kp(11)*h(12)-(Kp(11)*dz-Kp(10)*dz);
y(12)=-Kp(11)*h(11)+(Kp(11)+Kp(12)+Ke(12))*h(12)-Kp(12)*h(13)-(Kp(12)*dz-Kp(11)*dz);
y(13)=-Kp(12)*h(12)+(Kp(12)+Kp(13)+Ke(13))*h(13)-Kp(13)*h(14)-(Kp(13)*dz-Kp(12)*dz);
y(14)=-Kp(13)*h(13)+(Kp(13)+Kp(14)+Ke(14))*h(14)-Kp(14)*h(15)-(Kp(14)*dz-Kp(13)*dz);
J=jacobian(y,q);
dq=-J\f;
q=q+dq;
end