INTERPOLACIÓN POLINOMIAL Y AJUSTE DE CURVAS

Hay dos temas que se tratarán aquí, es decir, interpolación y ajuste de curvas. La interpolación
es conectar los puntos de datos discretos de manera plausible para que uno pueda obtener
estimaciones razonables de puntos de datos entre los puntos dados (nodos).
La curva de interpolación pasa a través de todos los puntos de datos. El ajuste de curvas, por
otra parte, es encontrar una curva que mejor se podría indicar la tendencia de un conjunto dado
de datos. La curva no tiene que pasar por los puntos de datos. En algunos casos, los datos
pueden tener diferente exactitud / fiabilidad / incertidumbre y necesitamos la curva de ajuste
ponderado para procesar dichos datos por mínimos cuadrados.
Interpolación es una herramienta importante cuando no podemos evaluar rápidamente la función
en puntos intermedios. La interpolación consiste en estimar valores intermedios entre cada par
de puntos o nodos que componen una serie de datos.
En MATLAB podemos interpolar en más de una dimensión. El comando interp1 interpola datos
e una dimensión, el comando interp2 lo hace en dos dimensiones (es decir, si z=f(x,y)),
podemos interpolar entre valores x y y a fin de encontrar los valores intermedios de z); el
comando interp3 interpola datos en tres dimensiones, y el comando interpn interpola datos en
mayor dimensión. Con el comando help investiguemos cómo usar los comandos de
interpolación anteriores.
Consideremos, como ejemplo práctico, el siguiente conjunto de mediciones de temperaturas
tomadas de un equipo agroindustrial durante un tiempo de funcionamiento:
Tiempo, s Temperatura, °C
0 -17.77
1 -6.66
2 15.55
3 20
4 25
5 43.33

En MATLAB:
>> x=0:5;
>> y=[-17.77 -6.66 15.55 20 25 43.33]; % Podemos usar el comando interp1 para interpolar
(interpolación lineal)
% una temperatura que corresponda a cualquier instante entre 0.0 y 5.0 segundos. Si escribimos:
>> y1=interp1(x,y,1.5) % y
>> y2=interp1(x,y,4.5)
obtendremos 4.4450 y 34.1650, respectivamente.
Si el segundo argumento de interp1 es un matriz, el comando devuelve un vector fila con el
mismo número de columnas (cada valor devuelto interpola la columna de datos
correspondiente). Por ejemplo, supongamos que medimos las temperaturas en tres puntos
distintos de un equipo agroindustrial, obteniendo el conjunto de datos siguiente:
Tiempo, s Temp1, °C Temp2, °C Temp3, °C
0 -17.77 -17.77 -17.77
1 -6.66 -3.88 11.11
2 15.55 16.66 32.22
3 20 19.44 32.77
4 25 27.77 33.88
5 43.33 39.44 35.55

Así, podemos escribir:

x(:,1)=(0:5)';
y(:,1)=[-17.77 -6.66 15.55 20 25 43.33]';
y(:,2)=[-17.77 -3.88 16.66 19.44 27.77 39.44]';
y(:,3)=[-17.77 11.11 32.22 32.77 33.88 35.55]';
Si deseamos determinar valores de temperaturas interpolados en los tres puntos del equipo a los
1.5 segundos, escribimos:
>> temps=interp1(x,y,1.5)
temps =
4.4450 6.3900 21.6650

INTERPOLACIÓN POLINÓMICA DE LAGRANGE

function [l,L] = lagranp(x,y)

% Entrada : x = [x0 x1 ... xN], y = [y0 y1 ... yN]
%Salida: l = coeficientes de los polinomios de Lagrange de grado N
% L = coeficiente polinómico de Lagrange
N = length(x)-1; % el grado de polinomio
l = 0;
for m = 1:N + 1
P = 1;
for k = 1:N + 1
if k ~= m, P = conv(P,[1 -x(k)])/(x(m)-x(k)); end
end
L(m,:) = P; % coeficiente polinómico de Lagrange
l = l + y(m)*P; % polinomio de Lagrange
end

% En Command Window: -----------------------------------------

%do_lagranp.m
x = [-2 -1 1 2]; y = [-6 0 0 6]; % puntos de datos dados
l = lagranp(x,y) % encontrar el polinomio Lagrange
xx = [-2: 0.02 : 2]; yy = polyval(l,xx); %interpolar para [-2,2]
clf, plot(xx,yy,'b', x,y,'*') %trazar la gráfica
l=
1 0 -1 0

