Matlab之Gauss-Newton法

haidaowang

这一题就是已知初始点，然后求某一点到三个圆的距离的和的最小值。使用的就是大名鼎鼎的Gauss-Newton方法。这里要求雅克比行列式，由函数 jacobian可以求出。所以这就减少了许多不必要的麻烦。
3 R) r8 W9 k% o0 a4 i4 z

代码如下：         
6 m  c# F0 U8 a* \6 R+ V4 {

2 c' c$ M: `/ ]* j+ I4 {

• 
  9 h7 X" S" P- a- j
  % using Gauss-Newton method
  
  
• 
  
  % initial point(0,0)
  , F: g% H  B7 o4 h
  0 ]: m& Y$ y7 w6 h; j' \
• 
  2 g1 V5 C& n) M% N$ m0 Q8 b5 d
  . d$ k+ t* r3 K. _+ R3 T  N7 P% 1、with circul point (0,1),(1,1),(0,-1)
  
  
• 
  
  0 A# {# A( U" T" O% 2、with circul point (-1,0),(1,1),(1,-1)
  
  
• 
  ) w4 e1 }; e5 N, b+ a
  + P+ [6 m% v, E3 Q% its radius is 1
    a. h1 g4 H3 v
  
• 
    h; U/ S& e9 l$ E
  / S( }0 c' R/ A+ ]5 ]/ c. H4 e% output :the sum of distance is the minist
  
  
• 
  
  . L& Q3 r, d1 `& H7 d2 A/ nfunction page_237_1
  : D9 f3 w5 \8 b5 @$ U7 _6 ?
  
• 
  
  8 }, m5 w4 `  Q/ Tformat long
  
  1 k$ c' y- c2 Z2 l2 F% m: Q
• 
  1 y6 T* y; A4 q) K) J" l. A
  + a- Y- n% Z+ \7 a* Fsyms x y
  / S! u" E9 Z/ E  ]. ?. Q; R" B' ]$ {
  
• 
  ' V& ~* i( J2 u: q! Y# @
  ' H; [, i* h* b& u3 C5 Y9 Rr1      = sqrt(x^2+(y-1)^2) - 1;
  
  
• 
  
  , J% G4 @$ O1 I3 I7 |/ xr2      = sqrt((x-1)^2+(y-1)^2) - 1;
  ; e- ?/ g+ h# f& n& A* ]: N
  8 d" P' Z6 m: E# V
• 
  
  9 C* D) k" X( j  K. P  Y8 s( a, Ar3      = sqrt(x^2+(y+1)^2) - 1;
  
  ! e$ Q3 [- [& ?, {9 ^
• 
  
  ' Z( y2 G; H% N, g3 p% orxk     = [r1;r2;r3];
  1 U, t1 f, ~! A  K: Z( y
  ; Y" O: A" j& `  \0 m2 A9 U
• 
  
  Drx     = jacobian([r1;r2;r3]);
  2 h' D0 r6 i9 G3 o8 {7 j
  ' D& T: b/ Q0 `& C( L9 ?
• 
  ( z8 X5 \  B+ `
  , M) i9 v( I$ q; ]6 l- @DrxT    = Drx';
  
  + S( N# a7 L( {2 l1 q$ g# e  d" y
• 
  
  * J: c3 M7 X1 l0 A  n/ ~DrxTdrx = DrxT*Drx;                 %Dr(xk)TDr(xk)
  
  ! E: h- o+ J1 m: W+ \
• 
  1 v$ l* W9 _% l+ b' g3 I& d- }' I1 w/ h
  DrxTrxk = -1*DrxT*rxk;              %Dr(xk)Tr(xk)
  
  8 U8 A$ h/ v6 ~" [) ^9 y1 W5 r
• 
  * D( q) |1 b; j9 ^
  %vk     = DrxTrxk\DrxTdrx;
  
  8 w8 R: b( }, h( J2 ~5 T
• 
  
  8 p! R& c7 p: U# {* K" d" q3 WDrxTdrx_= subs(DrxTdrx,{x,y},{0,0});
  , P. F: v- D8 ?+ R
  , q6 w4 P3 J6 o" E5 y5 H2 b
• 
  
  / ~# g# j6 ~& I6 a( xDrxTrxk_= subs(DrxTrxk,{x,y},{0,0});
  
  . T1 v' E: @1 k! d3 b- P7 K1 e
• 
  ! i% P: R. z( {2 [( U; a
  ) `* \% ]( ?! V6 X9 z3 z  ~xy      = zeros(2,1);
  ' b0 ~# ^8 G0 w. h- G( W
  
• 
  7 I6 g" `% Z9 _7 q6 X$ h8 O
  for k = 1:15
  
  
• 
  
  ! H  ?+ |+ B* D! R0 |    v       = DrxTdrx_\DrxTrxk_;
    ~3 e( j7 R  [
  ) }3 ^( D, z: x" @# a* b
• 
  
      xy      = xy+v;      
  * D/ J) ]1 |# ]6 a) m3 x& G
  
• 
  
      DrxTdrx_= subs(DrxTdrx,{x,y},{xy(1),xy(2)});
    X7 X* ?2 p2 H0 z
  " t. I; P; Y* H
• 
  
  2 I! g$ Q" M2 _" {; t9 o5 p1 o    DrxTrxk_= subs(DrxTrxk,{x,y},{xy(1),xy(2)});
  # Z6 G. W! m7 Y7 M% K& E
  4 y6 @& \6 h+ |  K" J
• 
  
  , M. g. `! ~) N) v3 b3 Iend
  
  
• 
  ' m: N7 s2 ?# X" d% }2 y& z2 V8 f/ d
  disp('a with circular initial vector(0,0),Centers(0,1), (1,1), (0,?1)and all radii 1,the point is');
  9 F( H9 g) i; j
  
• 
  
  vpa(xy,6)
  
  
• 
  
  

NingW

Matlab之Gauss-Newton法