clear
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N=10;           % Anzahl der Messpunkte

a=0; b=10;      % Messintervall



lambda_1=0.5;        % wahre Parameter
lambda_2=0.2; 


M=20;           % Anzahl der Iterationen beim Newton-Verfahren


lambda_Newton=zeros(2,M+1); % Gauß-Newton-Iterierte 
lambda_Newton(1,1)=1; % Startwerte für die Iteration
lambda_Newton(2,1)=1;


warning('')   %Warnung auf "leer" setzen



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Wahl der Ansatzfunktion
ansatzfunk='linear'
ansatzfunk='exponential'

t=linspace(a+(b-a)/N,b,N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
% Datenerzeugung
% 
% 
% f=zeros(1,N);
% 
% 
% switch ansatzfunk
%         case 'linear'
%              f = lambda_1 + lambda_2*t;
%         case 'exponential'
%              f = lambda_1*exp(lambda_2*t);
%              
%              f = lambda_1 + lambda_2*t;
%              
% end    
% 
% 
% 
% % % Gleichverteilte Zufallsvariable auf (0,1)
% GV=rand(1,N);
% delta=0.2;
% % Gleichverteilte Zufallsvariable auf [-delta,delta]
% stoerung=zeros(1,N); 
% stoerung=delta*[2*(GV-0.5)];   
% 
% f_delta = zeros(1,N);
% f_delta = f + stoerung;



% Direkte Vorgabe der Werte für die Exponentialfunktion 
% f_delta = [0.74  0.64  1.03  1.01  1.53  1.60  1.90 2.37  3.07  3.68]; %Exponentialfunktion
   f_delta = [0.56 1.01 1.02 1.31 1.36 1.72 1.80 2.16 2.37 2.59]    %  Lineare Funktion

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%% 
% Gauß-Newton-Verfahren %
%%%%%%%%%%%%%%%%%%%%%%%%%


d=zeros(1,2);   % Richtungsvektor für Gauß-Newton-Iteration (Schrittweite = 1)
F=zeros(1,N);   % Fehlerquadrate der Komponenten
J_F=zeros(2,N); % Jacobi-Matrix der Zielfunktion im jeweiligen Iterationspunkt

 




if strcmp(ansatzfunk,'linear')
 % Lineare Ansatzfunktion --> Jacobi-Matrix ist unabhängig von lambda
            J_F(1,:) = 1;
            J_F(2,:) = t;
end





for m=1:M       
   
    
        
    switch ansatzfunk
        case 'linear'
            F = lambda_Newton(1,m) +  lambda_Newton(2,m)*t - f_delta;
        
         case 'exponential'
            % Exponentielle Ansatzfunktion
             F =  lambda_Newton(1,m)*exp( lambda_Newton(2,m)*t) - f_delta;
  
             J_F(1,:) = exp(lambda_Newton(2,m)*t);
             J_F(2,:) =  lambda_Newton(1,m)*t.*exp( lambda_Newton(2,m)*t);         
            
    end
      
      % Abbruch,  falls die Kondition der MAtrix (J_F*J_F') zu schlecht ist
    if  strcmp(lastwarn,'')
          lambda_Newton(:,m+1) =lambda_Newton(:,m)- (J_F*J_F')\(J_F*F');
    else   
          lambda_Newton(:,m+1:M+1)=repmat(lambda_Newton(:,m),1,M+1-m)
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for m=1:M+1
f_Newton=zeros(1,N);
switch ansatzfunk
    case 'linear'
        f_Newton = lambda_Newton(1,m)+ lambda_Newton(2,m)*t;
    case 'exponential'  
        f_Newton = lambda_Newton(1,m)*exp(lambda_Newton(2,m)*t); 
end



% Fehlerquadratsumme 
Fehlerquadratsumme(m) = 0.5*norm(f_Newton-f_delta,2)^2;
Fehlersumme(m)=sum(abs(f_Newton-f_delta));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
figure
plot(t,f_delta,'ro');  hold on;
plot(t,f_Newton,'g-'); 
%plot(t,f); % nur falls wahre Fuktion bekannt
title('Daten (rot), Funktion für die berechneten Parameter (grün)'); 
saveas(gcf,strcat('figure1gausnewton.fig'))  

figure;
plot(0:M, lambda_Newton(1,:),'b-o'); hold on;
plot(0:M, lambda_Newton(2,:),'r-o'); % title('Komponenten der Iterierten');
title('Werte von \lambda_1 und \lambda_2'); 
saveas(gcf,strcat('figure2gausnewton.fig'))  


if  ~strcmp(lastwarn,'')
    display('Achtung die Matrix war fast singulär, das Ergebnis ist vermutlich unbrauchbar. Wählen Sie einen anderen Startwert! \n')
end


display(['Die Summe des punktweisen Fehlers beträgt nach ', num2str(M),' Iterationen: ',num2str(Fehlersumme(M+1))])
display(['Die Fehlerquadratsumme  beträgt nach ', num2str(M),' Iterationen: ',num2str(Fehlerquadratsumme(M+1))])