Se ejecuta la rutina "lagranp ()" para encontrar el l3 polinomio de tercer grado (x), que coincide
con los cuatro puntos dados {(-2, -6), (-1, 0 ), (1, 0), (2, 6)} y para comprobar si la gráfica de l3
(x) realmente pasa los cuatro puntos. Los resultados de la ejecución de este programa se
representan en la Fig. anterior
% significa l3(x) = 1 · x3 + 0 · x2 − 1 · x + 0

Otro programa: polLAGRANGE.m

clear;clc;
fprintf('\n*** POLINOMIO DE LAGRANGE ***\n\n')
%x=[5 -7 -6 0];
%y=[1 -23 -54 -954];
x=input('Vector X= ');
y=input('Vector Y= ');
p=0;
for k=1:length(x)
den=1;
L=[];% término de lagrange
for r=1:length(x)
if r~=k
xp=x; % hallando el numerador de L
xp(k)=[]; % elimino el término k-esimo
num=poly(xp);% polinomio de x
d1=(x(k)-x(r)); % denominador de L
den=den*d1;
end

end
L=num/den;
p=p+y(k)*L; % término de lagrange
end
Pl=poly2sym(p);% cambio de vector a expresión.
fprintf('\nPl(x)=')
pretty(Pl)% mostrando ordenadamente
fprintf('\n\n')
% graficando
Z=input('¿Deseas Graficar [Y]/[N]? : ','s');
switch Z
case 'Y';
hold on
xo=min(x):0.1:max(x);
yo=polyval(p,xo);
plot(x,y,'ro',xo,yo,'-r')
grid on;
disp('Fin del Programa')
case 'N';
disp('Fin del Programa')
end

POLINOMIO DE INTERPOLACIÓN DE NEWTON

Tabla. Tabla de diferencia dividida.

function [n,DD] = newtonp(x,y) % Interpolación
%Input:x= [x0 x1 ... xN]
% y = [y0 y1 ... yN]
%Output: n = coeficientes del polinomio de Newton grado N
N = length(x)-1;
DD = zeros(N + 1,N + 1); DD(1:N + 1,1) = y';
for k= 2:N+1
for m= 1:N+2-k % tabla de diferencias divididas
DD(m,k) = (DD(m + 1,k - 1) - DD(m,k - 1))/(x(m + k - 1)- x(m));
end
end
a = DD(1,:);
n = a(N+1);
for k = N:-1:1
n = [n a(k)] - [0 n*x(k)]; %n(x)*(x - x(k - 1))+a_k - 1
end

Ejemplo: supongamos que queremos encontrar un polinomio de Newton con los siguientes
puntos de datos {(-2, -6), (-1, 0), (1, 0), (2, 6), (4, 60)}
A partir de estos puntos de datos, construimos la tabla de diferencias divididos de la Tabla
siguiente y luego utilizar esta tabla, para obtener el polinomio de Newton

En Command Window:
%do_newtonp.m
x = [-2 -1 1 2 4]; y = [-6 0 0 6 60];
n = newtonp(x,y) %l = lagranp(x,y) para comparación
x = [-1 -2 1 2 4]; y = [ 0 -6 0 6 60];
n1 = newtonp(x,y) % con el orden de los datos modificados para la comparación
xx = [-2:0.02: 2]; yy = polyval(n,xx);
clf, plot(xx,yy,'b-',x,y,'*')

Otro ejemplo:(Nieves Hurtado. Métodos numéricos aplicados a la ingeniería. Pág. 382-384)

x=[-2 -1 0 2 3 6];
fx=[-18 -5 -2 -2 7 142];
M=6;N=M-1;
for i=1:N
T(i,1)=(fx(i+1)-fx(i))/(x(i+1)-x(i));
end
for j=2:N
for i=j:N
T(i,j)=(T(i,j-1)-T(i-1,j-1))/(x(i+1)-x(i-j+1));
end
end
T
%Respuesta:
T =
13 0 0 0 0
3 -5 0 0 0
0 -1 1 0 0
9 3 1 0 0
45 9 1 0 0

