Newton interpolation method algorithm in Matlab

Newton interpolation method algorithm in Matlab :

function [yy,E,d,n_iter]=interpola_newton(x,y,xx,error) % Syntaxis: [yy,d,n_iter]=interpola_newton(x,y,xx,error) % ---------------------------------------------------------------------- % Devuelve el valor interpolado yy en un determinado punto xx solicitado % mediante el metodo de las diferencia divididas de Newton, empleando (n-1) % nodos para la interpolación y un nodo para determinar el error de la % interpolacion. % % Entrada: % x --> Vector que contiene las abcisas de los nodos de interpolacion % y --> Vector que contiene las ordenadas de los nodos de interpolacion % xx --> Abcisa del nodo que se desea incorporar % error --> Tolerancia de la interpolacion {default | 0.05} % % % Salida: % yy --> Valor interpolado para el nodo con abcisa xx % E --> Error real de la interpolacion % d --> Tabla de Diferencias Divididas % n_iter --> Numero de iteraciones necesitadas para verificar la condicion % de parada % % % Ejemplo: % % x=[1.00, 1.05, 1.10, 1.15, 1.20, 1.25, 1.30]; % y=[1.00000, 0.97350, 0.95135, 0.93304, 0.91817, 0.90640, 0.89747]; % xx=1.17; error=0.00005; % [yy,E,d,n_iter] = interpola_newton(x,y,xx,error) % % yy = 0.9267 % % E = 2.0944e-05 % % d = % % -0.5300 0.8700 -0.6800 0.7333 -1.3333 6.2222 % -0.4430 0.7680 -0.5333 0.4000 0.5333 0 % -0.3662 0.6880 -0.4533 0.5333 0 0 % -0.2974 0.6200 -0.3467 0 0 0 % -0.2354 0.5680 0 0 0 0 % -0.1786 0 0 0 0 0 % % n_iter = 4 % % % Ver tambien: interpola_lagrange % % ---------------------------------------------------------------------- % % Elaborado por: % % Msc. Alexeis Companioni Guerra % Lic. Vianka Orovio Cobo % % Revision: 1.0 % Fecha: 30/09/2007 %% Validacion de las entradas if nargin<3 error('Verifique los datos de entrada.'); elseif nargin<4 error = 0.05; end if min(size(x))>1 || min(size(y))>1 error('Tanto x como y deben ser vectores.'); elseif max(size(x)) ~= max(size(y)) error('Verifique la dimension de los vectores.'); else x = x(:); y = y(:); n = length(x); end if sum(size(xx))~=2 error('La abcisa del nodo de entrada debe ser un escalar.'); end if sum(size(error))~=2 error('El error maximo de la interpolacion debe ser un escalar.'); end %% Conformacion de las Diferencias Divididas d = diff(y)./diff(x); for j=2:n-1 for k=1:n-j d(k,j) = ( d(k+1,j-1) - d(k,j-1) ) / ( x(k+j) - x(k) ); end end %% Seleccion de la primera fila de la matriz de Diferencias Divididas p = [y(1),d(1,:)]; %% Conformacion iterativa de la solucion yy = p(1); producto = xx-x(1); E = p(2)*producto; i = 2; while abs(E)>error && i<n yy = yy + E; producto = producto*(xx-x(i)); i = i + 1; E = p(i)*producto; end n_iter = i - 1;