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

lianghuan

我写了一个关于牛拉法解非线性方程的程序，进行故障测距的。但结果总是发散，哪位同行前辈或高手帮忙看看，到底哪出错了呀？
万分感激
' `$ e3 s% `3 ?8 M- p! h/ u0 I# o
8 a6 }; s, T! Z9 c9 G0 ofunction nlf2 =nlf2(Um1,Um2,Un1,Un2,Im1,Im2,In1,In2,r1,Zc1)
3 U  h  g9 H0 q- t; l7 @x=150%迭代过程没问题
L=300;
" k+ s7 x% I8 G( ukk=0;
e=1;
! \2 H8 ]) Z1 X2 f8 }# O: P%使用牛顿-拉夫逊法迭代求解
while (kk<10000&&e>0.001)
: F- k& K5 Q% j! |%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)));
. O. T. J, C3 h- C1 o7 `( A4 L   %对A展开整理得(sinh(x))'=cosh(x) 　　(cosh(x))'=sinh(x)
5 ]0 Z, |, _) a% M  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)))...
3 v2 V4 A7 E1 r* C) s   +(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);
    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))...
1 o; y6 f6 s- `1 ~; j9 U) O    +(Um2*In1*Zc1-Um1*In2*Zc1)*(r1*sinh(r1*x)*sinh(r1*L-r1*x)-r1*cosh(r1*x)*cosh(r1*L-r1*x))...
, A* |) @3 p. K7 v/ n6 }& n    +(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))
' b" {) g5 P/ ^* D4 _1 V       a11=real(B1);
       a21=imag(B1);
: Y7 b' V0 s0 {1 s  g" A  %A对L求偏导数
" O2 V/ M. [( c3 v; P" x  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))...
1 r- ^' F: h; p/ k$ a0 C: r      +(Im2*Un1*Zc1-Un2*Im1*Zc1)*(r1*sinh(r1*x)*sinh(r1*L-r1*x))...
' {1 H( b4 D/ _1 {      +(Im1*In2*(Zc1.^2)-Im2*In1*(Zc1.^2))*(r1*cosh(r1*L-r1*x)*sinh(r1*x))
      a12 = real(B2);
* [, m6 U1 X( l3 y/ O      a22 = imag(B2);
: Z8 M7 {* h+ R
     R = [a11, a12; a21, a22]
    RR=inv(R);
    detb=inv(R) * f
    %detb =R\f;
1 ^3 Y, l+ D- w$ l$ e- L" a     L=L-detb(1);
     x=x-detb(2);
; {8 M+ c; Y# v' ^- y& c4 n% U     e1=abs(detb(1));
7 `' h* x- W# {; e3 W     e2=abs(detb(2));
) z$ N, V+ y' o7 T2 T: }3 a     e=max(e1,e2)
3 o3 L! j: f# b: C* t     kk = kk + 1
end
$ o4 M! z" `% P0 R5 O6 E& Onlf2 = [L,x]