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

lianghuan

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

function nlf2 =nlf2(Um1,Um2,Un1,Un2,Im1,Im2,In1,In2,r1,Zc1)
" H( [/ Y/ X$ r1 m7 u+ J5 Tx=150%迭代过程没问题
L=300;
8 b( ^- o6 S9 f$ u8 E* Gkk=0;
2 o# S0 d; D& A/ |2 ]" Ue=1;
! o; {0 u) V2 J7 j# ^6 i%使用牛顿-拉夫逊法迭代求解
while (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)))...
%  -(Um2*cosh(r1*x)-Im2*Zc1*sinh(r1*x))*(Un1*cosh(r1*(L-x))-In1*Zc1*sinh(r1*(L-x)));
7 c! i& q: K: C/ C* D   %对A展开整理得(sinh(x))'=cosh(x) 　　(cosh(x))'=sinh(x)
  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)))...
: V, e7 L( |! g& C; y. t$ N8 u5 v- f   +(Im2*Un1*Zc1-Un2*Im1*Zc1)*(sinh(r1*x)*cosh(r1*(L-x)))...
0 U; ~/ A6 d* x, D4 @% L) b   +(Im1*In2*(Zc1.^2)-Im2*In1*(Zc1.^2))*(sinh(r1*x)*sinh(r1*(L-x)))
    f1 = real(D);
, o; L. W* M3 Z$ `0 w    f2 = imag(D);
    f = [f1,f2].'
" D0 [/ O! s& @3 O; q( i" Y, L8 I  %A对x求偏导数
  B1=(Um1*Un2-Um2*Un1)*(r1*sinh(r1*x)*cosh(r1*L-r1*x)-r1*cosh(r1*x)*sinh(r1*L-r1*x))...
$ d$ Q7 \6 b* I& G0 g    +(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))...
    +(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))
; Z+ s# j: l( k% ^) w, |       a11=real(B1);
- D. a$ i+ M& F, K& _7 g       a21=imag(B1);
  %A对L求偏导数
, s/ R# ]" Q+ ?* L. ]( `  B2=(Um1*Un2-Um2*Un1)*(r1*cosh(r1*x)*sinh(r1*L-r1*x))...
! \) N% J/ ~: m  K      +(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))...
0 \& W) q3 Q) ?* a+ v4 q8 g      +(Im1*In2*(Zc1.^2)-Im2*In1*(Zc1.^2))*(r1*cosh(r1*L-r1*x)*sinh(r1*x))
      a12 = real(B2);
. j( k/ F( S7 Z  S      a22 = imag(B2);
- J, v( \4 Z$ E7 E2 V8 H
     R = [a11, a12; a21, a22]
% o+ S' x1 g) w. J0 T- k    RR=inv(R);
3 @; o/ L; a  W3 W- a    detb=inv(R) * f
2 T0 `/ w$ R  r4 M2 }4 c    %detb =R\f;
     L=L-detb(1);
" u. q" @+ ^1 `! }0 t& L     x=x-detb(2);
' l! B+ o, y% _/ ^$ _     e1=abs(detb(1));
% V4 b! u5 z5 y( ?4 n     e2=abs(detb(2));
     e=max(e1,e2)
" s+ }. ^+ g+ v8 l+ `. b     kk = kk + 1
end
nlf2 = [L,x]