www.pudn.com > GaussNewton.rar > GaussNewton.m, change:2015-10-19,size:2265b

function [unknowns,steps,S] = GaussNewton() 
%GaussNewton- uses the Gauss-Newton method to perform a non-linear least 
%squares approximation for the origin of a circle of points 
% Takes as input a row vector of x and y values, and a column vector of 
% initial guesses. partial derivatives for the jacobian must entered below 
% in the df function 
format long 
tol = 1e-8; %set a value for the accuracy 
maxstep = 30; %set maximum number of steps to run for 
z = 2*pi*rand(100,1); %create 100 random points 
p = cos(z); %create a circle with those points 
q = sin(z); 
d = ones(100,1); %set distance to origin as 1 for all points 
a = [0.1;0.1]; %set initial guess for origin 
m=length(p); %determine number of functions 
n=length(a); %determine number of unkowns 
aold = a; 
for k=1:maxstep %iterate through process 
    S = 0; 
    for i=1:m 
         for j=1:n 
             J(i,j) = df(p(i),q(i),a(1,1),a(2,1),j); %calculate Jacobian 
             JT(j,i) = J(i,j); %and its trnaspose 
         Jz = -JT*J; %multiply Jacobian and 
 %negative transpose 
    for i=1:m 
        r(i,1) = d(i) - sqrt((a(1,1)-p(i))^2+(a(2,1)-q(i))^2); %calculate r 
        S = S + r(i,1)^2; %calculate sum of squares of residuals 
   g = Jz\JT; %mulitply Jz inverse by J transpose 
   a = aold-g*r; %calculate new approximation 
   unknowns = a; 
   err(k) = a(1,1)-aold(1,1); %calculate error 
   if (abs(err(k)) <= tol); %if less than tolerance break 
      aold = a; %set aold to a 
steps = k; 
hold all 
plot(p,q, 'r*' ) %plot the data points 
plot(a(1,1),a(2,1), 'b*' ) %plot the approximation of the origin(expect it to be 0,0) 
title('Gauss-Newton Approximation of Origin of Circular Data Points' ) %set axis lables, title and legend 
xlabel('X' ) 
ylabel('Y' ) 
legend('Data Points' , 'Gauss-Newton Approximation of Origin' ) 
hold off 
errratio(3:k) = err(2:k-1)./err(3:k); %determine rate of convergence 
function value = df(p,q,a1,a2,index) %calculate partials 
         switch index 
                case 1 
                     value = (2*a1 - 2*p)*0.5*((a1-p)^2+(a2-q)^2)^(-0.5); 
                case 2 
                    value = (2*a2 - 2*q)*0.5*((a1-p)^2+(a2-q)^2)^(-0.5);