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

lianghuan

我写了一个关于牛拉法解非线性方程的程序，进行故障测距的。但结果总是发散，哪位同行前辈或高手帮忙看看，到底哪出错了呀？
万分感激

+ G6 C9 k% }6 k: e8 d$ [function nlf2 =nlf2(Um1,Um2,Un1,Un2,Im1,Im2,In1,In2,r1,Zc1)
x=150%迭代过程没问题
L=300;
4 E  C7 X1 L: T% ~; k  ]$ pkk=0;
: i$ T* m, Q; v2 ^6 i7 J, ?5 fe=1;
%使用牛顿-拉夫逊法迭代求解
while (kk<10000&&e>0.001)
: `0 ~7 N$ W! @%A=(Um1*cosh(r1*x)-Im1*Zc1*sinh(r1*x))*(Un2*cosh(r1*(L-x))-In2*Zc1*sinh(r1*(L-x)))...
%  -(Um2*cosh(r1*x)-Im2*Zc1*sinh(r1*x))*(Un1*cosh(r1*(L-x))-In1*Zc1*sinh(r1*(L-x)));
   %对A展开整理得(sinh(x))'=cosh(x) 　　(cosh(x))'=sinh(x)
  D=(Um1*Un2-Um2*Un1)*(cosh(r1*x)*cosh(r1*(L-x)))...
* v+ v: ~* E7 Z   +(Um2*In1*Zc1-Um1*In2*Zc1)*(cosh(r1*x)*sinh(r1*(L-x)))...
   +(Im2*Un1*Zc1-Un2*Im1*Zc1)*(sinh(r1*x)*cosh(r1*(L-x)))...
( o' r4 @  E7 x( {$ I" C3 p8 P   +(Im1*In2*(Zc1.^2)-Im2*In1*(Zc1.^2))*(sinh(r1*x)*sinh(r1*(L-x)))
, o- @9 R" M. }  c5 Y" u    f1 = real(D);
    f2 = imag(D);
) o$ f, @; f5 q8 @$ O( P! G. \    f = [f1,f2].'
  %A对x求偏导数
  }, \: I" a" o- g% f  B1=(Um1*Un2-Um2*Un1)*(r1*sinh(r1*x)*cosh(r1*L-r1*x)-r1*cosh(r1*x)*sinh(r1*L-r1*x))...
. ]. N) V0 c: @' S    +(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))...
" W6 t6 b: b7 ~* Z    +(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))
  v0 k( n" o2 U+ C& w% T       a11=real(B1);
       a21=imag(B1);
9 ?) o) U9 f0 F% B, J# b  %A对L求偏导数
% h6 |8 o1 |3 x2 X  B2=(Um1*Un2-Um2*Un1)*(r1*cosh(r1*x)*sinh(r1*L-r1*x))...
6 m3 y, e' Y% r' g2 l  S& Y0 v      +(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))...
      +(Im1*In2*(Zc1.^2)-Im2*In1*(Zc1.^2))*(r1*cosh(r1*L-r1*x)*sinh(r1*x))
4 z: M0 g9 {3 D8 k      a12 = real(B2);
( U! M$ n+ w5 n3 J8 L      a22 = imag(B2);
( D5 W  J3 ~' V! ^; l+ V5 A. _
     R = [a11, a12; a21, a22]
    RR=inv(R);
5 n; o! S8 W0 c    detb=inv(R) * f
& Q+ j4 [! l- z/ k- r    %detb =R\f;
     L=L-detb(1);
2 R8 ?# U  N; ]' U: W     x=x-detb(2);
7 n1 f4 c0 W& w% q. m- I, W     e1=abs(detb(1));
$ {+ g) A' u3 H     e2=abs(detb(2));
     e=max(e1,e2)
6 Q8 ^2 N/ G2 W! i4 w     kk = kk + 1
end
. F9 H7 y0 ^, X- mnlf2 = [L,x]