Diferencas divididas

Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original

// Programa para calcular as diferenças divididas (da ordem 0 até a ordem n) de uma função f tabelada
// Vou armazenar na matriz F os valores de x(i); os valores de f[x(i)]=y(i) e todas a outras diferenças divididas (para i=1,..,n+1). A primeira coluna armazeno x, na segunda as diferenças divididas de ordem zero, na coluna 3 as diferenças divididas de ordem 2,...
clf
clear
// Iniciando contagem de tempo
tic()
//Definindo n (numero de linhas da matriz F será n+1, eo de colunas será n+2)
n=10;
// Definindo os n+1 valores de x
F(1,1)= 0;
F(2,1)= 1;
F(3,1)= 5;
F(4,1)=4;
F(5,1)=-1;
F(6,1)=7;
F(7,1)= 3;
F(8,1)=2;
F(9,1)=0.5;
F(10,1)=6;
F(11,1)=1.3;
// Definindo do valores de y
F(1,2)=1;
F(2,2)=-1;
F(3,2)=3;
F(4,2)=1;
F(5,2)= 5;
F(6,2)= 6;
F(7,2)= 7;
F(8,2)=2;
F(9,2)=0;
F(10,2)=9;
F(11,2)=3;
// Construindo diferenças divididas
for j=3:n+2
 for i=1:n+3-j
 F(i,j)=(F(i+1, j-1)-F(i, j-1))/(F(i+j-2, 1)-F(i, 1));
 end
end
//Mostrando diferenças divididas
disp(F, 'Diferencas divididas');
x=poly(0,"x");
// Definindo vetor de produtórios de x-xi para i=0..n (ou 1..n+1 pois aqui não posso usar o indice 0)
for l=1:n+1 
 q(l)=x-F(l,1); 
end
A(1)=1;
for r=2:n+1
 A(r)=A(r-1)*q(r-1);
end
// Montando polinomio interpolador
soma=F(1,2);
for m=2:n+1
 soma=soma+F(1,m+1)*A(m);
end
pn=soma; 
//Mostrando polinomio interpolador
disp(pn, 'Polinomio interpolador pn(x)=');
// Mostrando Grafico
plot2d(F(:,1),F(:,2), -4);
M1=max(F(:,1)); 
M2=min(F(:,1));
h=abs(M1-M2)/50;
function y=polinomio(w), y=horner(pn,w), endfunction
X=[M2:h:M1]';
plot2d(X,polinomio(X), 4);
//Encerando contagem de tempo
toc()
t=toc()