nanewton1.m

牛顿迭代法 若高阶非线性方程组: u ( x , y) = 0 v ( x , y) = 0 可以用迭代公式

function  [x,y] = nanewton1( xd,yd )
%NEWTON Summary of this function goes here
%  Detailed explanation goes here

%x0 为初始给定的一点 ，e 为迭代精度  N为迭代上限
% f为f( x, y )函数，dfx为f( x, y )对x所求的导函数，dfy为f( x, y )对y所求的导函数，
% g为g( x, y )函数，dgx为g( x, y )对x所求的导函数，dfy为g( x, y )对y所求的导函数，
% k1 = 0.13 ;
% k2 = 0.11 ;
% xd = 0.45 ;
% yd = 0.39 ;
err = 0.00000001;
k1 = 0.080659002 ;
k2 = 0.0059059998 ;

x0 = xd ;
y0 = yd ;


x = x0;
y = y0;

cha =1;

while ( cha > err )
r = sqrt ( x^2 + y^2 );
k =  1 + k1 * r^2 + k2 * r^4 ;
f = ( 1 + k1 * r^2 + k2 * r^4 ) * x - x0 ;
g = ( 1 + k1 * r^2 + k2 * r^4 ) * y - y0 ;

dfx = 1 + k1 * r^2 + k2 * r^4 + 2 * k1 * x^2 + 4 * k2 * x^2 *r^2;
dfy = 2 * k1 * x * y / r + 4 * k2 * x * y * r^2;
dgx = 2 * k1 * x * y / r + 4 * k2 * x * y * r^2;
dgy = 1 + k1 * r^2 + k2 * r^4 + 2 * k1 * y^2 + 4 * k2 * y^2 * r^2;

Jacobi = det ( [dfx , dfy ; dgx , dgy] );

  if Jacobi == 0   
     continue
  else
  Jx = det ( [ dfy , f ; dgy , g] ) / Jacobi;
  Jy = det ( [ f , dfx ; g ,dgx ] ) / Jacobi;

  xn = x + Jx ;
  yn = y + Jy ;
 
  cha =abs(max( xn - x , yn - y ));

  end 
x = xn ;
y = yn ;
end
x = xn ;
y = yn ;

f = ( 1 + k1 * r^2 + k2 * r^4 ) * x - x0 ;
g = ( 1 + k1 * r^2 + k2 * r^4 ) * y - y0 ;
%figure, imshow ( x,y) ;

% 
% if nargin <5,N=500;end
% if nargin <4,e=1e-4;end
% x=x0;x0=x+2*e;k=0;
% while abs(x0-x)>e&k<N,
%     k=k+1;
%     x0=x;x=x0-feval(fname,x0)/feval(dfname,x0);
%     disp(x)
% end
% if k==N,warning('已经达到跌代次数上限');end
% gradient
% diag
% det