% aufgabe 9, numerik 1, Bibiana 5600816, Marina 5622011, Thomas 5610641
%% eingabe: n x-koordinaten als vektor
%%          n y-koordinaten als vektor
%% ausgabe: koeffizentenvektor den man so in polyval einsetzen kann.
function p=interpol(xi,yi);
c=1:length(yi);
%% zunaechst einmal berechnen wir die newtonkoeffizienten des polynomes.
%% anders als in der vorlesung vorgelesen wird hier aber nicht die untere matrix
%% berechnet, sondern die obere
%%
%%   x - x - x - x
%%   x / x / x /
%%   x / x / 
%%   x /
%% die newtonkoeffizienten stehen dann oben.
%%
%% jeder neue eintrag wird aus den dividierten differenzen gebildet.
%% ueber dem bruchstrich bildet man die differenz der beiden funktionswerte die
%% 1. in der gleichen zeile stehen  g(j)
%% 2. in der zeile darunter stehen  g(j+1)
%% und teilt das durch 
%% 1. die x-koordinate der zeile xi(j)
%% 2. minus die x-koordinate, i zeilen darunter xi(j+i)
%% dadurch entsteht eine matrix, weil es aber ausreicht die spalte vorher (g)
%% zubetrachten, konzentriert sich diese implementierung nur auf den aktuellen
%% spaltenvektor (h).
g=yi; %% in der ersten spalte stehen die uebergebenen funktionswerte
c(1)=g(1);
for i=1:length(yi)-1
	h=1:length(g)-1;
	for j=1:length(g)-1
		h(j)=(-g(j+1)+g(j))/(-xi(j+i)+xi(j));
	end
	c(i+1)=h(1);  %% der neue newtonkoeffizient steht ganz am anfang.
	g=h;
end
c
%% um aus der newtondarstellung eine vernuenftige koeffizientendarstellung
%% zu gewinnen muss man sich im prinzip so eine matrix bauen:
%% 1        0          0    0       
%% 1       -a          0    0      a
%% 1     -a-b         ab    0      b
%% 1   -a-b-c ab-(-a-b)c -abc      c
%% also links alles einsen, und dann die zeile darueber plus einen hoch und links mal das neue x.
%% dann multipliziert man c0 mit der hauptdiagonalen, addiert das ganze, und
%% hat a0 der ueblichen polynomdarstellung. 
%% dann c1 mit der diagonalen darunter, addieren, und schon hat man den 
%% koeffizeiten fuer x^1.
%% dann c2 mit der diagonalen DARUNTER, addieren, koeffizeiten fuer x^2.
%% usw. war ganz schoen heftig, sich das auszudenken. aber klappt ;-)
a=zeros(length(c),length(c));
for i=1:length(a)
	a(i,1)=1;	
end
for i=2:length(xi)
	for j=1:length(xi)-1
		a(i,j+1)=a(i-1,j+1)-xi(i-1)*a(i-1,j);
	end
end
%% und jetzt den richtigen koeffizientenvektor berechnen
p=zeros(1,length(a));
for i=1:length(a)
	for j=i:-1:1
		p(length(a)-i+j)=p(length(a)-i+j)+a(i,j)*c(i);
	end 
end