%
% -----------------------------------------------------
%
%   Approximation  im  Wintersemester 2013/14
%   H.J. Oberle
%
%   L_2 Approximation ebener Daten durch einen Kreis                
%   H. Spaeth: Computing 57, 1790-185, 1996
%
%   Aufgabe 26 i), ii)                 
%
% -----------------------------------------------------
%
clear all
close all
%
beispiel = 1;
if beispiel == 1
%    
%  Daten i)
    X = [5 10 7 9 2 0 -1 -1];
    Y = [-3 3 5 9 8 3 -1 7];
else
%    
%  Daten ii)
    X = [1 2 5 7 9 3];
    Y = [7 6 8 7 5 7];
end    
%
n1 = length(X);
n  = n1-1;
%
%   Startdaten (nach Spaeth)
%
a(1,1) = 2*sum(X.*X);
a(1,2) = 2*sum(X.*Y);
a(1,3) = 2*sum(X);
a(2,1) = a(1,2);
a(2,2) = 2*sum(Y.*Y);
a(2,3) = 2*sum(Y);
a(3,1) = a(1,3);
a(3,2) = a(2,3);
a(3,3) = 2*n1;
%
c(1) = sum(X.*(X.*X+Y.*Y));
c(2) = sum(Y.*(X.*X+Y.*Y));
c(3) = sum(X.*X+Y.*Y);
%
c  = c';                                                      %'
d  = a\c;
xm = [d(1); d(2)]; 
xn = xm;
%
%  Newton-Iteration
%  iteil=1:    Funktionsauswertung,
%  iteil=2,3:  Numerische Differentiation 
%
iter  = 0;
konv  = 0;
itmax = 20;
%
format long e
format compact
fprintf(' Newton-Iteration: \n')
fprintf(' \n')
fprintf(' iter     x(1)                x(2)                r                ||f|| \n')
fprintf(' \n')
%
while and(konv == 0,iter < itmax)
%
for iteil = 1:3
%
sn = (Y-xm(2))./(sqrt((X-xm(1)).^2+(Y-xm(2)).^2));
cn = (X-xm(1))./(sqrt((X-xm(1)).^2+(Y-xm(2)).^2));
%
e1 = sum(X);
e2 = sum(Y);
e3 = sum(X.*cn+Y.*sn);
u1 = sum(cn);
u2 = sum(sn);
r1 = (n1*e3-u1*e1-u2*e2)/(n1*n1 - u1*u1 - u2*u2);
xn(1) = (e1-r1*u1)/n1;
xn(2) = (e2-r1*u2)/n1; 
%
if iteil == 1
   x  = xm;
   x0 = xn;
   f  = xm - xn; 
%
fprintf(' %2.0f   %1.10e  %1.10e  %1.10e  %1.5e\n',...
iter, x(1), x(2), r1, norm(f)) 
%
   delta = abs(xm(1))*1e-8;
   if delta < 1e-8  
   delta = 1e-8;
   end
   xm(1) = xm(1) + delta;
end % ----------------------    iteil=1
%
if iteil == 2
   J(1,1) = 1 - (xn(1) - x0(1))/delta;
   J(2,1) = -(xn(2) - x0(2))/delta;   
   xm = x;
   delta = abs(xm(2))*1e-8;
   if delta < 1e-8  
   delta = 1e-8;
   end
   xm(2) = xm(2) + delta;
end % -----------------------   iteil=2
%
if iteil == 3
   J(1,2) = -(xn(1) - x0(1))/delta;
   J(2,2) = 1 - (xn(2) - x0(2))/delta;   
%
   dx = J\f;
   xm = x - dx; 
   iter = iter + 1;
   if norm(dx) < 1e-10
   konv = 1;
   end 
%
end % -----------------------  iteil=3
%
end
%
end
%
%    Berechnung des Residuums
%
res = norm((X - (x(1)+r1*cn)).^2 + (Y - (x(2)+r1*sn)).^2);
fprintf(' \n')
fprintf(' Residuum:  %1.5e\n',res) 
%
t  = [0:0.01:6.4];
x1 = x(1) + r1*cos(t);
y1 = x(2) + r1*sin(t);
plot(x1,y1,'b','Linewidth',1.2)
hold
%
if beispiel == 1
   xlabel('Aufgabe 26 i)','FontSize',14)
   axis([-4 11 -4 11])
   axis('square')
   set(gca,'FontSize',12)
   plot(x(1),x(2),'+','Linewidth',1.2)
   plot(X,Y,'k.','Markersize',20)
   print -depsc aufg26i
%
else
   xlabel('Aufgabe 26 ii)','FontSize',14)
   axis([-1 10 -2 9])
   axis('square')
   set(gca,'FontSize',12)
   plot(x(1),x(2),'+','Linewidth',1.2)
   plot(X,Y,'k.','Markersize',20)
   print -depsc aufg26ii
end
%
    



