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
end
end
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
end
S
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
break
end
aold = a; %set aold to a
end
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
end
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);
end
end```