Diferencias divididas de sumas y productos¶
Apartado 1¶
Crea una función que:
- Recibe como argumentos dos listas de valores:
- una lista de ordenadas
xs, de longitud n - otra lista
ysde la misma longitud (que contiene el resultado de evaluar una cierta función $f$ en cada punto dexs)
- una lista de ordenadas
- Devuelve un array
Tfde dimensiones $n\times n$ con todas las diferencias divididas del tipo $f[x_{\nu },\ldots ,x_{{\nu +j}}]$, almacenadas enTfde la forma siguiente:- En la entrada $(i,j)$, si $i>j$, dejamos $T_f[i,j]=0$.
- En la diagonal principal aparecen los valores $T_f[j,j]=f(x_j)=y_j$.
- Sobre la diagonal $i<j$ aparecen diferencias divididas en puntos consecutivos: $$ T_f[i,j] = f[x_i,x_{i+1},\dots,x_j] $$
En otras palabras $$T_f=\left(\begin{array}{ccccc} f [x_0] & f [x_0, x_1] & f [x_0, x_1, x_2] & \cdots & f [x_0, x_1, x_2, \ldots, x_n]\\ & f [x_1] & f [x_1, x_2] & \cdots & f [x_1, x_2, \ldots, x_n]\\ & & f [x_2] & & f [x_2, \ldots, x_n]\\ & & & \ddots & \vdots\\ & & & & f [x_n] \end{array}\right)$$
In [3]:
function M = DD(xs, ys)
n = length(xs);
M = zeros(n);
%Primero rellenamos la diagonal
for i =[1:n]
M(i,i) = ys(i);
end
%Ahora construimos diferencias divididas usando la fórmula
% recursiva
for j = [1:n]
for i = [1:n-j]
M(i,i+j) = (M(i+1,i+j) - M(i,i+j-1))/(xs(i+j)-xs(i));
end
end
end
In [4]:
% Ejemplo del Burden Faires capitulo 3.3
xs = linspace(1,2.2,5);
ys = [0.7651977, 0.6200860, 0.4554022, 0.2818186, 0.1103623];
DD(xs, ys)
- Otro ejemplo: las diferencias divididas de un polinomio de grado 2 son no nulas si tienen a lo sumo tres puntos, ya que $$f[x_i,x_{i+1},x_{i+2}] = \frac{f''(\xi)}{2!}$$
- para un punto $\xi$ desconocido.
In [5]:
xs = [0:4]
p = @(x)(x.^2 + x + 2);
ys = p(xs);
DD(xs, ys)
Apartado 2¶
- Dadas dos funciones $f$ y $g$, y una lista de ordenadas $x_0,\dots,x_n$, creamos las matrices $T_f$ y $T_g$ definidas antes.
- Creamos también la matriz $T_{fg}$ que corresponde a la función producto $x\rightarrow f(x)g(x)$ y la matriz $T_{f+g}$ que corresponde a la función suma $x\rightarrow f(x) + g(x)$.
- Afirmamos que $T_f + T_g = T_{f + g}$.
- Afirmamos que $T_f\cdot T_g = T_{fg}$.
Comprueba empíricamente las dos afirmaciones anteriores.
In [7]:
%Tomamos dos funciones cualesquiera, y construimos sus matrices Tf
f = @(x)(2.^x);
xs = linspace(0,3,4);
Tf = DD(xs, f(xs))
g = @(x)(x.^2);
xs = linspace(0,3,4);
Tg = DD(xs, g(xs))
- Comprobamos que $T_f + T_g = T_{f + g}$
In [8]:
Tf + Tg
Tg + Tf
DD(xs, g(xs) + f(xs))
- Comprobamos que $T_f\cdot T_g = T_{fg}$
In [9]:
Tf * Tg
Tg * Tf
DD(xs, g(xs) .* f(xs))
Apartado 3¶
Para un polinomio: $$ P(x) = a_0+a_1 x +\dots+a_n x^n $$ las dos propiedades anteriores implican $$ T_P = a_0\cdot Id + a_1\cdot T_x+a_2\cdot (T_x)^2+\dots+a_n\cdot (T_x)^n. $$
- Utiliza esta identidad para encontrar los coeficientes del polinomio $Q$ que interpola a otro polinomio $P$ dado por la expresión $$ P(x) = 1 + 2 x + 3x^2 + 4x^3 $$ en los puntos $x_1=1, x_2=2, x_2=3$.
- Comprueba tu solución construyendo el polinomio interpolador por otra técnica diferente.
$$ P(x)=f[x_{0}]+f[x_{0},x_{1}](x-x_{0})+\cdots +f[x_{0},\ldots ,x_{n}](x-x_{0})(x-x_{1})\cdots (x-x_{{n-1}}). $$
- Las diferencias divididas de la forma $f[x_{0},\ldots ,x_{k}] $ son todo lo que necesitamos para construir el polinomio interpolador en forma de Newton:
- En la matriz anterior, las diferencias que necesitamos aparecen en la primera fila.
In [32]:
%P es un polinomio de grado 4
coefs_P = [1:4];
%Interpolamos en 4 puntos => Q es de grado 3
xs = linspace(-1,1,3);
Tx = DD(xs, xs);
d = length(coefs_P) - 1;
T_P = zeros(d);
for j=[0:d]
T_P = T_P + coefs_P(j+1)*(Tx^(d-j));
end
T_P
coefs_Q = T_P(1,:)
In [50]:
ys = polyval(coefs_P, xs);
T_P = DD(xs, ys)
coefs_Q = T_P(1,:)
- Usamos evaluación anidada para evaluar el polinomio interpolador $$P(x)=f[x_{0}] + (x-x_{0})\Big( f[x_{0},x_{1}]+(x-x_{1})\big( f[x_{0},x_{1},x_{2}] + \cdots +(x-x_{{k-1}})f[x_{0},\ldots ,x_{k}] \big)\Big) $$
In [52]:
function y_eval = eval_Newton(xs, coefs, x_eval)
n = length(xs);
neval = length(x_eval);
y_eval = zeros(1,neval);
%Evaluacion anidada:
for i=[0:n-1]
y_eval = y_eval.*(x_eval-xs(n-i)) + coefs(n-i);
end
end
In [53]:
x_eval = linspace(min(xs), max(xs) ,100);
y_P = polyval(coefs_P, x_eval);
y_Q = eval_Newton(xs, coefs_Q, x_eval);
hold on;
plot(xs, ys, 'o');
plot(x_eval, y_P);
plot(x_eval, y_Q);
hold off;
legend('puntos', 'polinomio original', 'polinomio interpolador')
In [58]:
x_eval = linspace(min(xs), max(xs) ,100);
y_P = polyval(coefs_P, x_eval);
y_Q = eval_Newton(xs, coefs_Q, x_eval);
coefs_polyfit = polyfit(xs, ys, length(xs)-1);
y_polyfit = polyval(coefs_polyfit, x_eval);
hold on;
plot(xs, ys, 'o');
plot(x_eval, y_P);
plot(x_eval, y_Q, 'g--');
plot(x_eval, y_polyfit, 'r:');
hold off;
legend('puntos', 'polinomio original', 'polinomio interpolador', 'polyfit')
