Trabajo Final Metodos Numericos (Crout-Jacobi)

Universidad de Cuenca
Facultad de Ingeniería
Escuela Civil
Materia: Métodos Numéricos
Profesor: Ing. Fernando Zalamea
Tema: Factorización Código en Fortran y Python.
• METODO CROUT OPTIMIZADO.

• METODO ITERATIVO DE JACOBI
• METODO ITERATIVO DE GAUSS-SEIDEL
Ciclo: Septiembre-Febrero 2014
Cuenca, 24 de enero del 2014.

OBJETIVOS:
➢ Realizar el código en Fortran de la factorización de una matriz A en dos matrices L

(triangular inferior) y U (triangular superior).
➢ Utilizar el método de Crout para encontrar las matrices L y U.
➢ Resolver el sistema mediante sustitución.
➢ Realizar el código de Jacobi y Gauss-Seidel en distintos lenguajes de
programación.
DESARROLLO:
Ejercicio 1:
Realice un programa en Fortran y Python , que obtenga las
soluciones de un sistema [A]*[X]=[B] ,considerando a la matriz
[A]nxn que es dE tipo banda pentadimensional (5 diagonales de
banda),mediante el método de Crout optimizado para este tipo de
matriz.
a) PROGRAMA EN FORTRAN.
program crout_optimizado
integer::n,i,j
real::a(100,100),u(100,100),l(100,100),b(100),x(100),z(100)
write(*,*)"metodo de crout para una matriz banda de ancho 5"
write(*,*)"ingrese la dimension de la matriz"
read(*,*)n
! no permite el ingreso de matrices con dimension menor a 6

do while (n.lt.6)
write(*,*)"ingrese nuevamente la dimension de la matriz"
read(*,*)n
end do
!creacion de la matriz de ceros

do i=1,n
do j=1,n
a(i,j)=0.0 ! llenado de la matriz a con ceros
u(i,j)=0.0 ! llenado de la matriz u con ceros
l(i,j)=0.0 ! llenado de la matriz l con ceros
end do
end do
!entrada de datos
do i=1,n
if (i.le.3) then
do j=1,(2+i)
write(*,*)"ingrese el valor de la matriz "
read(*,*)a(i,j)
end do
else
if ((2+i).le.n) then
do j=i-2,2+i
read(*,*)a(i,j)
end do
else
do j=i-2,n
read(*,*)a(i,j)
end do
end if
end if
end do
!ingreso de los terminos independientes

do i=1,n
write(*,*)"ingrese el termino independiente del sistema"
read (*,*)b(i)
end do
do i=1,n
if (i.le.3) then
sum1=0.0
u(i,i)=1.0
do k=1,i
sum1=sum1+l(i,k)*u(k,i)
end do
l(i,i)=a(i,i)-sum1
else
sum1=0.0
u(i,i)=1.0
do k=i-2,i
sum1=sum1+l(i,k)*u(k,i)
end do
l(i,i)=a(i,i)-sum1
end if
if (abs (l(i,i)).lt.0.0001)then
write(*,*)"la factorizacion no se puede crear"
stop
end if
if (i+3.ge.n)then
do j=i+1,i+3
if (abs(i-j).eq.2)then
u(i,j)=a(i,j)/l(i,i)
l(j,i)=a(j,i)
else
k=i-1
sum1=0.0
sum2=0.0
sum1=sum1+l(i,k)*u(k,j)
u(i,j)=(a(i,j)-sum1)/l(i,i)
sum2=sum2+l(j,k)*u(k,i)
l(j,i)=(a(j,i)-sum2)
end if
end do
else
do j=i+1,n !para no tomar en cuenta los ceros y disminuir
operaciones
if (abs(i-j).eq.2)then
u(i,j)=a(i,j)/l(i,i)
l(j,i)=a(j,i)
else
k=i-1
sum1=0.0
sum2=0.0
sum1=sum1+l(i,k)*u(k,j)
u(i,j)=(a(i,j)-sum1)/l(i,i)
sum2=sum2+l(j,k)*u(k,i)
l(j,i)=(a(j,i)-sum2)
end if
end do
end if
end do
! para encontrar los valores de z y los valores solucion x
do i=1,n
x(i)=0.0
z(i)=0.0
end do
z(1)=b(1)/l(1,1)
z(2)=(b(2)-l(2,1)*z(1))/l(2,2)
do i=3,n
suma=0.0
do j=i-2,i
suma=suma+l(i,j)*z(j)
end do
z(i)=(b(i)-suma)/l(i,i)
end do
x(n-1)=z(n-1)
x(n-2)=z(n-2)-(u(n-2,n-1)*x(n-1))
do i=n-3,-1,-1
sum=0.0
do j=i+1,i+3
sum=sum+u(i,j)*x(j)
end do
x(i)=z(i)-sum
end do
write(*,*)"la matriz a factorizar es:"
do i=1,n
write(*,*)(a(i,j),j=1,n)
end do
write(*,*)"la matriz l es:"
do i=1,n
write(*,*)(l(i,j),j=1,n)
end do
write(*,*)"la matriz u es:"
do i=1,n
write(*,*)(u(i,j),j=1,n)
end do
write(*,*)"el valor de la incognita es: "
do i=1,n
write(*,*)x(i)
end do
read(*,*)
end program
b) PROGRAMA EN PYTHON.
from numpy import*
import numpy
import sys
#la funcion realiza el llenado de la matriz y los terminos independientes