Hicimos el programa "do_newtonp1.m" MATLAB para hacer este trabajo y trazar la gráfica de
las funciones polinómicas, junto con la gráfica de la función verdadera f(x) y sus funciones de
error por separado para la comparación como se representa en la Fig., donde se han omitido las
partes para n8 (x) y n10 (x) para proporcionar a los lectores algo de espacio para la práctica.
x = [-1 -0.5 0 0.5 1.0]; y = (1./(1+8*x.^2)); % f(x) es la función verdadera.
n = newtonp(x,y)
xx = [-1:0.02: 1]; % el intervalo superior
yy =spline(x,y,xx); % gráfica de la función verdadera
yy1 = polyval(n,xx); % gráfica de la función polinómica aproximada
subplot(221), plot(xx,yy,'k-', x,y,'o', xx,yy1,'b')
subplot(222), plot(xx,yy1-yy,'r') % gráfico de la función de error

n=
2.3704 -0.0000 -3.2593 -0.0000 1.0000
Otro programa: polinomioNEWTON.m

clear;clc
x=[1 2 3 4 5 6]; y=[1 3 -1 0 3 2]; % entrada de datos
%x =[-4 -2 0 2 4]; y =[0.7568 -0.9093 0 0.9093 -0.7568];

xa=x;ya=y;
% Formacion de las diferencias divididas
d=zeros(length(y));
d(:,1)=y';
for k=2:length(x)
for j=1:length(x)+1-k
d(j,k)=(d(j+1,k-1)-d(j,k-1))/(x(j+k-1)-x(j));
end
end
% Formacion del polinomio
for w=1:length(x)
ds=num2str(abs(d(1,w)));
if w>1
if x(w-1)<0
sg1='+';
else
sg1='-';
end
end
if d(1,w)<0
sg2='-';
else
sg2='+';
end
if w==1
acum=num2str(d(1,1));
elseif w==2
polact=['(x' sg1 num2str(abs(x(w-1))) ')' ];
actual=[ds '*' polact];
acum=[acum sg2 actual];
else
polact=[polact '.*' '(x' sg1 num2str(abs(x(w-1))) ')' ];
actual=[ds '*' polact];
acum=[acum sg2 actual];
end
end

% Presentacion de resultados
fprintf(' n Valores de X y Y n ');
disp(xa);
disp(ya);
fprintf('\n Polinomio interpolación Newton : %s \n',acum);
x=input(' X interp = ');
if x>max(xa)|x<min(xa)
fprintf('\nPunto fuera de rango. El resultado puede ser equivocado \n');
end
xinterp=x;
yinterp=eval(acum);
fprintf(' Y(%g) = %g \n',x,yinterp);
% Grafica de los puntos
fprintf(' Pulse cualquier tecla para ver la grafica de los puntos\n');
pause
xg=linspace(min(xa),max(xa));
x=xg;yg=eval(acum);
plot(xg,yg,xa,ya,'.r',xinterp,yinterp,'or');
grid on

Función polyfit de MATLAB

MATLAB se refiere a la guarnición con un polinomio como curva de "regresión polinómica."
La función polyfit devuelve un vector de m + 1 coeficientes que representan el polinomio de
ajuste óptimo de grado m para el (xi, yi) un conjunto de puntos de datos. El orden coeficiente
corresponde a la disminución de potencias de x, es decir,
y c =a1 x m +a 2 x m−1++ a3 x m−2+ …+am x +am +1
El vector A = (a1, a2, ..., am, am + 1) que contiene los coeficientes se obtiene mediante la ejecución
de A = polyfit (X, Y, m), donde X = (x1, x2, ..., xn) e Y = (y1, y2, ..., yn).
MATLAB proporciona la polyval (a, x) función para calcular los valores aproximados y c,i en
(x1, x2, ..., xn) para los coeficientes empotrados (a1, a2, ..., am + 1):
y c ,i=a1 x mi + a2 x m−1
i ++a3 x m−1
i + …+a m x i +a m+1
Para calcular el cuadrado medio del error (MSE), polyval se usa en el siguiente ejemplo.
Ejemplo: Se tiene los siguientes puntos de datos, encontrar el mejor ajuste polinomial de
segundo, tercero, cuarto y quinto grado. Como también calcular el cuadrado medio del error.
(x,y)=(-10,980),(-8,-620),(-6,-70),(-4,80),(-2,100),(0,90),(2,0),(4,-80),(6,-90),(8,10),(10,220)

