%Aufgabe 32
%Natürliche Spline Interpolation mit äquidistanten Stützstellen
clear all
n = 41;      % Anzahl der Knoten
a0 =-5;      % Intervallanfang
b0 = 5;      % Intervallende
u = linspace(a0,b0,n);  % Knoten 
f = 1./(1+u.*u);        % Zu approximierende Funktion
h = (b0-a0)/(n-1);
%Momenten Gleichung
q = linspace(0.5,0.5,n-3);
A = 2*eye(n-2)+diag(q,-1)+diag(q,1);
r = (6/(2*h));
for i = 1:n-2
    D(i) = r*(f(i+2)-2*f(i+1)+f(i))/h;
end
M = A\D';
M = [0;M;0]';
% Koeffizienten des Splines
a = f(1:n-1)
b = (f(2:n)-f(1:n-1))/h-(h/6)*(2*M(1:n-1)+M(2:n))
c = 0.5*M(1:n-1)
d = (1/(6*h))*(M(2:n)-M(1:n-1))
%Auswertungen
for i = 1:n-1
    x(i,:) = a0+(i-1)*h + linspace(0,h,n);
end
for j = 1:n-1
    F(j,:) = 1./(1+x(j,:).*x(j,:));
    P(j,:) = [d(j) c(j) b(j) a(j)];
    S(j,:) = polyval(P(j,:),x(j,:)-x(j));
end    
for k = 1:n-1
    subplot(2,2,1)
    plot(x(k,:),F(k,:),'r-',x(k,:),S(k,:),'g-',x(k,:),F(k,:)-S(k,:),'b-');
    hold on
end
for k = 1:n-1
    subplot(2,2,2)
    for j = 1:n-1
        G(j,:)= -2*x(j,:).*(F(j,:).^2);
        Q(j,:) = [3*d(j) 2*c(j) b(j)];
        R(j,:) = polyval(Q(j,:),x(j,1:10)-x(j));
    end
        plot(x(k,:),G(k,:),'r-',x(k,1:10), R(k,:),'g-')
    hold on
end
    for k =1:n-1
        subplot(2,2,3)
        for j = 1:n-1
           H(j,:) = (-2+6*x(j,:).*x(j,:))./((1+x(j,:).^2).^3);
           DDS(j,:) = M(j)*(x(j,n)-x(j,:))/h + M(j+1)*(x(j,:)-x(j,1))/h;
        end
        plot(x(k,:),H(k,:),'r',x(k,:),DDS(k,:),'g')
        hold on
    end
    for k =1:n-1
        subplot(2,2,4)
        for j = 1:n-1
            K(j,:) = (24*x(j,:)-24*x(j,:).^3)./((1+x(j,:).^2).^4);
            DDDS(j,1:n) = -M(j)+M(j+1);
        end
        plot(x(k,:),K(k,:),'r', x(k,:),DDDS(k,:),'g')
        hold on
    end