Matlab之Gauss-Newton法

haidaowang

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

  q) i: ~) A0 w8 S  x% E4 m

代码如下：         

• 
  9 \# i, k; f" N% l( J0 f7 r1 y
  8 B! u. [% S+ A# i# r: j3 e% using Gauss-Newton method
  . _/ }6 U, H; A5 }: Z
  / z& V2 d- L, C% C8 I9 s
• 
  
  8 V. q  U* L+ K" U% initial point(0,0)
  
  
• 
  
  % 1、with circul point (0,1),(1,1),(0,-1)
  , Z) [3 K% @4 T% g5 a
  
• 
  9 Q: N- p1 B6 s; x: [5 r
  % 2、with circul point (-1,0),(1,1),(1,-1)
  ; }7 i5 Y) o' T
  $ \, ~/ a5 d1 `3 H; D
• 
  
  ' Z2 i% s6 V9 S* ^% its radius is 1
  7 q( E: w$ q# H% h) i5 f
  ! g* ?; h- f1 K  R
• 
  
  " S6 Q! u" E' G% ], \9 i; w% output :the sum of distance is the minist
  
  : x  ?4 {0 K8 V3 [  \
• 
  % W. Y7 o$ I% ^
  . j- @2 a) K+ x  u5 lfunction page_237_1
  : q) o, M( Y4 x2 K( ~$ ?
  
• 
  6 M9 o8 a8 l% H) N- ^4 T- X
  format long
  0 t+ w% N6 N. R7 s- d
  ( O% ~) Z4 W- [2 b+ A
• 
  
  3 Z/ T9 b' p2 E9 a' y8 x/ Vsyms x y
  + ~" e0 c. @# y7 R5 w
  ( S, W$ e6 x# d& N) j
• 
  
  r1      = sqrt(x^2+(y-1)^2) - 1;
  
  
• 
  ( p' Z( R& S4 J( U
  5 r6 e1 q& R0 H& Lr2      = sqrt((x-1)^2+(y-1)^2) - 1;
  + h1 F( |- `; o8 F- ?( F/ y
  ; q2 ~: r" p% n) p8 w
• 
  9 A1 y+ \, T7 i- m" S; F7 o
  , q# k: ~) P5 |4 G. b- }r3      = sqrt(x^2+(y+1)^2) - 1;
  
  & c: o, X% i- Q; ~
• 
  1 S( p( z8 u: ^) _
  # F# L1 ~/ `* B7 Orxk     = [r1;r2;r3];
  ( s- }0 F) Y9 F% f) u
  8 H0 Z& r1 Y7 t% a
• 
  0 x' y$ ~" X( [
  5 Z# u" X$ B! @. r- hDrx     = jacobian([r1;r2;r3]);
  
  
• 
  
  1 @4 J) T4 @0 b" p/ ^9 I% q: {; oDrxT    = Drx';
  ( K8 K: }; o! X
  
• 
  1 ?6 i) W9 \+ C) _  m
  9 r  O$ w5 f) U! }, ADrxTdrx = DrxT*Drx;                 %Dr(xk)TDr(xk)
  6 M* y8 Q# d" }0 D
  
• 
  
  , t$ @: w2 g  ]* h1 f1 VDrxTrxk = -1*DrxT*rxk;              %Dr(xk)Tr(xk)
  6 v3 `4 W! z# N
  
• 
  
  & i, m/ t4 H' `, y5 o# e%vk     = DrxTrxk\DrxTdrx;
  
  
• 
  9 s  F/ P/ E  A# v4 l- @/ H9 |
  8 r2 R( ], h, s7 r  |+ ADrxTdrx_= subs(DrxTdrx,{x,y},{0,0});
  
  
• 
  
  DrxTrxk_= subs(DrxTrxk,{x,y},{0,0});
  
  # }+ _$ ]& U: T7 R1 O# }5 a! s
• 
  
  xy      = zeros(2,1);
  8 V4 R4 R% A( l6 o; g! T% c* n: e
  
• 
  
  8 R# [1 R/ l6 o9 ^7 pfor k = 1:15
  
  
• 
  
  ' T7 n+ M- G3 V" b1 w' I    v       = DrxTdrx_\DrxTrxk_;
  
  
• 
  ' y9 b. _% y9 O+ I" O/ A( W$ v. U: ^
  ; v8 s9 @2 L; E$ o' N    xy      = xy+v;      
  
  0 j/ _9 ?8 C$ C$ Z
• 
  : T2 [+ x5 f9 Y
  0 w2 k, ?# U9 G% j6 v' K    DrxTdrx_= subs(DrxTdrx,{x,y},{xy(1),xy(2)});
  * ?' w, u6 z' [# t  D
  - C3 j- i: x2 l5 p, A- Q* D
• 
  
  3 {( V4 {  {% s5 F    DrxTrxk_= subs(DrxTrxk,{x,y},{xy(1),xy(2)});
  - v0 |- w$ ^, E; M
  
• 
  / x! O! q4 F0 Z
  + ]) ~& v% z. T; f. oend
  ' @3 u: o! q* N  D
    j, b. n  [: R4 \' S" R& d
• 
  - Y# H6 Q; @8 O7 f
  disp('a with circular initial vector(0,0),Centers(0,1), (1,1), (0,?1)and all radii 1,the point is');
  ( E  U. L* e$ @! w
  1 D7 z: r# t4 O  L) }
• 
  * B  z) X  S& ~: S
  vpa(xy,6)
  + F3 a, U) v$ b% _
  : R) c5 D8 r: j* l3 x2 Z  d' E
• 
  % I* D2 T$ E3 Q# m& G$ E3 j
  ! t( I7 m* N& L  }) j* r# o

& F" }6 d7 m# E0 E; P

NingW

Matlab之Gauss-Newton法