%aufgabe 4 Serie 5 marco weber

function gn
clear all; 
format compact; 
format long;
%Set start value 


%für 4i)
%x0=[-0.5;-2.0;1.75;1.2;0.8]
%für 4ii)
x0=[-1;-1;1;1;1]



%Set Tolerence values 
rtol=1e-8; epsilon=1e-12;
%Without minimisation
withmini=0;
[iter,x,resid] = gaussnewton(x0,rtol,epsilon,withmini,@funnew)

%With minimisation
withmini=1;

[iter,x,resid] = gaussnewton(x0,rtol,epsilon,withmini,@funnew)


function [iter,x,resid] = gaussnewton(x0,rtol,epsilon,withmini,f)

% Gauss-Newton for nonlinear overdetermined system, 
% optimized for Numerische Mathematik, D-MAVT FS 08, K.Nipp, Serie 5,
% exercise 4
%
% in:   x0       starting values for iteration
%       rtol     stop if ||x-xold|| <= rtol*||x||
%       epsilon  stop minimization if F(y) < F(x)+epsilon
%       withmini do NO minimization if NOT set to 1
%       f        function containing the evaluation of f(x)-z (string)
%                and its Jacobian. One output argument gives
%                just the residual
% out:  iter     number of iterations needed till convergence
%       x        solution in least square sense
%       resid    2-norm of resulting residual (for final x)

x=x0; 
xold=x0+1; 
iter = 0;
while (norm(x-xold)>rtol*norm(x) & iter<1000)
  iter = iter + 1; xold = x;
  % solve C*xi=b (C=C(x), b=-d=-(f(x)-z)) in least square sense
  [rx,C]=feval(f,x);        % residual at x
  xi = C \ -rx;  % C*xi=-rx
  if withmini==1  % minimization [Schwarz sec. 7.4.2]
    tt = 1; Fx = rx'*rx;
    while (1)
      y = x + tt*xi; ry = feval(f,y); Fy = ry'*ry; 
      if (Fy < Fx+epsilon) x = y; break; % accept y as next x
      else tt = tt/2; % reduce tt and try again 
      end
    end
  else x = x + xi; % 'normal' newton (no minimization, i.e. tt=1)
  end
end
resid=norm(feval(f,x));
if (iter==1000 | resid>10) error('no convergence'); end
%-------------------------------------------------------------------------
function [F,J]=funnew(x)
z=[3.85 2.95 2.63 2.33 2.24 2.05 1.82 1.8 1.75]';
t=[0 0.5 1 1.5 2 3 5 8 10]';
F = [x(3)+x(4).*exp(x(1)*t)+x(5).*exp(x(2)*t)];
F=F-z;
if nargout > 1
o=ones(size(t));
J = [x(4)*t.*exp(x(1)*t) x(5)*t.*exp(x(2)*t) o exp(x(1)*t) exp(x(2)*t)];
end


%%%lÖSUNG 4 i:
%a)
% x0 =
%   -0.500000000000000
%   -2.000000000000000
%    1.750000000000000
%    1.200000000000000
%    0.800000000000000
% iter =
%     13


% x =
%   -0.555250292810919
%   -3.383579794262234
%    1.757739464441553
%    1.421016220019200
%    0.670663941706279
% resid =
%    0.077097085229299
% iter =
%     13


%b)
% x =
%   -0.555250292810919
%   -3.383579794262234
%    1.757739464441553
%    1.421016220019200
%    0.670663941706279
% resid =
%    0.077097085229299



%%%lÖSUNG 4 ii:


% ??? Error: File: Aufgabe_4.m Line: 23 Column: 1
% Unbalanced or unexpected parenthesis or bracket.
% 
% x0 =
%     -1
%     -1
%      1
%      1
%      1
% Warning: Rank deficient, rank = 3,  tol =   5.9952e-015.
% > In Aufgabe_4>gaussnewton at 51
%   In Aufgabe_4 at 19
% Warning: Rank deficient, rank = 0,  tol =   1.#INFe+000.
% > In Aufgabe_4>gaussnewton at 51
%   In Aufgabe_4 at 19
% ??? Error using ==> Aufgabe_4>gaussnewton at 64
% no convergence
% 
% Error in ==> Aufgabe_4 at 19
% [iter,x,resid] = gaussnewton(x0,rtol,epsilon,withmini,@funnew)