def matriz(n):
matriz=zeros((n,n),float)
for i in range(0,n):
if i<3:
for j in range(0,3+i):
matriz[i][j]=float(input("*Ingrese el valor de la matriz
en la fila "+ str(i+1)+" columna "+ str(j+1)+":"))
else:
if (3+i)<=n:
for j in range(i-2,3+i):
matriz[i][j]=float(input("*Ingrese el valor de la
matriz en la fila "+ str(i+1)+" columna "+ str(j+1)+":"))
else:
for j in range(i-2,n):
matriz[i][j]=float(input("*Ingrese el valor de la
matriz en la fila "+ str(i+1)+" columna "+ str(j+1)+":"))
ind=[]
for i in range(n):
b=float(input("*Ingrese el termino independiente "+str(i+1)+"del
sistema [A][x]=[b]:"))
ind.append(b)
return (matriz,ind)
#esta funcion nos devuelve las matrices l y u para utilizarlas en la

solucion del sistema
def Crout(a,n,b):
u=zeros((n,n),float)
l=zeros((n,n),float)
for i in range(n):
if i<3:
sum1=0.0
u[i][i]=1.0
for k in range(i):
sum1=sum1+l[i][k]*u[k][i]
l[i][i]=a[i][i]-sum1
else:
sum1=0.0
u[i][i]=1.0
for k in range(i-2,i):
sum1=sum1+l[i][k]*u[k][i]
l[i][i]=a[i][i]-sum1
if abs (l[i][i])< 0.00001:
print"**La factorizacion no se puede crear**"
sys.exit(" ")
if i+3<=n:
for j in range(i+1,i+3):
if abs(i-j)==2:
u[i][j]=(a[i][j])/l[i][i]
l[j][i]=(a[j][i])
else:
k=i-1
sum1=0.0
sum2=0.0
sum1=sum1+(l[i][k]*u[k][j])
u[i][j]=((a[i][j]-sum1)/l[i][i])
sum2=sum2+(l[j][k]*u[k][i])
l[j][i]=(a[j][i]-sum2)
else:
for j in range(i+1,n): #para no tomar en cuenta los ceros y
no realizar operaciones innecesarias
if abs(i-j)==2:
u[i][j]=(a[i][j])/l[i][i]
l[j][i]=(a[j][i])
else:
k=i-1
sum1=0.0
sum2=0.0
sum1=sum1+(l[i][k]*u[k][j])
u[i][j]=((a[i][j]-sum1)/l[i][i])
sum2=sum2+(l[j][k]*u[k][i])
l[j][i]=(a[j][i]-sum2)
print
print "La matriz A Factorizar es:"
print a
print '*************************'
print 'El termino independiente es:'
print b
print'**************************'
print "La matriz L es: "
print l
print'**************************'
print "La matriz U es: "
print u
return (l,u)
#resolucion del sitema aplicando las dos formulas
#1.[L][Z]=[b]
#2.[U][X]=[z]
def sistema(l,u,b,n):
x=[]
z=[]
for i in range(n):
x.append(0.0)
z.append(0.0)
z[0]=b[0]/l[0][0]
z[1]=(b[1]-l[1][0]*z[0])/l[1][1]
for i in range(2,n):
suma=0.0
for j in range(i-2,i):
suma=suma+l[i][j]*z[j]
z[i]=(b[i]-suma)/l[i][i]
x[n-1]=z[n-1]
x[n-2]=z[n-2]-(u[n-2][n-1]*x[n-1])
for i in range(n-3,-1,-1):
sum=0.0
for j in range (i+1,i+3):
sum=sum+u[i][j]*x[j]
x[i]=z[i]-sum
return x
# llamadas de las funciones y obtencion de datos de las mismas
print "METODO DE CROUT PARA UNA MATRIZ BANDA DE ANCHO 5"

print'Factorizacion y resolucion de un sistema [A][x]=[b]'
n=int(input("Ingrese la dimension de la matriz banda A : "))
while n<6:
n=int(input("Ingrese nuevamente la dimension de la matriz banda A :
"))
valormat=matriz(n)
a=valormat[0]
ter_ind=valormat[1]
lu=Crout(a,n,ter_ind)
ml=lu[0]
mu=lu[1]
resultado=sistema(ml,mu,ter_ind,n)
print '***************'
for i in range(len(resultado)):
print 'El valor de la incognita x'+str(i+1)+' es :'+'
'+str(resultado[i])
Ejercicio 2:
− Realice un programa en Fortran del Metodo Iterativo de Jacobi
para la resolución de sistemas [A][X]=[B]
− Realice un programa en Python del Metodo Iterativo de Gauss-
Seidel para la resolución de sistemas [A][X]=[B]
a)PROGRAMA EN FORTRAN M.JACOBI.

program jacobi
real a(100,100),b(100),x(100),x0(100),y(100)
integer n,i,j,k
real sum1,sum2,sum3,sum4,error,nor,norx
!INGRESO DEL TAMANO DE LA MATRIZ

!INGRESO DE LOS VALORES DE LA MATRIZ
!INGRESO DEL TERMINO INDEPENDIENTE
write(*,*)"ingrese el tamano de la matriz"
read(*,*)n
do i=1,n
do j=1,n
write(*,*)"ingrese los coeficientes de la matriz",i,j
read(*,*)a(i,j)
end do
end do
do i=1,n
b(i)=0.0
x0(i)=0.0
x(i)=0.0
y(i)=0.0
end do
do i=1,n
write(*,*)"ingrese el valor del termino independiente "
read (*,*)b(i)
end do
write(*,*)"La matriz A ingresada es :"
write(*,*)
do i=1,n
write(*,*)(a(i,j),j=1,n)
end do
write(*,*)'Los terminos independientes son:'
do i=1,n
write(*,*)b(i)
end do
!METODO DE JACOBI
k=1
do while (k.le.1000)
do i=1,n
sum1=0.0
sum2=0.0
!REALIZACION DE LAS SUMAS PARA LA OBTENCION DE Xi
do j=1,n
if (j.ne.i) then
sum1=sum1+(a(i,j)*x(j))
end if
end do
if (a(i,i).eq.0) then
write(*,*)'No se puede resolver porque a(i,i)=0'
read(*,*)
stop
end if
x(i)=(-sum1+b(i))/a(i,i)
end do
!OBTENCION DE EL ERROR EN BASE A LA NORMA
do i=1,n
y(i)=x(i)-x0(i)
end do
sum3=0.0
sum4=0.0
do s=1,n
sum3=sum3+(y(s)**2)
sum4=sum4+(x(s)**2)
end do
nor=sum3**0.5
norx=sum4**0.5
! ERROR OBTENIDO SEGUN LAS NORMAS
error=nor/norx
if (error.le.0.0001) then
write(*,*)'La solucion es :'
do i=1,n
write(*,*)x(i)
end do
read(*,*)
stop
else
k=k+1
do i=1,n
x0(i)=x(i)
end do
end if
end do
write(*,*)'Se ha sobrepasado las iteraciones'
stop
end program
b)PROGRAMA EN PYTHON M.GAUSS-SEIDEL
from numpy import*
import numpy
import math
import sys
#funcion que realiza el ingreso de [A] y [b]

def matriz(n):
a=zeros((n,n),float)
b=[]
for i in range (n):
for j in range (n):
a[i][j]=float(input('Ingrese el elemento
'+str(i+1)+';'+str(j+1)+ ' de la matriz A :' ))
for i in range(n):
ind=float(input("*Ingrese el termino independiente "+str(i+1)+"
del sistema [A][x]=[b]:"))
b.append(ind)
return (a,b)
#funcion de el metodo de jacobi, devuelve la solucion

def gauss_seidel(a,b,n,apro):
x0=[]
x=[]
y=[]
itera=1000
tol=0.0001
for i in range(n):
x0.append(apro)
x.append(0.0)
y.append(0.0)
k=1
while k<=itera:
for i in range(n):
sum1=0.0
sum2=0.0
for j in range(0,i):
sum1=sum1+(a[i][j]*x[j])
for h in range(i+1,n):
sum2=sum2+(a[i][h]*x0[h])
if a[i][i]==0:
sys.exit( "No se puede resolver el sistema porque
a[i][i]=0" )
x[i]=(-sum1-sum2+b[i])/a[i][i]
#obtenemos el valor del error relativo segun las normas X ; X0
for i in range(n):
y[i]=x[i]-x0[i]
sum3=0.0
sum4=0.0
for s in range(n):
sum4=sum4+(x[s]**2)
sum3=sum3+(y[s]**2)
nor=abs(math.sqrt(sum3))
norx=abs(math.sqrt(sum4))
norma=nor/norx
# la norma del vector solucion debe ser menor a una tolerancia
if (norma)<=tol:
print 'LA SOLUCION CON ',str(k-1),'ITERACIONES ES :'
for i in range(n):
print 'x',str(i+1),' es ',str(x[i])
sys.exit(' ')
k+=1
for i in range(n):
x0[i]=x[i]
print 'Se ha sobrepasado las ',str(itera),' iteraciones'
return(x)
print 'RESOLUCION DE SISTEMA [A]nxn*[x]n=[b]n '

print 'METODO DE GAUSS-SEIDEL con presicion de 0.0001'
n=int(input('Ingrese el numero de incognitas y ecuaciones: '))
datos=matriz(n)
matri=datos[0]
ter_inde=datos[1]
print 'Sistema a resolver es:'
print 'La Matriz [A] es:'
print matri
print 'el termino independiente [b] es:'
print ter_inde
apro=float(input('Ingrese la aproximacion inicial de Xi: '))
gauss=gauss_seidel(matri,ter_inde,n,apro)
Conclusiones
El método de Crout es uno de los mejores métodos que se puede aplicar a matrices
generales, ya que mediante el mismo podemos esperar resultados bastante aproximados,
existe una pequeña variación en cuanto a los resultados, el programa creado fue optimizado
para la resolución de una matriz banda pentagonal sin hacer operaciones innecesarias.
Por otra parte, si hablamos del método iterativo para la resolución de sistemas lineales
existen dos que se parecen mucho, los cuales son: M. de Jacobi y M. Gauss-Seidel, los
cuales convergen rápidamente a la respuesta mediante k iteraciones, estos métodos son de
mucha utilidad para matrices con muchos ceros.
Bibliografia
− Numerical Analysis – Burden Richard,Faires Douglas,Reynolds Albert - 9th

Edition.
− Analisis-Numerico -Richard-L-Burden -7ma Edicion

Trabajo Final Metodos Numericos (Crout-Jacobi)

Cargado por

Copyright:

Formatos disponibles

Trabajo Final Metodos Numericos (Crout-Jacobi)

Cargado por

Información del documento

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Trabajo Final Metodos Numericos (Crout-Jacobi)

Cargado por

Copyright:

Formatos disponibles

Universidad de Cuenca

Materia: Métodos Numéricos

Profesor: Ing. Fernando Zalamea

Tema: Factorización Código en Fortran y Python.

• METODO CROUT OPTIMIZADO.

Ciclo: Septiembre-Febrero 2014

Cuenca, 24 de enero del 2014.

➢ Realizar el código en Fortran de la factorización de una matriz A en dos matrices L

! no permite el ingreso de matrices con dimension menor a 6

!creacion de la matriz de ceros

!ingreso de los terminos independientes

#la funcion realiza el llenado de la matriz y los terminos independientes

#esta funcion nos devuelve las matrices l y u para utilizarlas en la

print "METODO DE CROUT PARA UNA MATRIZ BANDA DE ANCHO 5"

a)PROGRAMA EN FORTRAN M.JACOBI.

!INGRESO DEL TAMANO DE LA MATRIZ

#funcion que realiza el ingreso de [A] y [b]

#funcion de el metodo de jacobi, devuelve la solucion

print 'RESOLUCION DE SISTEMA [A]nxn*[x]n=[b]n '

− Numerical Analysis – Burden Richard,Faires Douglas,Reynolds Albert - 9th

También podría gustarte