% Este programa determina la mejor curva de ajuste polinomio pasa

% a través de un conjunto dado de puntos de datos para polinomio
% de grados 2 a 5.
clear; clc;
% Datos suministrados
X = [ -10 -8 -6 -4 -2 0 2 4 6 8 10 ];
Y = [-980 -620 -70 80 100 90 0 -80 -90 10 220 ];
% Ajustar a los polinomios de grado m:
m=2:5;
% Crear un segundo conjunto de puntos x en el mismo intervalo
% con el fin de trazar la curva ajustada
X2 = -10:0.5:10;
% Vector para realizar un seguimiento de MSE (para imprimir al final)
MSE=zeros(length(m));
for i=1:length(m)
fprintf('m = %i \n',m(i));
% un polinomio de orden de grado m tiene m+1 coeficientes:
A = zeros(m(i)+1);
% hacer el ajuste
A = polyfit(X,Y,m(i));
% calcular el MSE (cuadrado medio del error)
Yc = polyval(A,X);
MSE(i) = (1/length(X)) * sum((Y-Yc).^2);
% imprimir tabla
fprintf(' x y yc \n');
fprintf('-----------------------------\n');
for j=1:length(X)
fprintf('%5.0f %5.0f %8.2f \n', ...
X(j),Y(j),Yc(j));
end
fprintf('\n\n');
% plot
subplot(2,2,i);
plot(X2,polyval(A,X2),X,Y,'o');
xlabel('x'); ylabel('y'); grid;
axis([ -10 10 -1500 500]);
title(sprintf(' Ajuste polinomial %d grado ',m(i)));
end
% impresión de MSE calculado
fprintf(' m MSE \n')
fprintf('-------------------\n');
for i=1:length(m)
fprintf(' %g %10.2f \n',m(i),MSE(i))
end

Salida del programa:

m MSE
------------------
2 32842.37
3 2659.97
4 2342.11
5 1502.95

El MSE dismuye a medida que se incrementa el ajuste de polinomio

Aproximación Funcional e Interpolación:

Ejemplo: A continuación, se presentan las presiones de vapor del cloruro de magnesio:
Puntos: 0 1 2 3 4 5 6 7
P (mm Hg): 10 20 40 60 100 200 400 700
T (°C): 930 988 1050 1088 1142 1316 1223 1418
Calcule la presión de vapor correspondiente a T = 1000 °C.

Solución:
A continuación, se muestra la forma de resolver este ejemplo con MATLAB. Para el cálculo
de Pint se aplica el método de Lagrange. El último renglón de código ingresado muestra la
forma de usar la interpolación con spline cúbico.
>> T=[930 988 1050 1088 1142 1316 1223 1418];
>> P=[10 20 40 60 100 200 400 700];
>> Tint=1000;
>> Pint=interp1(T,P,Tint)
Pint =
23.8710
>> PintCS=interp1(T,P,Tint,'spline')
PintCS =
22.1490

La función interp1(T,P,Tint) realiza la interpolación lineal, mientras que

interp1(T,P,Tint,'spline') realiza una interpolación con spline cúbico.

Ejemplo: La siguiente tabla muestra las viscosidades del isopentano a 59 °F y a diferentes

presiones.
Presión (psia) Viscosidad (centipoisese)
426.690 2468
483.297 2482
497.805 2483
568.920 2498
995.610 2584
1422.300 2672
2133.450 2811
3555.750 3094
2466.900 3236
7111.500 3807
Elabore un programa para aproximar el valor de la viscosidad a las presiones de 355.575,
711.150, 2844.600, 5689.200 y 8533.801 psia, utilizando la aproximación cúbica segmentaria
de Bessel.

Solución:
% Aproximación spline cúbica segmentaria de Bessel.
P=[426.690 483.297 497.805 568.920 995.610 1422.300 2133.450 3555.750 ...
2466.900 7111.500];
Visc=[2468 2482 2483 2498 2584 2672 2811 3094 3236 3807];
PP=[355.575 711.150 2844.600 5689.200 8533.801]; % La mitad del número % de elementos
de los vectores.
VV=spline(P,Visc,PP);
for i=1:5
fprintf(' Visc(%8.2f) = %8.2f\n',PP(i),VV(i))
end
% Los vectores deben tener el mismo número de elementos, el primero y
% el último valor en Y se utilizan como las pistas de extremo para el
% spline cúbico.

