Condicionamento de Matrizes e Aproximação de Funções

Vista previa del material en texto

Elementos de Matemática Aplicada 2011 1
Laboratorio 3: Condicionamiento de matrices
1. Una matriz molesta
Considere la siguiente matriz
H5 =

1 1/2 1/3 1/4 1/5
1/2 1/3 1/4 1/5 1/6
1/3 1/4 1/5 1/6 1/7
1/4 1/5 1/6 1/7 1/8
1/5 1/6 1/7 1/8 1/9

una matriz de apariencia inocente, simétrica y definida positiva. ¿Dije defi-
nida positiva? Calculen los autovalores para verificar.
¿Cuál es el condicionamiento de esta matriz?
¿Eso es malo?
Al final, tan buenita no era. Vamos a experimentar un poco. Para meter
estas matrices en MatLab hay un comando hilb(5) o hilb(10) para armar
una similar pero de 10× 10, o hilb(17) para 17× 17 y aśı.
Vamos a hacer un cálculo sencillo: calculemos la inversa de estas matrices,
y multipliquemos la inversa por la matriz original ¿Qué tiene que dar?
En MatLab, para calcular la inversa, está el comando inv. Por ejemplo
H5=hilb(5);
H5*inv(H5)
funciona bastante bien
Ahora hagan lo mismo para la matriz de 10× 10... Bueno, ya no da tan
bien, pero cuatro d́ıgitos de presición son suficientes para varias cosas.
Ahora de nuevo, con la de 15× 15...
¿Qué paso?
Parece que estas matrices están tan mal condicionadas que no permiten
ni siquiera calcular la inversa de manera correcta.
Otro śıntoma, es que la factoriación de Choleski falla para la matriz de
15× 15, porque los errores de redondo hacen que el autovalor mas pequeño
resulte negativo cuando debeŕıa ser positivo (chequeenlo).
Elementos de Matemática Aplicada 2011 2
Aśı que śı, son unas matrices muy molestas. Más valdŕıa evitarlas. Porque
después de todo ¿a quién se le podŕıa ocurrir armar un sistema de ecuaciones
con una de estas matrices?
2. Un problema de aplicación
Si queremos aproximar una función f mediante un polinomio P , una
posibilidad que tenemos a nuestra disposición es el polinomio de Taylor.
Esto nos va a dar una buena aproximación cerca de un punto.
¿Y si queremos una aproximación que sirva en todo un intervalo? Por
ejemplo en el intervalo (0, 1).
Tomense un rato para evaluar posibles alternativas antes de seguir ade-
lante.
La posibilidad que a mı́ se me ocurre es la siguiente: calcular el polinomio
de Taylor en el punto medio del inervalo.
Veamos qué tal funciona. Probemos con la función
f(x) = e10x
Hagamos el polinomio de Taylor de grado 4. Calcule las primeas 4 deri-
vadas, y grafique la función y el polinomio. Un consejo, para que se vea algo
en el gráfico va a hacer flata cambiar la escala de los ejes a mano mediante
el comando axis. Por ejemplo
h=0.01;
x=0:h:1;
f=exp(10*x);
P=exp(5)*(1+10*(x-.5)+50*(x-.5).^2+1000/6*(x-.5).^3+10000/24*(x-.5).^4);
plot(x,f,’b’,x,P,’--g’);
axis([0 1 -1000 24000]);
¿Qué tal quedó el gráfico? Como era de esperar, seguro que aproxima
bastante bien en el centro del intervalo, y bastante mal en los borde. ¿No
habrá alguna forma de mejorar esto?
Una primera idea (que en general funciona) es aumentarle el grado al
polinomio... Pero ¿no habrá alguna manera de mejorarlo con otro polinomio
de grado 4?
Bueno, pensemos, Taylor se basa en derivadas, que dan información cerca
del punto donde estoy parado. Si quiero tener algo que funcione bien en todo
un intervalo necesitamos algo que integre toda la información del intervalo.
¿Ya se les ocurrió que usar?
Lo que queremos es que P − f sea chico.
Usar
∫ 1
0 P − f no es buena idea, porque en esa integral los errores nega-
tivos se cancelan con los errores positivos.
Elementos de Matemática Aplicada 2011 3
La integral
∫ 1
0 |P − f | es una idea mejor. Podŕıamos intentar buscar P
para que esta integral sea mı́mina. El problema es que el valor absoluto es
dif́ıcil de derivar.
Aśı que pasamos a ∫ 1
0
(P (x)− f(x))2 dx
y queremos hallar el mı́nimo de esta cosa. ¿Mı́nimo respecto a qué variables?
Escribamos lo mismo de manera un poco más expĺıcita. Hacemos P (x) =
c1 + c2x + c3x
2 + c4x
3 + c5x
4 y queda
F (c1, c2, c3, c4, c5) =
∫ 1
0
(
c1 + c2x + c3x
2 + c4x
3 + c5x
4 − f(x)
)2
dx
Ahora lo que queremos es el mı́nimo de esto. Para eso, derivamos respecto
a ci e igualamos a cero. Háganlo, quedan cosas como∫ 1
0
(c1 + c2x + c3x
2 + c4x
3 + c5x
4)dx =
∫ 1
0
f(x)dx
∫ 1
0
(c1x + c2x
2 + c3x
3 + c4x
4 + c5x
5)dx =
∫ 1
0
xf(x)dx
y aśı las otras.
Del lado izquierdo, reparta las integrales, calcule las
∫ 1
0 x
k = 1/(k + 1),
aśı queda un sistema de ecuaciones lineales para calcular los coeficientes ci.
Y la matŕız de este sistema es...
A pesar de que haya quedado una matriz mal condicionada, podemos
seguir adelante, y resover el sistema de ecuaciones. Hay que calcular los
términos independientes, que son integrales no muy complicadas
∫ 1
0
f(x) =
∫ 1
0
e10x =
1
10
e10x
∣∣∣∣1
0
= (e10 − 1)/10∫ 1
0
xe10x = x
1
10
e10x
∣∣∣∣1
0
−
∫ 1
0
1
10
f(x) = e10/10− (e10 − 1)/100∫ 1
0
x2e10x =< int. por partes, bla, bla >= e10/10− 2e10/100 + 2(e10 − 1)/1000
y aśı siguiendo.
Ahora resuelva el sistema y arme el nuevo polinomio. Compare la gráfica
de este nuevo polinomio con la de la función original, si todo salió bien, la
aproximación debeŕıa ser bastante mejor que la de Taylor.
Elementos de Matemática Aplicada 2011 4
3. Algunos desaf́ıos
Para los que quieran seguir adelante con este problema, acá van algunos
desaf́ıos
Podemos comparar las dos formas de aproximar para polinomios de
grado más grande. En estos casos ambas aproximaciones funcionan
bastante bien, aśı que no alcanza un gráfico para notar cuál aproxima-
ción es mejor. Piensen alguna manera de controlar qué aproximación
funciona mejor.
Al aumentar más el grado del polinomio, el condicionamiento de la
matriz de Hilbert va empeorando, llega un momento en que los errores
de redondeo son tan grandes que hacen que la aproximación por Taylor
pase a ser mejor. Busquen para qué grados pasa esto.
Una forma de intentar zafar de los errores de redondeo es usar el
siguiente programa que calcula la matriz inversa de Hilbert.
function HI=hilbinv(n)
k=n;
for i=1:n
if i>1
k=((n-i+1)*k*(n+i-1)) / (i-1)^2;
end
l=k*k;
HI(i,i)=l/(2*i-1);
for j=(i+1):n
l=-((n-j+1)*l*(n+j-1)) / (j-1)^2;
HI(j,i)= l/(i+j-1);
HI(i,j)= HI(j,i);
end
end
¿Alguien se anima a investigar por qué este programa funciona?
Pruebe a ver si con este programa puede mejorarse los cálculos de
los polinomios en los casos en que los errores de redondeo causaban
problemas.