clear
close all


n=20;  % n = Grad des Approximationspolynoms
m=100;  % Anzahl der Berechnungsdatenpunkte




%Auswertugnspunkte
t=linspace(-2,4,m)';





% Startwert des Newton-Verfahrens

   
x0=-0.3;    % --> Konvergenz gegen x = -0.599
%x0=0.4;    % --> Konvergenz gegen x = 0.3199
%x0=1;    % --> Konvergenz gegen x ~= 0.315740980638456
%x0=2    ;% --> Konvergenz gegen x ~= 2.257845892125693
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Funktion = @(x) atan(x) .* exp(-(x-1).^2) .* (x-1).^3;
  
% Alternative Testfunktionen  
% Funktion = @(x) (x-0.3).^3;
% Funktion = @(x) atan(x)*(x-0.5).^3;
% Funktion = @(x) (x-1).^2;











%% Approximation 


f=zeros(m,1);
f=Funktion(t);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Aufstellen der Größen für das lineare Programm %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
    alpha=zeros(n+1,1); % gesuchter Koeffizientenvektor
    c=zeros(n+1,1); c(n+1)=1;

    b=zeros(2*m,1);
    b(1:m)=-f;
    b(m+1:2*m)=f;

    Aapp=zeros(2*m,n+1);
    A=zeros(m,n)    
   
   for k=1:m
        for j=1:n
            A(k,j) = t(k)^(j);
        end
   end
    
    Aapp(1:m,1:n)=-A;
    Aapp(m+1:2*m,:)=-Aapp(1:m,:);
    Aapp(:,n+1)=-ones(2*m,1);

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


options = optimset('LargeScale','off','Simplex','on','Display','iter');  %nutze Simplex



% Bestimmung des Koeffizientenvektors
alpha = linprog(c,Aapp,b,[],[],[],[],[],options);
alpha=alpha(1:n);




%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Newton-Verfahren mit Interpolationspolynom und Regula Falsi %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





% Ableitungsmatrix   (es gibt keinen konstanten Term im
% Interpolationspolynom)
    B=zeros(m,n);
    for k=1:m
        for j=1:n
            B(k,j) = j*t(k)^(j-1) ;
        end
    end

    
    
% 2. Ableitungsmatrix
    C=zeros(m,n-1);
    for i=1:m
        for j=1:n-1
            C(i,j) = (j+1)*j*t(i)^(j-1);
        end
    end


    
 % Koeffizienten   der 2.Ableitung
    beta_2=zeros(n-1,1);
    for j=1:n-1
        beta_2(j) = alpha(j+1);
    end

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







M=50; % Anzahl der Iterationen

c=zeros(n-1,1); % Vektor, mit dem die 1. Ableitung des Polynoms in x(k) berechnet wird
d=zeros(n-2,1); % Vektor, mit dem die 2. Ableitung des Polynoms in x(k) berechnet wird

x=zeros(M,1);
x(1)=x0;

for m=2:M
    
    for j=1:n
        c(j) = j*x(m-1)^(j-1);       % p'(x(m)) = c'*beta
    end
    
    for j=1:n-1
        d(j) = (j+1)*j*x(m-1)^(j-1); % p''(x(m)) = d'*beta
    end
    
    x(m) = x(m-1) - (c'*alpha) / (d'*beta_2);
    
end








%% Plots


figure;
plot(t,f,'*r'); 
hold on
plot(t,A*alpha); title('Ausgangspunkte und Interpolierendes Polynom');
%saveas(gcf,strcat('figure1polynom.fig'))


figure; 
plot(t,B*alpha); title('Ableitung des interpolierenden Polynoms'); hold on; plot(t,0);
%saveas(gcf,strcat('figure2ersteableitung.fig'))


figure;
plot(t,C*beta_2); title('Zweite Ableitung des interpolierenden Polynoms'); hold on ; 
%saveas(gcf,strcat('figure3zweiteableitung.fig'))


figure;
plot(x,'o-'); title('Berechneter Iterationspunkt pro Iterationsschritt');
%saveas(gcf,strcat('figure4konvergenzbeiminus03.fig'))