>> viscosidad
Visc( 355.57) = 2350.20
Visc( 711.15) = 2529.30
Visc( 2844.60) = 3479.55
Visc( 5689.20) = 644.29
Visc( 8533.80) = 16288.23

Ejemplo: El calor específico Cp (cal/K gmol) del Mn3O4 varía con la temperatura de acuerdo
con la siguiente tabla:
T (K): 280 650 1000 1200 1500 1700
Cp (cal/K gmol): 32.7 45.4 52.15 53.7 52.9 50.3
Aproxime esta información con un polinomio por el método de mínimos cuadrados.

Solución:
…
Ejemplo: Use la aproximación polinomial de segundo grado obtenida en el ejemplo anterior
para aproximar el calor específico del Mn3O4 a una temperatura de 800 K.
Solución:
T =[280 650 1000 1200 1500 1700];
Cp =[32.7 45.4 52.15 53.7 52.9 50.3];
a=polyfit(T,Cp,2);
fprintf('a0= %8.5f a1= %9.6f a2= %9.6f\n',a(3),a(2),a(1))
Tint=800;
Cpint=a(3)+a(2)*Tint+a(1)*Tint^2;
fprintf(' Cp(%4.0f) = %6.2f\n',Tint,Cpint);

>> cp_interpola
a0=19.29544 a1= 0.053728 a2=-0.000021
Cp( 800) = 48.92

Ejemplo: Dada la tabla

Puntos: 0 1 2 3 4
xi: 1.00 1.35 1.70 1.90 3.00
f(xi ): 0.00000 0.30010 0.53063 0.64185 1.09861

Construya una tabla de diferencias divididas para aproximar f(x) en x = 1.50; utilice un
polinomio de Newton de segundo grado.

Solución:
X=[1.35 1.70 1.90];
FX=[0.30010 0.53063 0.64185];
Xint=1.50;
T(1,1)=(FX(2)-FX(1))/(X(2)-X(1));
T(2,1)=(FX(3)-FX(2))/(X(3)-X(2));
T(2,2)=(T(2,1)-T(1,1))/(X(3)-X(1));
Fxint=FX(1)+T(1,1)*(Xint-X(1))+T(2,2)*(Xint-X(1))*(Xint-X)

Fxint =
0.3947 0.4045 0.4101 %La respuesta es 0.4045

Interpolación spline cúbica con condiciones externas (csape).

Determina el spline cúbico, natural o sujeto para un conjunto de datos (x,y).
csape(x,y,'variational') natural; csape(x,y, 'clamped'),[a,b] sujeto
f'(x0)=a, f'(xn)=b
Ejemplo: Los trazadores se usan a menudo en software de gráfica para trazar una curva suave a
través de un conjunto de puntos de datos. En la siguiente tabla se proporcionan puntos de cierta
curva en coordenadas (x,y).
i 1 2 3 4 5 6 7 8 9 10
xi 0.15 0.20 0.40 0.50 0.50 0.30 0.30 0.40 0.25 0.15
yi 0.50 0.80 1.00 1.00 0.80 0.60 0.40 0.30 0.20 0.50
y no es función de x ni viceversa, sin embargo, es posible generar una curva suave haciendo que
x y y dependan de un parámetro común t cuyos valores para el conjunto de datos son ti=i. Por
medio de trazadores naturales ajustar x y y a t, para valores de t desde 1 a 10 en incrementos de
0.05 luego grafique.

Solución:
x=[0.15 0.20 0.40 0.50 0.50 0.30 0.30 0.40 0.25 0.15];
y=[0.50 0.80 1.00 1.00 0.80 0.60 0.40 0.30 0.20 0.50];
t=1:10;s1=csape(t,x,'variational');s2=csape(t,y,'variational');
tt=1:0.05:10;xx=ppval(s1,tt);yy=ppval(s2,tt);
plot(x,y,'o',xx,yy),grid on

Ejemplo: En una reacción gaseosa de expansión a volumen constante, se observa que la presión
del reactor aumenta con el tiempo de reacción según se muestra en la tabla.
P (atm) 1.0000 1.0631 1.2097 1.3875 1.7232 2.0000 2.9100

