Pràctiques amb Octave: consultes
Consultes d'interés sobre les pràctiques d'Octave per l'Àlgebra lineal, Enginyeria Química de la UPC, curs 2003-4.
Les consultes s'aniran posant en una pàgina única per facilitar la cerca de paraules clau, ordenades per pràctiques:
Quan tot falla, és que ha arribat l'hora de mirar el manual d'instruccions.
(Pregunta del 8/1/2004) Tenim un dubte sobre la condicio.m, si escribim acció per acció a l'Octave, ens surten be els càlculs. En canvi quan apliquem la funció,diu que hi ha un error, cada ordre que li donem ha de seguir un ordre lògic? o potser al posar una vintena d'ordres es produeixen errors? hi ha alguna manera destalviar-se posar tot el carro d'ordres? a veure si axi ens queda tot una mica més clar.
Resposta (8/1/2004):
Sense saber quin missatge d'error us dona Octave fa de mal dir quin es el problema.
L'ordre de les instruccions en la funcio es el mateix en el que les executeu a Octave, la idea de la funcio es que fa exactament el mateix pero ho fa ella sola per a no haber de picar la sequencia cada cop.
1) El problema mes tipic amb funcions es error en la lectura o en el retorn de dades. La vostra funcio comenc,a per la linea
function [cond_1,cond_2,cond_inf]=condicio(A)
Aixo vol dir que la matriu de partida a la que heu de mirar els nombres de condicio es diu A, i que heu de definir unes variables anomenades exactament cond_1,cond_2,cond_inf abans de la linea final
return;
Si teniu aquest problema se us queixara de que li falta A, o que no sap que retornar com cond_1,2,inf
2) El segon problema mes tipic es error en cridar la funcio. Per aixo us vaig adjuntar l'exemple amb la descripcio de la funcio: vosaltres definiu una matriu fent per exemple:
A=[-1,2,0;1,3,5;-1,4,4];
Aleshores crideu a la funcio:
[a,b,c]=condicio(A)
I els nombres que la vostra funcio retorna com a cond_1, cond_2, cond_inf es guarden a a,b,c respectivament (i a mes els ensenya per pantalla).
3) Un altre problema habitual en els inicis es d'error en el format
de la funcio: cal que el fitxer que conte la funcio es digui exactament
condicio.m, que estigui al directori de treball de Octave (en sistemes
Windows sol ser C:\Archivos de Programa\GNU Octave 2.1.36\octave_files
) i que sigui un fitxer de texte que contingui exactament les instruccions
que us han funcionat en octave (compte amb el Word, perque fins i tot si
feu 'Exportar en format de texte' us pot colocar algun caracter sorpresa).
(dubte sobre el càlcul de les normes matricials i el nombre de condició, cap al 5/1/2004):
Resposta (5/1/2004):
Si, la frase
El nombre de condicio de A per cada norma es la norma de A multiplicada per la de la seva inversa inv(A).
vol dir que el nombre de condicio de A per cada norma es aquesta norma de la matriu A, multiplicada per aquesta mateixa norma de la matriu inversa de A.
Per tant la questio es pensar com calcular la norma 1,2, infinit d'una matriu, i aplicar el calcul a A i a inv(A).
Per a sumar les columnes de A: ho podeu fer amb un bucle doble, o podeu fer el seg"uent:
Si A es n x n, creeu un vector 1 x n de la forma [1,1,1,...,1] (tot de uns), aixo ho podeu fer aixi
[n,n]=size(A);
v=ones(1,n);
Acte seguit multipliqueu u=v*abs(A), us queda un vector fila que te components les sumes dels valors absoluts de cada columna.
Per a saber quin es el nombre mes gran en aquest vector, es simplement
n_inf=max(u);
Per la norma sub 1 podeu fer el mateix: creeu un vector columna (n x 1) tot de uns, posem que es diu w, multipliqueu z=abs(A)*w; i ara el vector z es un vector columna amb components les sumes dels valors absoluts de cada fila. mireu ara max(z) i teniu el maxim de la suma de les files.
Finalment, la norma sub 2: la instruccio eig es pot fer anar de dues maneres:
- la que varem explicar a classe: [C,D]=eig(A); retorna C= matriu de canvi de base de VEPs a canonica i D=forma diagonal de l'aplicacio.
- una altra: feu L=eig(A); i L es una llista dels VAPs de A. Aleshores max(abs(L)) es el VAP amb valor absolut (o modul, si es complex) maxim.
Ara el que suggereixo que feu es: probeu dins del programa Octave amb
aquestes instruccions (o d'alguna altra manera que us agradi mes, als manuals
mes la 'Breu
introduccio' hi surten mes ordres del llenguatge i variants...), i
quan us calculin el nombre de condicio que toca poseu aquestes instruccions
dins del fitxer condicio.m.
(dubte general sobre càlcul de formes de Jordan i el funcionament de la funció eig, cap al 30/12/2003):
Resposta (31/12/2003):
La funcio eig es va explicar a la practica que varem fer a l'aula d'ordenadors, es la funcio standard per diagonalitzar d'Octave. El problema que te es que no se n'entera de quan la matriu no diagonalitza, i en aquests casos dona solucions incorrectes.
Per fer la funcio jordan2 l'esquema suggerit es:
1) comenceu fent
[C0,D0]=eig(A);
per a que Octave 'diagonalitzi' la matriu.
2) Si els dos VAPs que surten son diferents aleshores A diagonalitzava i la solucio que ha trobat Octave es correcta. En aquest cas les solucions que heu de posar son
C=C0;
Jor=D0;
i ja esta.
3) En canvi si els dos VAPs han sortit iguals poden passar dues coses:
- o be el subespai de VEPs te dimensio 2, i aleshores la diagonalitzacio
de l'Octave es correcta, C=C0; Jor=D0;
- o be el subespai de VEPs te dimensio 1, A no diagonalitza, cal forma
de Jordan, i la solucio C0,D0 de Octave es incorrecta.
Per saber si esteu en el cas 2) o 3) heu de comparar els dos VAPs trobats a A. Com son els termes de la diagonal en D0 ho podeu fer aixi:
if (abs(D0(1,1)-D0(2,2))<tol)
% aqui va el que toca fer en el cas 2: retornar C=C0 i Jor=D0
else
% aqui va el que toca fer en el cas 3
endif;
(com Octave no calcula exactament sino amb arrodoniment enlloc de demanar que els 2 vaps siguin iguals demanem que la diferencia sigui mes petita que el marge d'error tol que ens hemprefixat).
Quan arribeu al cas 3: estem en el cas VAP doble, i cal decidir si hi han suficients VEPs o no. Per aixo en teniu prou amb mirar el rang de A-(VAP)Id. Per a fer-ho amb el marge d'error tol ho podeu fer aixi:
% el que toca fer en el cas 3:
if (rank(A-D0(1,1)*eye(2),tol)>0)
% esteu en el cas que no diagonalitza, les C0,D0 calculades per
Octave son incorrectes
% heu de trobar la base de Jordan
else
% esteu en el cas en que diagonalitza, i la solucio C0,D0 calculades
son correctes
% podeu retornar C=C0; Jor=D0;
end;
El que s'ha de fer en el cas en que no diagonalitza es el metode standard per calcular forma i base de Jordan 2x2, esta per exemple a l'exemple 5.48 i la proposicio 5.47 dels apunts de l'assignatura ( http://www-ma1.upc.es/~amoros/algquim.pdf ):
La forma de Jordan te el VAP doble D0(1,1) a la diagonal i un 1 a sota.
Per a trobar el primer vector de la base podeu probar per exemple amb u1=[1;0], i si resulta que aquest vector es VEP del VAP doble, aleshores podem assegurar per la teoria que u1=[0;1] si que serveix com a primer vector, i el vector u2 ve donat per A*u1-D0(1,1)*u1.
Si aneu malament de temps: en la funcio de Jordan limiteu-vos a detectar si esteu en el cas 1, 2 o 3 i escriviu per pantalla quin es el cas. Amb aixo no teniu la practica completa, pero si la majoria de punts.
La instruccio per escriure un texte per pantalla:
printf('Aquest texte sortira per pantalla');
(pregunta del 2004/1/9) Hem estat intentant fer el programa "quocient.m" i tenim alguns dubtes, un dels quals és com eliminar files o columnes de zeros d'una matriu, per aconseguir una base d'aquesta, per a qualsevol matriu esglaonada.
Per a aconseguir una base del subespai generat per les files o per les columnes d'una matriu A amb Octave no cal esglaonar: la instruccio
basecols=orth(A);
retorna una matriu basecols que te per columnes una base del subespai generat per les columnes de A. Per a trobar una base de les files podeu fer anar la mateixa instruccio pero transposant. Per exemple
basefils=orth(A')';
retorna una matriu basefils que te per files una base del subespai generat per les files de A (recordo que els ' son transposicions de matrius).
Una altra instruccio que us pot convenir es la null: per exemple
base=null(A);
retorna una matriu base que te per columnes una base del nucli de A (es a dir, una base del subespai amb equacions (A|0)).
Aquestes instruccions surten a la 'Breu introduccio' a l'Octave, a les notes de la practica que vam fer, i als exemples de funcions per passar de generadors a equacions, a base, o de equacions a base per un subespai. Tot aixo esta a la pagina
http://www-ma1.upc.es/~amoros/octave.html
Si encara insistiu en eliminar files o columnes: heu de tallar la matriu en dos trossos i despres enganxar-los de nou com vam fer a la practica (notes a la web op. cit.). Per exemple, per treure la columna j d'una matriu m x n A fariem:
B(:,1:j-1)=A(:,1:j-1);
B(:,j:n-1)=A(:,(j+1):n);
i ja tenim la matriu menys la columna que no voliem.
(pregunta general sobre splines cúbics del 5/1/2004):
Resposta (5/1/2004):
Els splines cubics no cal explicar-los a teoria, son com un problema de la llista! (un problema de splines cubics ha estat a punt de caure a l'examen algun any...). La formacio ja la teniu:
Buscar un polinomi de grau 3 p(x)=a_3 x^3+ a_2 x^2+ a_1 x +a_0 tal que
p(x_1)=y_1
p(x_2)=y_2
p(x_3)=y_3
es resoldre un sistema lineal amb 4 incognites els coeficients a_3,a_2,a_1, a_0, i 3 equacions donades per aquestes igualtats p(x_i)=y_i per i=1,2,3.
Aquest problema d'interpolacio per polinomis ja va sortir a classe de problemes, quan al tema de determinants vareu parlar del determinant de Van der Monde, que es el que apareix en aquests problemes.
El sistema te rang 3 i quatre incognites, aixo vol dir que podeu imposar una 4a condicio (equacio) per a tenir un sistema compatible determinat. Aquesta condicio extra en el cas dels splines cubics es fixar el pendent de la recta tangent en x_1. Com ja heu vist a Calcul 1 (i abans) el pendent de la recta tangent es la derivada, es a dir que vosaltres teniu una pendent inicial m_1 fixada i heu de demanar
p'(x_1)=m_1
Aquesta es la quarta equacio, i amb Octave podeu resoldre en un moment el sistema 4x4.
Feu aixo amb els punts 1,2,3 de la llista, amb pendent inicial m_1=0. Aixi calculeu el primer polinomi, aleshores li mireu la pendent en el punt (x_3,y_3) (es m_3=p'(x_3)), aquesta es la pendent inicial per trobar el polinomi que passi pels punts (x_3,y_3),(x_4,y_4), (x_5,y_5). Repetiu aquest proces fins que s'acabin els punts (cada cop agafeu el darrer punt i dos de nous), i tindreu una colla de polinomis de manera que al fer les grafiques de cada polinomi en els seus 3 punts aquestes grafiques enganxaran de manera llisa, sense punxes (C^1 que diuen a calcul 1).
Per a no fer-ho tot de cop: primer agafeu tres punts (x_1,y_1),(x_2,y_2),(x_3,y_3) i busqueu directament amb octave el polinomi de grau 3 y=p(x) que passa per ells i te pendent 0 en x_1. Despres repetiu el proces imposant una pendent inicial m_1 que no sigui necessariament zero, i quan tingueu aixo nomes us faltara posar-ho dins d'un bucle for - end que vagi calculant el polinomi per cada terna de punts.
Actualitzada el 2004/1/9.
Tornar a l'índex
de la web de Jaume Amorós.