

#algorytm Newtona vs algorytm roznic dzielonych (z def rekurencyjnej RD)
# znajdowania roznic dzielonych czyli wielomianu interpolacyjnego Lagrange'a w 
#bazie Newtona

0==0;

function pv=horneru(rd,x,q,k);
#uogolniony alg Hornera dla 
#bazy Newtona dla wezlow z x 
#INPUT
#rd - wspolczynniki w bazie dla wezlow x(1),...,x(k)
#x - wezly x(1),...,x(k)
#q - wektor liczb w ktorych chemy polcizyc wartosc iwielomianu
#k - stopien wielomianu
#OUTPUT
#y wektor wartosci wiel w q (tego samego wymiaru)
 pv=rd(k+1);
 for j=k:-1:1,
  pv=pv.*(q-x(j))+rd(j);
 endfor
end

function rd=newton(x,y)
#znajduje roznice dzielone dla wezlow z x i wartosci w y
#metoda newtona
 n=length(x);
 nbp=rd=zeros(n,1)';
 rd(1)=y(1);
 for k=2:n,
  nbp(k)=1;
  rd(k)=(y(k)-horneru(rd,x(1:k-2),x(k),k-2))/horneru(nbp,x(1:k-1),x(k),k-1);
  nbp(k)=0;
 endfor

endfunction


function rd=algrozdziel(x,y)
#znajduje roznice dzielone dla wezlow z x i wartosci w y
#alg roznic dzielonych (z definicji rekurencyjnej)
 n=length(x);
 rd=y;
 for k=2:n,
  rd(k:n)=(rd(k:n)-rd(k-1:n-1))./(x(k:n)-x(1:n-k+1));
 end
endfunction


#testujemy
x=-1:4;
y=x.^5+x.*x-10;
n=length(x)-1;#stopien wielomianu
rd=newton(x,y) #czy rd(6)=1?
rd1=algrozdziel(x,y) #czy to samo?
rd-rd1
s=linspace(min(x)-0.5,max(x)+0.5,100);
plot(x,y,"+",s,horneru(rd,x,s,n),";czy + na wykresie?;");pause(1);

#przetestujmy dla sin
y=sin(2*x);
rd=newton(x,y) #czy rd(6)=1?
rd1=algrozdziel(x,y) #czy to samo?
rd-rd1
s=linspace(min(x)-0.5,max(x)+0.5,100);
plot(x,y,"+",s,horneru(rd,x,s,n),s,sin(2*s),";sin(2x);");

#przetestujmy dla duzej ilosci wezlow i wysokich stopni
#interpolacje na wezlach rownoodleglych na [0,1] dla losowych wartosci y
N=2;M=45;
erin=zeros(M,3);
errd=zeros(M,1);
for k=1:M,
 x=linspace(0,1,N);
 y=randn(size(x));
 rd1=newton(x,y); #czy rd(6)=1?
 rd2=algrozdziel(x,y); 
 p=polyfit(x,y,N-1);
 errd(k)=norm(rd1-rd2,'inf');
 erin(k,1)=norm(y-horneru(rd1,x,x,N-1),'inf');
 erin(k,2)=norm(y-horneru(rd2,x,x,N-1),'inf');
 erin(k,3)=norm(y-polyval(p,x),'inf');
 printf("St wiel =%2d norma sup w wezlach f-L_n f=%f\n",N-1,erin(k,2));
 N++;
endfor
semilogy(1:M,erin,1:M,ones(M,1))
#0ba alg dla roznic sa wyraznie lepsze niz polyfit ze wzgledu na zaburzenia numeryczne