t (min) 0.0 0.1 0.3 0.5 0.8 1.0 1.5

a) Determinar el mejor interpolante polinómico, de grado a lo más 3, para ajustar los datos
aplicando el método de los mínimos cuadrados.
b) Graficar los datos y la curva de ajuste.
c) Determine P para t=0.7 minutos.
Solución:
% a)
t=[0, 0.1, 0.3, 0.5, 0.8, 1, 1.5];
P=[1, 1.0631, 1.2097, 1.3875, 1.7232, 2, 2.91];
p1=polyfit(t,P,1); %Polinomio Lineal
p2=polyfit(t,P,2); p3=polyfit(t,P,3);
E1=sum((P-polyval(p1,t)).^2); %Error de p1
E2=sum((P-polyval(p2,t)).^2);
E3=sum((P-polyval(p3,t)).^2);
[E1; E2; E3] %El mejor ajuste el de menor error
ans =
0.096983631503322
0.000415487678758
0.000000031060139

% b)
tt=0:0.005:1.5; Pp3=polyval(p3,tt); plot(t,P,'o',tt,Pp3); grid
xlabel('Tiempo(min)'), ylabel('Presión(atm)'), legend('Datos', 'Ajuste Cúbico')
>> pi=polyval(p3,0.7)
pi =
1.601416680484238

Ejemplo: Construyamos mediante splines cúbicas el contorno de la citroneta de la fotografía:

Solución:

i) Primero, se realiza la toma de datos usando una malla cuadriculada sobre la

fotografía:

ii) Luego, los datos se tabulan. Los datos del contorno superior son:
 k : 0 1 2 3 4 5 6 7 8 9 10 11 12
xk: 2 2.7 3.8 6 8 10 13 16 18 21 25 30 36
yk: 5 7.8 9 10 10.2 10.3 10.4 14.5 15 15.4 15.5 14 5

iii) Ahora se calculan los elementos de la matriz, considerando las definiciones,hasta k=3, como
empieza de cero, debería realizarse hasta 11:
hk = xk+1 − xk
λk =(yk+1 – yk) / hk
μk = 3(λk − λk−1)

k 0 1 2
3
hk 0.7 1.1 2.2 2
λk 2.8/0.7 = 4 1.2/1.1 = 12/11 1/2.2 = 5/11 0.2/2 = 1/10
μk 3 × ( 12/11 − 4) = −96/11 3 × ( 5/11 – 12/11) = −21/11 3×(1/10–5/11) =
−117/110

iv) Ahora graficando con Matlab.

Xk=[ 2 2.7 3.8 6 8 10 13 16 18 21 25 30 36];
Yk= [5 7.8 9 10 10.2 10.3 10.4 14.5 15 15.4 15.5 14 5]; %(Xk,Yk) puntos de la Curva
xk1=2:0.01:36; % discretización de [2,36]
spl1=spline(Xk,Yk,xk1);
plot(xk1,spl1)
axis([0 37 0 17])

Tarea: Realice trazadores cúbicos sujetos o splines cúbicas de contorno de la imagen de una
fruta u hotaliza.
Comparación de interpolación:
Observamos que obtenemos la misma respuesta con la interpolación lineal, sino una respuesta
diferente a la interpolación cúbica. Aquí la interpolación cúbica hace mucho mejor trabajo de
montaje. Usted puede ver por qué en esta figura. La interpolación cúbica es mucho más cerca de
la verdadera función de la mayor parte del rango. Aunque tenga en cuenta que sobreestima en
algunas regiones y subestima en otras regiones.
t = [0.5 1 3 6];
f = [0.6065 0.3679 0.0498 0.0025];
figure
hold all
x = linspace(0.25, 6);
plot(t,f,'ko ')
plot(x, interp1(t,f,x),'r-')
plot(x, interp1(t,f,x,'pchip'),'b-')
plot(x,exp(-x),'k--')
xlabel('t')
ylabel('f(t)')
legend('datos en bruto (nodo)','interpolación lineal','interpolación cúbica','exp(-t)')
% exp(-t) es la función verdadera.
% Polinomio de interpolación cúbica de Hermite a trozos(PCHIP:
% Piecewise Cubic Hermite Interpolating Polynomial).

Bibliografía:
http://www.sc.ehu.es/sbweb/energias-renovables/MATLAB/numerico/diferencial/diferencial.html