Factorizando la matriz de van der Monde¶
Apartado 1¶
Define una función que
- Recibe como argumentos:
- Una lista
xsde nodos connelementos $x_1,\dots,x_n$. - Un entero
kmenor quen
- Una lista
- Devuelve un array
Vde dimensiones $n\times k$ $$ {\displaystyle V={ \begin{bmatrix} 1&x _{1}&x _{1}^{2}&\dots &x _{1}^{k-1}\\ 1&x _{2}&x _{2}^{2}&\dots &x _{2}^{k-1}\\ \vdots &\vdots &\vdots &\ddots &\vdots \\ 1&x _{n}&x _{n}^{2}&\dots &x _{n}^{k-1} \end{bmatrix}}.} $$
In [16]:
function V = vdm(xs, k)
% Matriz de van der Monde dada una lista de puntos xi
%
% Tcrea una matriz NxN cuyo elemento M[i,j] es xi^j
%
% xs: lista de puntos
%
% Devuelve un array NxN
N = length(xs);
M = zeros(N, N);
for i = [1:N]
for j = [1:k]
V(i,j) = xs(i)^(j-1);
end
end
end
Otra versión vectorizada (no parece ser más rápida que la anterior)
In [17]:
function V = vdm2(xs,k)
% Matriz de van der Monde dada una lista de puntos xi
%
% Tcrea una matriz NxN cuyo elemento M[i,j] es xi^j
%
% xs: lista de puntos
%
% Devuelve un array NxN
N = length(xs);
M = zeros(N, N);
for j = [1:k]
V(:,j) = xs.^(j-1);
end
end
Otra versión doblemente vectorizada (no parece ser más rápida que la anterior)
In [36]:
function V = vdm3(xs, k)
% Matriz de van der Monde dada una lista de puntos xi
%
% Tcrea una matriz NxN cuyo elemento M[i,j] es xi^j
%
% xs: lista de puntos
%
% Devuelve un array NxN
V = xs'.^(0:(k-1));
end
In [37]:
xs = [-1,0,1,2,3];
vdm(xs,3)
vdm2(xs,3)
vdm3(xs,3)
In [ ]:
xs = [1:100];
callV = @()(vdm(xs, 30));
timeit(callV)
callV = @()(vdm2(xs, 30));
timeit(callV)
callV = @()(vdm3(xs, 30));
timeit(callV)
Apartado 2¶
La matriz $V$ está relacionada con el polinomio interpolador y con el polinomio aproximador.
En particular, los coeficientes del polinomio aproximador son la solución de $$ \left(V^T\cdot V\right)\cdot \mathbf{a} = V^T\cdot \mathbf{y}, $$
Se pide:
- Discuta qué factorizaciones admiten las matrices $V$ y $V^T\cdot V$, que puedan ser de ayuda para calcular polinomios interpoladores o aproximadores.
- Calcule estas factorizaciones para una lista de puntos
[1,2 ... n], dondenes un entero.
- La matriz V sólo es cuadrada cuando el grado del polinomio es igual al número de puntos menos 1. En este caso, hemos visto en clase de teoría que la matriz V es invertible siempre que los números $x_j$ sean todos distintos. Como V no es simétrica ni banda, sólo podemos aplicar una factorización V=PLU (o PV=LU, que es casi lo mismo), cuando V es cuadrada.
In [26]:
n = 5;
xs = linspace(1,n,n);
V = vdm3(xs,length(xs))
[L, U, P] = lu(V)
- El método
lutambién permite factorizar matrices que no son cuadradas, aunque no lo hemos estudiado en clase...
In [27]:
n = 5;
k = 3;
xs = linspace(1,n,n);
V = vdm3(xs,k)
[L, U, P] = lu(V)
- La matriz $V^T\cdot V$ es simétrica y definida positiva, así que admite una factorización de Cholesky.
- Comprobamos que hemos obtenido una matriz $L$ tal que $L^T\cdot L=V^T\cdot V$, salvo error numérico (observamos que $L^T\cdot L=V^T\cdot V$, a pesar de que $L$ y $V$ tienen distintas dimensiones).
In [30]:
k = 3;
V = vdm3(xs,k);
VTV = V'*V
L = chol(VTV);
L'*L
VTV - L'*L
Apartado 3¶
Utilice las factorizaciones del apartado anterior para encontrar los polinomios aproximadores de grados 2 y 4 con la siguiente serie de datos:
xs = [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19, 20]
ys = [ 0.071, 0.458, 0.908, 1.344, 1.641, 1.934, 1.451, 1.571,
1.282, 0.446, 0.261, -0.183, -0.372, -0.162, -0.181, -0.12 ,
-0.545, -0.114, 0.505, 0.474]
- Utilizamos la factorización de Cholesky $$ V^T\cdot V\cdot \mathbf{a} =L^T\cdot L\cdot \mathbf{a} = V^T\cdot \mathbf{y} $$ Ahora es suficiente con resolver dos sistemas de ecuaciones con matriz de coeficientes triangular: $$ L^T\cdot \mathbf{z} = V^T\cdot \mathbf{y} $$ $$ L\cdot \mathbf{a} = \mathbf{z} $$
In [45]:
xs = [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20];
ys = [ 0.071, 0.458, 0.908, 1.344, 1.641, 1.934, 1.451, 1.571, 1.282, 0.446, 0.261, -0.183, -0.372, -0.162, -0.181, -0.12 , -0.545, -0.114, 0.505, 0.474];
n = length(xs);
In [46]:
hold on;
deg = 4;
legends = cell(1,deg + 1);
l = 1;
for k = [1:deg]
%para ajustar un polinomio de grado k, llamamos a la función
% vdm con argumento k+1
V = vdm3(xs,k+1);
VTV = V'*V;
L = chol(VTV);
z = (L')\(V'*ys');
a = L\z;
x_eval = linspace(1,n,200);
y_eval = polyval(flip(a), x_eval);
plot(x_eval, y_eval)
legends{l} = sprintf('grado %d',k);
l = l + 1;
end
plot(xs, ys, '*')
legends{l} = 'puntos';
legend(legends);
hold off;
title('Polinomio aproximador')
- Comparamos el resultado con el obtenido con
polyfit
In [47]:
hold on;
deg = 4;
legends = cell(1,deg + 1);
l = 1;
for k = [1:deg]
cs = polyfit(xs, ys, k);
x_eval = linspace(1,n,200);
y_eval = polyval(cs, x_eval);
plot(x_eval, y_eval)
legends{l} = sprintf('grado %d',k);
l = l + 1;
end
plot(xs, ys, '*')
legends{l} = 'puntos';
legend(legends);
hold off;
title('Polinomio aproximador')
In [ ]: