牛拉法解非线性方程的程序

lianghuan

我写了一个关于牛拉法解非线性方程的程序，进行故障测距的。但结果总是发散，哪位同行前辈或高手帮忙看看，到底哪出错了呀？
万分感激
% l8 ^9 ^' U! `3 W
function nlf2 =nlf2(Um1,Um2,Un1,Un2,Im1,Im2,In1,In2,r1,Zc1)
x=150%迭代过程没问题
' x9 D5 ]/ v& [( jL=300;
kk=0;
" K3 \' _* U. e0 U% P4 Ze=1;
# r1 s  j0 Q  w, N9 Z, _) U1 `" ?0 K%使用牛顿-拉夫逊法迭代求解
) a- r" E% Z1 Ewhile (kk<10000&&e>0.001)
%A=(Um1*cosh(r1*x)-Im1*Zc1*sinh(r1*x))*(Un2*cosh(r1*(L-x))-In2*Zc1*sinh(r1*(L-x)))...
/ j+ k# [# O1 @/ K& J1 X- I%  -(Um2*cosh(r1*x)-Im2*Zc1*sinh(r1*x))*(Un1*cosh(r1*(L-x))-In1*Zc1*sinh(r1*(L-x)));
- s/ C7 i# ^$ v) c- W   %对A展开整理得(sinh(x))'=cosh(x) 　　(cosh(x))'=sinh(x)
: Y6 t( S2 |8 G  D=(Um1*Un2-Um2*Un1)*(cosh(r1*x)*cosh(r1*(L-x)))...
   +(Um2*In1*Zc1-Um1*In2*Zc1)*(cosh(r1*x)*sinh(r1*(L-x)))...
+ W# V% w" Q: ]7 W. i/ D6 @" N   +(Im2*Un1*Zc1-Un2*Im1*Zc1)*(sinh(r1*x)*cosh(r1*(L-x)))...
8 T7 y& R6 z) B) y3 L# @   +(Im1*In2*(Zc1.^2)-Im2*In1*(Zc1.^2))*(sinh(r1*x)*sinh(r1*(L-x)))
    f1 = real(D);
. _5 O: I) y* b- t+ ~( _    f2 = imag(D);
    f = [f1,f2].'
  %A对x求偏导数
3 w- o1 P# n, U" O1 Y* X/ N  B1=(Um1*Un2-Um2*Un1)*(r1*sinh(r1*x)*cosh(r1*L-r1*x)-r1*cosh(r1*x)*sinh(r1*L-r1*x))...
    +(Um2*In1*Zc1-Um1*In2*Zc1)*(r1*sinh(r1*x)*sinh(r1*L-r1*x)-r1*cosh(r1*x)*cosh(r1*L-r1*x))...
    +(Im2*Un1*Zc1-Un2*Im1*Zc1)*(r1*cosh(r1*x)*cosh(r1*L-r1*x)-r1*sinh(r1*x)*sinh(r1*L-r1*x))...
, C' `! X0 a$ Y. m" h, `    +(Im1*In2*(Zc1.^2)-Im2*In1*(Zc1.^2))*(r1*cosh(r1*x)*sinh(r1*L-r1*x)-r1*cosh(r1*L-r1*x)*sinh(r1*x))
3 R: w0 M, i2 M8 n5 @5 Q       a11=real(B1);
       a21=imag(B1);
  %A对L求偏导数
. I4 F( g' ~, H2 [  B2=(Um1*Un2-Um2*Un1)*(r1*cosh(r1*x)*sinh(r1*L-r1*x))...
5 i9 _  Y) n# d6 Q      +(Um2*In1*Zc1-Um1*In2*Zc1)*(r1*cosh(r1*x)*cosh(r1*L-r1*x))...
      +(Im2*Un1*Zc1-Un2*Im1*Zc1)*(r1*sinh(r1*x)*sinh(r1*L-r1*x))...
9 u  B3 I* l2 P6 [6 Z      +(Im1*In2*(Zc1.^2)-Im2*In1*(Zc1.^2))*(r1*cosh(r1*L-r1*x)*sinh(r1*x))
( @# K" X* ~" r* f1 @      a12 = real(B2);
      a22 = imag(B2);

% T' s6 v$ q; M% Y6 s4 k     R = [a11, a12; a21, a22]
/ O# r% ?7 i; |* S    RR=inv(R);
2 J3 X$ c8 P$ I( `6 n% m0 }    detb=inv(R) * f
+ y, y4 [8 k) G" s- Y8 m$ n& l" ]1 q    %detb =R\f;
) f' Q/ G9 Z; c) D, O     L=L-detb(1);
# J3 R" K0 E/ ^2 U$ U  n! @     x=x-detb(2);
     e1=abs(detb(1));
/ k; g) w7 h" W; C% K# _     e2=abs(detb(2));
     e=max(e1,e2)
/ W1 `, y+ J) q: H. E) o4 Q4 I+ v     kk = kk + 1
# m$ L7 T1 e, r0 yend
nlf2 = [L,x]