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

lianghuan

我写了一个关于牛拉法解非线性方程的程序，进行故障测距的。但结果总是发散，哪位同行前辈或高手帮忙看看，到底哪出错了呀？
万分感激
1 Y6 x0 o/ j$ I
function nlf2 =nlf2(Um1,Um2,Un1,Un2,Im1,Im2,In1,In2,r1,Zc1)
x=150%迭代过程没问题
L=300;
7 N# C1 f7 G' X" m) R* N0 r1 G; r: E) ]kk=0;
e=1;
' A) D/ G3 U  f%使用牛顿-拉夫逊法迭代求解
while (kk<10000&&e>0.001)
1 b; U$ G. e3 O8 H%A=(Um1*cosh(r1*x)-Im1*Zc1*sinh(r1*x))*(Un2*cosh(r1*(L-x))-In2*Zc1*sinh(r1*(L-x)))...
- Y8 N( X7 b0 D: k: R! p! R%  -(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)
7 r6 N, n8 e1 J3 Z! Y  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)))...
) `& Y5 B4 [& ^: X   +(Im2*Un1*Zc1-Un2*Im1*Zc1)*(sinh(r1*x)*cosh(r1*(L-x)))...
   +(Im1*In2*(Zc1.^2)-Im2*In1*(Zc1.^2))*(sinh(r1*x)*sinh(r1*(L-x)))
    f1 = real(D);
    f2 = imag(D);
2 |3 a+ e# S1 C# [" U0 d    f = [f1,f2].'
  %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))...
    +(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))...
! k0 ^4 E6 B) x8 ?    +(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))
0 G! Q$ Z; G& o1 [       a11=real(B1);
; L" V3 Y  X! i" N# Y       a21=imag(B1);
8 s* i! y8 W* J( z, f9 _  %A对L求偏导数
  B2=(Um1*Un2-Um2*Un1)*(r1*cosh(r1*x)*sinh(r1*L-r1*x))...
      +(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))...
! I* f% ^+ L9 z      +(Im1*In2*(Zc1.^2)-Im2*In1*(Zc1.^2))*(r1*cosh(r1*L-r1*x)*sinh(r1*x))
; K2 N# [+ Z& u7 i- c      a12 = real(B2);
      a22 = imag(B2);

     R = [a11, a12; a21, a22]
' `; P: P. U: ]. O1 L  u$ a) \# [    RR=inv(R);
    detb=inv(R) * f
: a/ J2 |) E' w* U! N8 V    %detb =R\f;
     L=L-detb(1);
" n; t* R0 w  g1 |: @4 ?1 x% K     x=x-detb(2);
     e1=abs(detb(1));
     e2=abs(detb(2));
) h. }2 w+ t4 `- ^% ^     e=max(e1,e2)
     kk = kk + 1
9 {& i3 J5 z! \% S( x9 Kend
# m5 [& W2 g5 Q0 B. x: vnlf2 = [L,x]