PROGRAMACION DE METODOS ITERATIVOS:
//METODO DE LAS POTENCIAS
A=[-2 1 3;-3 1 1;-4 2 2] //se debe modificar en cada ejercicio
X=ones([size(A),1]') //se debe de modificar para cada ejercicio
for n=1:5 //se modifica segun las iteracciones que requiera el ejerc
anterior=X;
X=A*anterior;
if(min(abs(X)))<>0;X=X./min(abs(X));end
end
autoval=X./anterior
autovaldom=X./max(abs(autoval))
autovec=X./min(abs(X'))
//METODO DE FRANCIS
A=[2 -1 0;1 2 3;-1 -2 2] // se modifica segn ejer
U=eye(A)
for n=1:10 // ''''''''''''''''''''''
[Q R]=qr(A);
U=U*R;
A=Q*U;
end
A
autovalores=(diag(A))
U
//METODO JACOBI
function [x,res,rad, BJ,CJ]=jacobi(A,b,n,x0)
// Metodo iterativo de Jacobi en forma matricial
// A=matriz del sistema
// b= vector de terminos independientes
// x0=vector inicial (columna)
// n= numero de iteraciones
// x= solucion tras las iteraciones
// res= vector residuo tras las iteraciones =Ax-b
// rad= radio espectral de la matriz de Jacobi
// BJ,CJ= matriz y vector del m'etodo iterativo
if argn(2)<2,disp('error, Jacobi no dispone de suficientes argumentos'),end
if argn(2)==2, n=100; x0=zeros(size(b));end
if argn(2)==3, x0=zeros(size(b)),end
D=diag(diag(A));L=tril(A)-D;R=triu(A)-D;ID=inv(D);
BJ=-ID*(L+R);CJ=ID*b;
x=x0;
for k=1:n; // Iteraciones
x=BJ*x+CJ;
end
if argn(1)>1,res=A*x-b;end
if argn(1)>2,rad=max(abs(spec(BJ)));end
endfunction
//METODO GAUSS-SEIDEL
function [x,res,rad,BGS,CGS]=gseidel(A,b,n,x0)
// Metodo iterativo de Gauss-Seidel en forma matricial
// A=matriz del sistema
// b= vector de terminos independientes
// x0=vector inicial (columna)
// n= numero de iteraciones
// x= solucion tras las iteraciones
// res= vector residuo tras las iteraciones =Ax-b
// rad= radio espectral de la matriz de Gauss-Seidel
// BGS,CGS= matriz y vector del m'etodo iterativo
if argn(2)<2,disp('error, Gauss-Seidel no dispone de suficientes argumentos'),end
if argn(2)==2, n=100; x0=zeros(size(b));end
if argn(2)==3, x0=zeros(size(b)),end
D=diag(diag(A));L=tril(A)-D;R=triu(A)-D;ID=inv(L+D);
BGS=-ID*R;CGS=ID*b;
x=x0;
for k=1:n; // Iteraciones
x=BGS*x+CGS;
end
if argn(1)>1,res=A*x-b;end
if argn(1)>2,rad=max(abs(spec(BGS)));end
endfunction
function [x,res,rad,BSOR,CSOR]=sor(A,b,w,n,x0)
// Metodo iterativo de Gauss-Seidel en forma matricial
// A=matriz del sistema
// b= vector de terminos independientes
// w=parametro de sobrerelajación
// x0=vector inicial (columna)
// n= numero de iteraciones
// x= solucion tras las iteraciones
// res= vector residuo tras las iteraciones =Ax-b
// rad= radio espectral de la matriz de Gauss-Seidel
// BSOR,CSOR= matriz y vector del m'etodo iterativo
if argn(2)<2,disp('error, SOR no dispone de suficientes argumentos'),end
if argn(2)==2, w=1.2;n=100;x0=zeros(size(b));end
if argn(2)==3, n=100; x0=zeros(size(b));end
if argn(2)==4, x0=zeros(size(b)),end
D=diag(diag(A));L=tril(A)-D;R=triu(A)-D;ID=inv(D+w*L);
BSOR=ID*((1-w)*D-w*R);CSOR=ID*w*b;
x=x0;
for k=1:n; // Iteraciones
x=BSOR*x+CSOR;
end
if argn(1)>1,res=A*x-b;end
if argn(1)>2,rad=max(abs(spec(BSOR)));end
endfunction
//RELACION DE EJERCICIOS SCILAB
//Ejercicio 1: dado el sistema Ax=b
//Datos
A=[1 10 -10;10 1 -10;10 -10 1]
b=[1 2 2]'
//Apartado A: resolver con Scilab
X=A\b
//Apartado B: hallar la norma 1, norma 2 y norma infinita del vector residuo en los 3 casos.
residuo=(A*X)-b
// residuo =
//
// 0.
// 0.
// - 6.661D-16
norm(residuo,1)
//ans =
//
// 6.661D-16
norm(residuo)
//ans =
//
// 6.661D-16
norm(residuo,'inf')
//ans =
//
// 6.661D-16
//Apartado C:calculla la escalonada de A
rref(A)
//ans =
//
// 1. 0. 0.
// 0. 1. 0.
// 0. 0. 1.
//Apartado D: calcular el determinate de A
det(A)
//ans =
//
// - 99.
//Apartado E: calcula la inversa de A
inv(A)
//ans =
//
// 1. - 0.9090909 0.9090909
// 1.1111111 - 1.020202 0.9090909
// 1.1111111 - 1.1111111 1.
//Ejercicio 2: dado el sistema Ax=b
//Datos
A=[1 4 2 -1;-4 2 2 1;2 4 3 4;-1 1 4 1]
//A =
//
// 1. 4. 2. - 1.
// - 4. 2. 2. 1.
// 2. 4. 3. 4.
// - 1. 1. 4. 1.
b=[0 1 2 3]'
//b =
//
// 0.
// 1.
// 2.
// 3.
//Apartado A: calcula norma 1 de l matriz A y del vector b
norm(A,1)
//ans =
//
// 11.
norm(b,1)
//ans =
//
// 6.
//Apartadi B:factoriza la matriz A por metodo QR y resolver a continuacion
[Q R]=qr(A)
//DESCOMPOSICION EN QR
//Matriz ortogonal
//R =
//
// - 4.6904158 - 0.6396021 0.8528029 - 0.4264014
// 0. - 6.049042 - 4.7190042 - 2.4346455
// 0. 0. - 3.1628668 - 1.5411759
// 0. 0. 0. - 3.242755
////Matriz diagonal superior
// Q =
//
// - 0.2132007 - 0.6387187 0.2631468 0.6908963
// 0.8528029 - 0.4208029 0.2254421 - 0.2117263
// - 0.4264014 - 0.6161757 - 0.1441416 - 0.6463223
// 0.2132007 - 0.1878584 - 0.9269050 0.2451567
//Resuelvo QR
x=R\Q'*b
//x =
//
// 0.0309278
// - 0.3402062
// 0.7835052
// 0.2371134
//Apartado C: calcular por Jacobi dando 2 iteracciones y calcular la norma 1 del residuo
[x,res,rad, BJ,CJ]=jacobi(A,b,2,zeros(4,1))
// CJ =
//
// []
// BJ =
//
// 0. - 4. - 2. 1.
// 2. 0. - 1. - 0.5
// - 0.6666667 - 1.3333333 0. - 1.3333333
// 1. - 1. - 4. 0.
// rad =
//
// 2.8698637
// res =
//
// 0.
// 0.
// 0.
// 0.
// x =
//
// 0.
// 0.
// 0.
// 0.
norm(res,1)
//ans =
//
// 0
//Ejercicio 3: dado el sistema Ax=b
//Datos
A=[0 1 2 3; 1 0 -3 0;0 2 3 4;2 3 0 0]
// A =
//
// 0. 1. 2. 3.
// 1. 0. - 3. 0.
// 0. 2. 3. 4.
// 2. 3. 0. 0.
b=[5 -1 4 8.5]'
//b =
//
// 5.
// - 1.
// 4.
// 8.5
//Apartado A:Hallar la matriz inversa por el metodo de gaussjordan y resolver por gaussjordan.
//Inversa
inv(A)
// ans =
//
// 4. - 0.3333333 - 3. 0.6666667
// - 2.6666667 0.2222222 2. - 0.1111111
// 1.3333333 - 0.4444444 - 1. 0.2222222
// 0.3333333 0.2222222 0. - 0.1111111
//Resolver por gauss-jordan
[x,e]=gaussjordan(A,b)
// !--error 27
//División por cero...
//at line 18 of function gaussjordan called by :
//gaussjordan(A,b)
rank(A)
// ans =
//
// 4.
rank([A b])
// ans =
//
// 4.
// compruevo los rangos para ver si el sistema es sobreterminado, si es sobredeterminado Scilab te devuelve error porque es un sistema incompatible y no puede resolverse por gauss-jordan o gauss.
//El comando gaussjordan debe presentar algun error en su programacion debido a que este tipo de sistema deberia resolverlo ya que no esta sobredeterminado, es decir es compatible determinado.
//EJERCICIO 4: dado un sistema Ax=b
//Datos
A=[2 -1 0 0;1 2 -1 0;0 1 2 -1;0 0 1 2]
b=[2 0 1 3]'
//Apartado A:Realizar dos iteracciones por jacobi partiendo del vector nulo
[x,res,rad, BJ,CJ]=jacobi(A,b,2,zeros(4,1))
// CJ =
//
// 1.
// 0.
// 0.5
// 1.5
// BJ =
//
// 0. 0.5 0. 0.
// - 0.5 0. 0.5 0.
// 0. - 0.5 0. 0.5
// 0. 0. - 0.5 0.
// rad =
//
// 0.8090170
// res =
//
// 0.25
// - 0.75
// 0.
// 0.75
// x =
//
// 1.
// - 0.25
// 1.25
// 1.25
//Apartado B: estudiar el error cometido por jacobi
norm(res)
// ans =
//
// 1.0897247
cond(BJ)
// ans =
//
// 2.618034
//tras las dos iteracciones es muy buena la solucion y con seguridad el sistema converge por que el radio espectral es menor de 1 y sieempre que es rad < 1 -->converge. A parte tambien le hago el modulo al vector residuo y me da casi 1 por tanto es otro factor que me indica que el sitema seguramente converge.tambien veo si ha sido bien condicionada y asi es por que el resultado es 1.Esto son propiedades de los metodos iterativos para evaluar los errores cometidos y asi ver si la solucion dada por el metodo es buena o no.
//Apartado C: estudiar la convergencia del metodo
[x,res,rad]=jacobi(A,b,10000000,zeros(4,1))
// rad =
//
// 0.8090170
// res =
//
// 0.
// 0.
// 0.
// 0.
// x =
//
// 1.
// 0.
// 1.
// 1.
//como podemos observar cuantas mas iteracciones el error es menor dando in vector residuo nulo y cuyo radio espectral < 1 ,esto implica que el metodo coverge seguro.
//EJERCICIO 5: resolver por el metodo QR y por Jacobi(100 iteracciones partiendo del vector nulo)de los siguientes sistemas y calcular sus residuos en norma 2
//SISTEMA A
//2x+3y-z=2 x+y-2y=-4 2x+y+z=6
A=[2 3 -1;1 1 -2;2 1 1]
// A =
//
// 2. 3. - 1.
// 1. 1. - 2.
// 2. 1. 1.
b1=[2 -4 6]'
//b1 =
//
// 2.
// - 4.
// 6.
//SISTEMA B
//2x+3y=5 x+y=2 x-y=1
B=[2 3;1 1;1 -1]
// B =
//
// 2. 3.
// 1. 1.
// 1. - 1.
b2=[5 2 1]'
// b2 =
//
// 5.
// 2.
// 1.
//Metodo QR
//Sistema A
[Q R]=qr(A)
// R =
//
// - 3. - 3. 0.6666667
// 0. 1.4142136 - 1.4142136
// 0. 0. 1.8856181
// Q =
//
// - 0.6666667 0.7071068 0.2357023
// - 0.3333333 - 1.388D-16 - 0.9428090
// - 0.6666667 - 0.7071068 0.2357023
//calculo el residuo de QR
X=R\(Q'*b1)
// X =
//
// 1.
// 1.
// 3
norm(X)
//ans =
//
// 3.3166248
//Sistema B
[Q R]=qr(B)
// R =
//
// - 2.4494897 - 2.4494897
// 0. 2.236068
// 0. 0.
// Q =
//
// - 0.8164966 0.4472136 0.3651484
// - 0.4082483 - 2.776D-17 - 0.9128709
// - 0.4082483 - 0.8944272 0.1825742
//calculo el residuo del QR
X=R\(Q'*b2)
// X =
//
// 1.5666667
// 0.6
norm(X)
// ans =
//
// 1.6776306
//Metodo Jacobi
//Sistema A
[x,res,rad, BJ,CJ]=jacobi(A,b1,100,zeros(3,1))
// CJ =
//
// 1.
// - 4.
// 6.
// BJ =
//
// 0. - 1.5 0.5
// - 1. 0. 2.
// - 2. - 1. 0.
// rad =
//
// 2.0152901
// res =
//
// 10^30 *
//
// 3.4611361
// - 9.8875619
// 16.62628
// x =
//
// 10^29 *
//
// 41.497177
// 8.7213618
// 74.547079
//calculo la norma 2 del residuo de jacobi
norm(res)
//ans =
//
// 1.965D+31
//Sistema B
[x,res,rad, BJ,CJ]=jacobi(B,b2,100,zeros(3,1))
//tengo que construir una matriz que sea cuadrada para que sea util este metodo
B=[2 3 0;1 1 0;1 -1 0;]
// B =
//
// 2. 3. 0.
// 1. 1. 0.
// 1. - 1. 0.
//construyo vector nulo para que noo afecte al sistema pero si pueda operar el metodo al ser ya cuadrada dicha matriz
[x,res,rad, BJ,CJ]=jacobi(B,b2,100,zeros(3,1))
// !--error 19
//El problema es singular.
//at line 15 of function jacobi called by :
//[x,res,rad, BJ,CJ]=jacobi(B,b2,100,zeros(3,1))
//tampoco sirve el metodo porque el metodo no esta bien programado y cuando algun elemento de la traza es 0 no puede calcular nada.
//Ejercicio 6:
//Datos
A=[2 2 1 0;1 2 2 0;0 1 2 2;0 1 1 4]
b=[2 1 0 -1]'
//Apartado A: estudiar la compatibilidad del sistema y explicar que tipo de sistema tenemos
det(A)
//ans =
//
// 6.
rank(A)
//ans =
//
// 4.
Ampliada=[A b]
//Ampliada =
//
// 2. 2. 1. 0. 2.
// 1. 2. 2. 0. 1.
// 0. 1. 2. 2. 0.
// 0. 1. 1. 4. - 1.
rank(ampliada)
//ans =
//
// 3.
//el rango(A)>rango(ampliada) ---> el sistema es compatible indeterminado segun el teorema de Rouche-Frobenius.
//Apartado B: resolver usando el comando de scilab y calcular la norma 2 del error cometido
X=A\b
// X =
//
// 1.6666667
// - 1.
// 0.6666667
// - 0.1666667
error=(A*X)-b
//error =
//
// 0.
// 0.
// - 5.551D-17
// 0.
norm(error,2)
// ans =
//
// 5.551D-17
//Apartado C: calcular por gauss y calcular la norma 2 del error cometido
x=gauss(A,b)
// x =
//
// 1.6666667
// - 1.
// 0.6666667
// - 0.1666667
x2=(A*x)-b
//x2 =
//
// 0.
// 0.
// - 5.551D-17
// 0.
norm(x2)
//ans =
//
// 5.551D-17
//Ejercicio 7: dado el sistema Ax=b
//Datos
A=[4 3 1;-1 4 1;-1 1 4]
// A =
//
// 4. 3. 1.
// - 1. 4. 1.
// - 1. 1. 4
b=[1 1 0]'
//b =
//
// 1.
// 1.
// 0.
x0=ones(3,1)
//x0 =
//
// 1.
// 1.
// 1
//Apartado A: ¿Que tipo de sistema tenemos?
[x,res,rad, BJ,CJ]=jacobi(A,b,10,x0)
// CJ =
//
// []
// BJ =
//
// 0. - 0.75 - 0.25
// 0.25 0. - 0.25
// 0.25 - 0.25 0.
// rad =
//
// 0.5
// res =
//
// 0.0098724
// 0.0014458
// 0.0014458
// x =
//
// 0.0018158
// 0.0006523
// 0.0006523
//el sistema converge por tanto es compatible determinado
//Apartado B:¿Es la matriz diagonal dominante por filas o por columnas?
//Por Filas
norm(A(1,[2 3]))
// ans =
//
// 3.1622777
norm(A(2,[1 3]))
//ans =
//
// 1.4142136
norm(A(3,[1 2]))
// ans =
//
// 1.4142136
//Es dominate por filas
//Por Columnas
norm(A([2 3],1))
// ans =
//
// 1.4142136
norm(A([1 3],2))
// ans =
//
// 3.1622777
norm(A([1 2],3))
// ans =
//
// 1.4142136
//Es dominate por columnas.
//Para que la matriz sea extrictamente dominate, tiene que ser dominante por filas y por columnas.
//Es dominate por filas cuando para todas las filas el elemento de la diagonal en valor absoluto es mayor que la suma del resto de elementos en valor absoluto de esa fila
//Es dominate por columnas cuando para todas las columnas el elemento de la diagonal en valor absoluto es mayor que la suma del resto de elementos en valor absoluto de esa columna.
//TEOREMA DE HADAMARD:Podemos enunciar a partir de esta definición el siguiente teorema de convergencia aplicable a los procesos iterativos de Jacobi y Gauss Seidel:
// Sea A una matriz cuadrada, si A es diagonal dominante, los métodos iterativos de Jacobi y Gauss-Seidel convergen a la solución del sistema de ecuaciones Ax=b.
//Apartado C: resolver por gauss y calcular la norma 2 del error cometido.
X=gauss(A,b)
//X =
//
// 0.0555556
// 0.2777778
// - 0.0555556
error=(A*X)-b
// error =
//
// 0.
// 0.
// 2.776D-17
norm(error)
// ans =
//
// 2.776D-17
//Apartado D:resolver con jacobi dando dos iteracciones y estudia la convergencia de jacobi y el error cometido.
[x,res,rad, BJ,CJ]=jacobi(A,b,2,zeros(3,1))
//CJ =
//
// []
// BJ =
//
// 0. - 0.75 - 0.25
// 0.25 0. - 0.25
// 0.25 - 0.25 0.
// rad =
//
// 0.5
// res =
//
// 0.
// 0.
// 0.
// x =
//
// 0.
// 0.
// 0.
//Estudio de la convergencia
//el sistema converge segun el metodo de jacobi porque el radio expectral < 1.
//Estudio el error cometido
//el error cometido es 0 como bien muestra el vector residuo que obtenemos tras iterar por Jacobi.
//Ejercicio 8: dado el sistema Ax=b resolver por gauss,QR,Jacobi(50 iteracciones partiendo del vector nulo).Calcula el error utilizando la norma infinita en los tres casos.
//Datos
A=[1 -1 3;-1 10 3;3 3 14]
// A =
//
// 1. - 1. 3.
// - 1. 10. 3.
// 3. 3. 14.
b=[-1 2 0]'
//b =
//
// - 1.
// 2.
// 0.
X=gauss(A,b)
// X =
//
// - 9.4444444
// - 1.4444444
// 2.3333333
x=(A*X)-b
//x =
//
// 0.
// 0.
// 0.
norm(x,'inf')
//ans =
//
// 0.
[Q R]=qr(A)
// R =
//
// - 3.3166248 0.6030227 - 12.663476
// 0. - 10.470738 - 7.3190983
// 0. 0. 0.2591605
// Q =
//
// - 0.3015113 0.0781398 - 0.9502553
// 0.3015113 - 0.9376781 - 0.1727737
// - 0.9045340 - 0.3386060 0.2591605
X=R\(Q'*b)
// X =
//
// - 9.4444444
// - 1.4444444
// 2.3333333
norm(X,'inf')
// ans =
//
// 9.4444444
[x,res,rad, BJ,CJ]=jacobi(A,b,50,zeros(3,1))
// CJ =
//
// []
// BJ =
//
// 0. 1. - 3.
// 0.1 0. - 0.3
// - 0.2142857 - 0.2142857 0.
// rad =
//
// 0.9694175
// res =
//
// 0.
// 0.
// 0.
// x =
//
// 0.
// 0.
// 0.
norm(res,'inf')
//ans =
//
// 0.
//Ejercicio 9: dada la matriz A y los vectores b1,b2,b3
//Datos
A=[1 1 2;2 5 -2;3 6 1]
//A =
//
// 1. 1. 2.
// 2. 5. - 2.
// 3. 6. 1.
b1=[1 2 0]'
//b1 =
//
// 1.
// 2.
// 0.
b2=[0 1 3]'
//b2 =
//
// 0.
// 1.
// 3.
b3=[1 2 1]'
// b3 =
//
// 1.
// 2.
// 1
//Apartado A:Hallar la matriz inversa(A) y resolver conjuntamente los sistemas Ax=b1; Ax=b2; Ax=b3 por un metodo directo
inv(A)
// ans =
//
// 5.6666667 3.6666667 - 4.
// - 2.6666667 - 1.6666667 2.
// - 1. - 1. 1.
X=gauss(A,[b1 b2 b3])
// X =
//
// 13. - 8.3333333 9.
// - 6. 4.3333333 - 4.
// - 3. 2. - 2.
//Apartado B:calcula la escalonada de A
rref(A)
// ans =
//
// 1. 0. 0.
// 0. 1. 0.
// 0. 0. 1.
//Ejercicio 10: dado el sistema
A=[-3 1 1 1;1 -3 1 1;1 1 -3 1;1 1 1 -3]
//A =
//
// - 3. 1. 1. 1.
// 1. - 3. 1. 1.
// 1. 1. - 3. 1.
// 1. 1. 1. - 3.
b=[0 8 -8 1]'
// b =
//
// 0.
// 8.
// - 8.
// 1.
//Apartado A: estudiar la compatibilidad del sistema
rank(A)
//ans =
//
// 3.
ampliada=rank([A b])
// ans =
//
// 4.
// rango(A)< rango(ampliada) --->Segun el teorema de ROUCHE-FROBENIUS el sistema es compatible indeterminado.
//Apartado B:
X=gauss(A,b)
// !--error 27
//División por cero...
//at line 7 of function sustreg called by :
//at line 29 of function gauss called by :
//X=gauss(A,b)
//el problema que te da es debido a que hay una fila que es linealmente dependiente de alguna de las otras o bien es combinacion lineal entre las otras o entre algunas de ellas.
[Q R]=qr(A)
// R =
//
//
// column 1 to 4
//
// 3.4641016 - 1.1547005 - 1.1547005 - 1.1547005
// 0. 3.2659863 - 1.6329932 - 1.6329932
// 0. 0. 2.8284271 - 2.8284271
// 0. 0. 0. 1.110D-15
//
// column 5
//
// 0.2886751
// - 9.3897107
// 6.363961
// 0.5
// Q =
//
// - 0.8660254 5.179D-17 5.179D-17 0.5
// 0.2886751 - 0.8164966 4.241D-17 0.5
// 0.2886751 0.4082483 - 0.7071068 0.5
// 0.2886751 0.4082483 0.7071068 0.5
[Q R]=qr([A b])
// R =
//
//
// column 1 to 4
//
// 3.4641016 - 1.1547005 - 1.1547005 - 1.1547005
// 0. 3.2659863 - 1.6329932 - 1.6329932
// 0. 0. 2.8284271 - 2.8284271
// 0. 0. 0. 1.110D-15
//
// column 5
//
// 0.2886751
// - 9.3897107
// 6.363961
// 0.5
// Q =
//
// - 0.8660254 5.179D-17 5.179D-17 0.5
// 0.2886751 - 0.8164966 4.241D-17 0.5
// 0.2886751 0.4082483 - 0.7071068 0.5
// 0.2886751 0.4082483 0.7071068 0.5
//el metodo es totalmente valido para calcular este sistema.
//Ejercicio 11: calcular el autovalor dominante(y su autovector unitario en norma infinita)por el metodo de la potencia, realizando 5 iteracciones partiendo del vector (1 1 1)'
//METODO DE LAS POTENCIAS
A=[-2 1 3;-3 1 1;-4 2 2]
X=ones([size(A),1])//se debe de modificar para cada ejercicio
for n=1:5
anterior=X;
X=A*anterior;
if(min(abs(X)))<>0;X=X./min(abs(X));end
end
autoval=X./anterior
autovaldom=X./max(abs(autoval))
autovec=X./min(abs(X'))
//SOLUCION:
//-3-> A=[-2 1 3;-3 1 1;-4 2 2]
// A =
//
// - 2. 1. 3.
// - 3. 1. 1.
// - 4. 2. 2.
//
//-3-> X=ones([size(A),1])
// X =
//
// 1. 1. 1.
//
// for n=1:5
// anterior=X;
// X=A*anterior;
// if(min(abs(X)))<>0;X=X./min(abs(X));end
// end
// !--error 10
// Multiplicación inconsistente.
//
//autoval=X./anterior
// autoval =
//
// 1. 1. 1.
//
//autovaldom=X./max(abs(autoval))
// autovaldom =
//
// 1. 1. 1.
//
// autovec=X./min(abs(X'))
// autovec =
//
// 1. 1. 1.
//Ejercicio 12: calcular autovalores y autovectores usando la orden de Scilab y el metodo de Francis partiendo del ventor (1 1 1)'
//Datos
A=[2 -1 0;1 2 3;-1 -2 2]
// A =
//
// 2. - 1. 0.
// 1. 2. 3.
// - 1. - 2. 2.
//Por el comando Scilab
spec(A)
//ans =
//
// 2.4181283
// 1.7909359 + 2.6704163i
// 1.7909359 - 2.6704163i
//Por el metodo de Francis
//METODO DE FRANCIS
A=[2 -1 0;1 2 3;-1 -2 2]
U=eye(A)
for n=1:10
[Q R]=qr(A);
U=U*R;
A=Q*U;
end
A
autovalores=(diag(A))
U
//SOLUCION:
//A=[2 -1 0;1 2 3;-1 -2 2]
// A =
//
// 2. - 1. 0.
// 1. 2. 3.
// - 1. - 2. 2.
//
//-3->U=eye(A)
// U =
//
// 1. 0. 0.
// 0. 1. 0.
// 0. 0. 1.
//
//for n=1:10
// [Q R]=qr(A);
// U=U*R;
// A=Q*U;
//end
//
//A
// A =
//
// 10^239 *
//
// - 1.314D-40 - 0.0005083 - 4.007D+39
// - 6.571D-41 - 0.0007188 5.036D+41
// 6.571D-41 0.0007188 4.073D+41
//
//autovalores=(diag(A))
// autovalores =
//
// 10^239 *
//
// - 1.314D-40
// - 0.0007188
// 4.073D+41
//
//U
// U =
//
// 10^239 *
//
// 1.610D-40 0.0010019 - 3.604D+40
// 0. 0.0005365 - 5.790D+40
// 0. 0. 6.441D+41
//control de practicas:ejercicio de repaso
//datos
A=[-0.001 0.2 0.2;0.333 8 -1;-1 0.5 -2]
b=[5 -1 2]'
//Apartado A: estudiar la compatibilidad del sistema y explicar que tipo de sistema es.
rank(A)
// ans =
//
// 3.
ampliada=[A b]
// ampliada =
//
// - 0.001 0.2 0.2 5.
// 0.333 8. - 1. - 1.
// - 1. 0.5 - 2. 2.
det(A)
rank(ampliada)
// ans =
//
// 3.
//Para una primera apreciacion de que tipo de sistema es si compatible o incompatible hago el determinate(A) y veo que es distinto de cero por tanto ya se que va a ser compatible determinado y atendiento al teorema de Rouche-frobenius como el rango(A)=rango(A|b)=3 y tiene el mismo numero de incognitas que sale al hacer el rango de ambas por tanto es un SISTEMA COMPATIBLE DETERMINADO.
//Aparatado B: Calcular el modulo de A y el modulo de b
norm(A)
//8.1037701
norm(b)
// 5.4772256
//Apartado C:
//metodo QR
[Q R]=qr(A)
// R =
No hay comentarios:
Publicar un comentario