Calcula polinomio de lagrange para interpolar data.

combinatoria.m

function res = combinatoria(n,r)
  % Esta es una combinatoria especial, dado que n es un
  % polinomio y r es un escalar (para el caso de newton
  % adelantado y atrasado)

  res = n; % Aca almacenaremos el resultado.
  
  for i=1:r-1
    res = res * (n - i);
  end
  
  % Finalmente dividimos por el factorial de r
  res = res / factorial(r);
end

interpolacion.m

function [sol] = interpolacion(f, pts, met)
  if strcmp(met, "lag")
    sol = interpolacion_lagrange(f, pts);
  end
  
  if strcmp(met, "reg")
    sol = interpolacion_datrasada(f, pts);
  end
  
  if strcmp(met, "prog")
    sol = interpolacion_dadelantada(f, pts);
  end
end

interpolacion_dadelantada.m

function [sol] = interpolacion_dadelantada(f, pts)
  symbols; % Inicializamos uso de simbolos.
  x = sym('x');

  [filas, columnas] = size(pts);
  n = columnas;
  puntos = zeros(2, n); % Matriz de 2 filas y n columnas, que contendra los x en una fila
                        % y los f(x) en la otra.
  
  % Generamos la matriz de puntos
  for i=1:n
    puntos(1,i) = pts(i);
    puntos(2,i) = f(pts(i));
  end
  
  matdelta = zeros(n,n); % Esta matriz contendra la matriz de diferencias
  
  % Generamos la primera columna, que son los f(x_i)
  for i=1:n
    matdelta(i,1) = puntos(2,i);
  end
  
  % Generamos el resto de las diferencias, como una matriz triangular inf.
  for i=2:n
    for j=2:i
      matdelta(i,j) = (matdelta(i,j-1) - matdelta(i-1,j-1)) / (puntos(1,i) - puntos(1,i-j+1));
    end
  end
  
  % Info necesaria para la sumatoria que viene
  h = abs(puntos(1,2) - puntos(1,1));
  %s = (x - puntos(1,n))/h;
  s = (x - puntos(1,1))/h;
  
  % Sumatoria, casi listos.
  suma = puntos(2,1); % Primer elemento de la sumatoria
  for i=1:(n-1)
    suma = suma + combinatoria(s, i) * factorial(i) * h^i * matdelta(i+1,i+1);
  end
  
  sol = obtener_coeficientes(suma, x, n);
end

interpolacion_datrasada.m

function [sol] = interpolacion_datrasada(f, pts)
  symbols; % Inicializamos uso de simbolos.
  x = sym('x');

  [filas, columnas] = size(pts);
  n = columnas;
  puntos = zeros(2, n); % Matriz de 2 filas y n columnas, que contendra los x en una fila
                        % y los f(x) en la otra.
  
  % Generamos la matriz de puntos
  for i=1:n
    puntos(1,i) = pts(i);
    puntos(2,i) = f(pts(i));
  end
  
  matdelta = zeros(n,n); % Esta matriz contendra la matriz de diferencias
  
  % Generamos la primera columna, que son los f(x_i)
  for i=1:n
    matdelta(i,1) = puntos(2,i);
  end
  
  % Generamos el resto de las diferencias, como una matriz triangular inf.
  for i=2:n
    for j=2:i
      matdelta(i,j) = (matdelta(i,j-1) - matdelta(i-1,j-1)) / (puntos(1,i) - puntos(1,i-j+1));
    end
  end
  
  % Info necesaria para la sumatoria que viene
  h = abs(puntos(1,2) - puntos(1,1));
  %s = (x - puntos(1,n))/h;
  s = (puntos(1,n) - x)/h; % Lo pongo asi porque al hacer "-s" abajo, tira error.
  
  % Sumatoria, casi listos.
  suma = puntos(2,n); % Primer elemento de la sumatoria
  for i=1:(n-1)
    %suma = suma + (-1)^i * combinatoria(-s, i) * factorial(i) * h^i * matdelta(n,i+1);
    suma = suma + (-1)^i * combinatoria(s, i) * factorial(i) * h^i * matdelta(n,i+1);
  end
  
  sol = obtener_coeficientes(suma, x, n);
end

interpolacion_lagrange.m

function [sol] = interpolacion_lagrange(f, pts)  
  symbols; % Con esto inicializamos el uso de simbolos (como 'x')

  [filas, columnas] = size(pts);
  n = columnas;
  puntos = zeros(2, n); % Matriz de 2 filas y n columnas, que contendra los x en una fila
                        % y los f(x) en la otra.
  
  % Generamos la matriz de puntos
  for i=1:n
    puntos(1,i) = pts(i);
    puntos(2,i) = f(pts(i));
  end
  
  polinomio = 0;
  
  x = sym('x'); % Necesitamos este simbolo para calcular el polinomio.

  for i=1:n
    l = 1; % Variable auxiliar para calcular el polinomio
    for j=1:n
      if i ~= j
        numerador = x - puntos(1,j);
				denominador = puntos(1,i) - puntos(1,j);
				l = l*(numerador/denominador);
      end
    end
    
    polinomio = polinomio+l*puntos(2,i);
  end
  
  polinomio = expand(polinomio);
  
  sol = obtener_coeficientes(polinomio, x, n); % Esto nos devuelve la matriz con los
                                               % coeficientes, a partir del polinomio.
end

obtener_coeficientes.m

function coefs = obtener_coeficientes(polinomio, simbolo, n)
  % La idea es simple, derivar y evaluar en cero n veces.
  % Lo complejo es que derivar un polinomio "destruye" los 
  % coeficientes.
  %
  % By Ian Murray
  
  coefs = zeros(1,n); % Aca guardaremos los coeficientes.
  
  for i=1:n
    % Evaluamos en cero, y el elemento tenemos que dividirlo
    % por el factorial de (i-1), para corregir el error que
    % genera derivar el polinomio
    % El indice n-i+1 es para que el vector de coeficientes
    % quede en el orden correcto.
    coefs(n-i+1) = to_double(subs(polinomio, simbolo, 0) / factorial(i-1));
    
    % Derivamos
    polinomio = differentiate(polinomio, simbolo, 1);
  end
  
  % Listo :D muy simple
end