
%% Numerik D-MAVT, RW (Nipp), SS06
%
% This is the M-file solution to the sixth exercise sheet.
% To run, save this file to your account and ensure that the "Current directory"
% in MATLAB is the same as that where the file is saved.

% Then to run the file type "s6" in the command window

%Aufgabe 1
function s6

display(['Aufgabe 1']);
format long;
a=[0.25,0.5,0.75];
tl=[0:10]';
il=[0 8 12 13 13 12 10 8 7 6 5]';

A = [exp(-tl*0.25) exp(-tl*0.5) exp(-tl*0.75) ] ;
xne = (A'*A)\(A'*il) % Normal Gleichung
[U,S,V]=svd(A); c=U'*il;t=S(1:3,1:3)\c(1:3); % SVD decomposition
xsvd=V*t

%%%%%%% Antwort %%%%%%%%%%%%%%%%
%xne =
%
%  63.37110732929723
% -85.29263723977263
%  21.92968792921701
%
%
%xsvd =
%
%  63.37110732929711
% -85.29263723977215
%  21.92968792921661



%Aufgabe 3 Gauss Newton with and without minimization
display(['Aufgabe 3']);
clear all; format compact; format long;
rtol=1e-8; epsilon=1e-12;
% (i) 'good' initial values
x0=[-.5 -2 1.75 1.2 .8]';
% no minimization
[iter,x,resid] = gaussnewton(x0,rtol,epsilon,0,@funnew); 
% iter =
%     13
% x =
%   -0.55525029281092
%   -3.38357979426224
%    1.75773946444155
%    1.42101622001920
%    0.67066394170628
% resid =
%    0.07709708522930
% with minimization
[iter,x,resid] = gaussnewton(x0,rtol,epsilon,0,@funnew) 
% iter =
%     13
% x =
%   -0.55525029281092
%   -3.38357979426224
%    1.75773946444155
%    1.42101622001920
%    0.67066394170628
% resid =
%    0.07709708522930
% (ii) 'bad' initial values
x0=[-1 -1 1 1 1]';
% without minimization (comment out to give result)
%[iter,x,resid] = gaussnewton(x0,rtol,epsilon,0,@funnew); 
% ??? Error using ==> gaussnewton
% no convergence
% with minimization
[iter,x,resid] = gaussnewton(x0,rtol,epsilon,1,@funnew),
% iter =
%     22
% x =
%   -0.55525029251988
%   -3.38357978824621
%    1.75773946439759
%    1.42101621929148
%    0.67066394246689
% resid =
%    0.07709708522930
%===============================================================================
%Aufgabe 3 using fsolve
format long
%good values
x0=[-0.5 -2 1.75 1.2 0.8]';
tt=1;
options=optimset('Display','iter','Jacobian','on','NonlEqnAlgorithm','gn',...
'TolX',1e-8,'TolFun',1e-4);
[x,F,exitflag,output,JAC]=fsolve(@funnew,x0,options);
display(['Computed Solution using fsolve '])
x
display(['Norm of residual using fsolve = ',num2str(norm(feval(@funnew,x)))])
%bad values
x0=[-1 -1 1 1 1]';
tt=1;
options=optimset('Display','iter','Jacobian','on','NonlEqnAlgorithm','gn',...
'TolX',1e-8,'TolFun',1e-4);
[x,F,exitflag,output,JAC]=fsolve(@funnew,x0,options);
display(['Computed Solution using fsolve '])
x
display(['Norm of residual using fsolve = ',num2str(norm(feval(@funnew,x)))])

%Iteration  Func-count    Residual     Step-size      derivative 
%     0           1       0.0172006
%     1           5      0.00919113        0.607       -4.86e-05
%     2          10      0.00600827         1.07        -8.7e-05
%     3          14      0.00594404         1.01       -7.13e-09
%     4          15      0.00594396         1.01       -1.16e-07
%Optimization terminated: directional derivative along
% search direction less than TolFun and infinity-norm of
% gradient less than 10*(TolFun+TolX).
%Computed Solution using fsolve 
%
%x =
%  -0.55458081047393
%  -3.36942834566893
%   1.75764315258730
%   1.41929005463807
%   0.67243228158153
%Norm of residual using fsolve = 0.077098
%                                                     Directional 
% Iteration  Func-count    Residual     Step-size      derivative 
%     0           1         6.54092
%     1           6        0.123048         1.39        -0.00149
%Conditioning of Gradient Poor - Switching To LM method
%     2          36        0.111271     1.49e-16        2.79e+13              0 
%     3          40        0.109105       0.0205        -0.00118              0 
%     4          45       0.0531572         2.37        -2.2e-06         4.9827 
%     5          50       0.0360449          1.5        2.58e-07        1.48046 
%     6          55       0.0326723         1.89       -8.58e-07       0.740232 
%     7          60       0.0252064         3.96       -7.28e-06       0.255971 
%     8          65       0.0105199         1.65       -3.54e-05       0.051538 
%     9          69      0.00789559         1.18       -6.02e-07      0.0194473 
%    10          74      0.00665012         3.09       -3.94e-05      0.0089406 
%    11          79      0.00610391         1.82       -2.53e-05     0.00209378 
%    12          84      0.00598031         1.31       -3.09e-07    0.000707674 
%    13          89      0.00594667         1.53       -4.41e-07    0.000306078 
%    14          93      0.00594412         1.17       -3.44e-10     0.00011988 
%    15          94      0.00594397            1       -2.47e-07    5.53187e-05 
%Optimization terminated: directional derivative along
% search direction less than TolFun and infinity-norm of
% gradient less than 10*(TolFun+TolX).
%Computed Solution using fsolve 
%x =
%  -3.36580274428498
%  -0.55465242396258
%   1.75768435060471
%   0.67237346514705
%   1.41919222123728
%Norm of residual using fsolve = 0.077098
%For the Bad intial values Fsolve converges to a different solution compared
%to the Hand-Made Gauss Newton Function. Both satisify the same residual
%requirements. By choosing different intial values for the Hand Made Gauss
%Newton solver a similar solution can also be obtained.eg
x0=[-3 -0.5 1 1 1]';
[iter,x,resid] = gaussnewton(x0,rtol,epsilon,1,@funnew),
%iter =
%    13
%x =
%  -3.38357979234395
%  -0.55525029271812
%   1.75773946442754
%   0.67066394194881
%   1.42101621978716
%resid =
%   0.07709708522930


function [iter,x,resid] = gaussnewton(x0,rtol,epsilon,withmini,f)
%
% Gauss-Newton for nonlinear overdetermined system, 
% optimized for Numerik D-MAVT SS 04, K.Nipp, Serie 6, A4
%
% 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)
% auth: ledger@math.ethz.ch

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
