EJ PRACTICA 2

610 days ago by u49725808

 
       
 
       
 
       
 
       
 
       
 
       
dig=4 A=matrix(QQ,4,4,lambda i,j: -3 if i==j else 1 if abs(i-j)==1 else 0) b=matrix(QQ,4,1,lambda i,j: i+1) S=(A\b).n(digits=dig) show(html("$\\text{solucion = }"+ latex(S.transpose()) +"$")) print "GAUSS con pivoteo distinto de cero" print "==================================" sol=gauss(A,b,digitos=dig) show(html("$A="+latex(A.n(digits=dig))+"; b="+latex(b.n(digits=dig))+"\\text{. sol=}"+latex(sol.n(digits=dig))+"$")) print "error abosluto = " + (S-sol).norm(oo).str() print "error relativo = " + (100*(S-sol).norm(oo)/S.norm(oo)).str()+"%" print "GAUSS con pivoteo parcial" print "=========================" sol=gaussparcial(A,b,digitos=dig) show(html("$A="+latex(A.n(digits=dig))+"; b="+latex(b.n(digits=dig))+"\\text{. sol=}"+latex(sol.n(digits=dig))+"$")) print "error abosluto = " + (S-sol).norm(oo).str() print "error relativo = " + (100*(S-sol).norm(oo)/S.norm(oo)).str()+"%" print "GAUSS con pivoteo escalado" print "=========================" sol=gaussescalado(A,b,digitos=dig) show(html("$A="+latex(A.n(digits=dig))+"; b="+latex(b.n(digits=dig))+"\\text{. sol=}"+latex(sol.n(digits=dig))+"$")) print "error abosluto = " + (S-sol).norm(oo).str() print "error relativo = " + (100*(S-sol).norm(oo)/S.norm(oo)).str()+"%" 
       

GAUSS con pivoteo distinto de cero
==================================
determinante: 55.00

error abosluto = 0.0001220703125
error relativo = 0.0053711462026%
GAUSS con pivoteo parcial
=========================
determinante: 55.00

error abosluto = 0.0001220703125
error relativo = 0.0053711462026%
GAUSS con pivoteo escalado
=========================

determinante: 55.00

error abosluto = 0.0001220703125
error relativo = 0.0053711462026%
GAUSS con pivoteo distinto de cero
==================================
determinante: 55.00

error abosluto = 0.0001220703125
error relativo = 0.0053711462026%
GAUSS con pivoteo parcial
=========================
determinante: 55.00

error abosluto = 0.0001220703125
error relativo = 0.0053711462026%
GAUSS con pivoteo escalado
=========================

determinante: 55.00

error abosluto = 0.0001220703125
error relativo = 0.0053711462026%
dig=4 A=matrix(QQ,20,20,lambda i,j: -3 if i==j else 1 if abs(i-j)==1 else 0) b=matrix(QQ,20,1,lambda i,j: i+1) S=(A\b).n(digits=dig) show(html("$\\text{solucion = }"+ latex(S.transpose()) +"$")) print "GAUSS con pivoteo distinto de cero" print "==================================" sol=gauss(A,b,digitos=dig) show(html("$A="+latex(A.n(digits=dig))+"; b="+latex(b.n(digits=dig))+"\\text{. sol=}"+latex(sol.n(digits=dig))+"$")) print "error abosluto = " + (S-sol).norm(oo).str() print "error relativo = " + (100*(S-sol).norm(oo)/S.norm(oo)).str()+"%" print "GAUSS con pivoteo parcial" print "=========================" sol=gaussparcial(A,b,digitos=dig) show(html("$A="+latex(A.n(digits=dig))+"; b="+latex(b.n(digits=dig))+"\\text{. sol=}"+latex(sol.n(digits=dig))+"$")) print "error abosluto = " + (S-sol).norm(oo).str() print "error relativo = " + (100*(S-sol).norm(oo)/S.norm(oo)).str()+"%" print "GAUSS con pivoteo escalado" print "=========================" sol=gaussescalado(A,b,digitos=dig) show(html("$A="+latex(A.n(digits=dig))+"; b="+latex(b.n(digits=dig))+"\\text{. sol=}"+latex(sol.n(digits=dig))+"$")) print "error abosluto = " + (S-sol).norm(oo).str() print "error relativo = " + (100*(S-sol).norm(oo)/S.norm(oo)).str()+"%" 
       

GAUSS con pivoteo distinto de cero
==================================
determinante: 2.679e8

error abosluto = 0.001953125
error relativo = 0.0116049669258%
GAUSS con pivoteo parcial
=========================
determinante: 2.679e8

error abosluto = 0.001953125
error relativo = 0.0116049669258%
GAUSS con pivoteo escalado
=========================

determinante: 2.679e8

error abosluto = 0.001953125
error relativo = 0.0116049669258%
GAUSS con pivoteo distinto de cero
==================================
determinante: 2.679e8

error abosluto = 0.001953125
error relativo = 0.0116049669258%
GAUSS con pivoteo parcial
=========================
determinante: 2.679e8

error abosluto = 0.001953125
error relativo = 0.0116049669258%
GAUSS con pivoteo escalado
=========================

determinante: 2.679e8

error abosluto = 0.001953125
error relativo = 0.0116049669258%
 
       
 
       
 
       
 
       
 
       
 
       
 
       
 
       
 
       
# Método de Gauss con pivoteo escalado por columanas para Sistemas de Ecuaciones Lineales # con matriz de los coeficientes cuadrada def gaussescalado(A,b,anillo=QQ,digitos=53): m,n = A.ncols(),A.nrows() tam = b.nrows() if ((m<>n) or (tam <> n)): print 'Las dimensiones del sistema no son válidas' return False #factores para escalado A.change_ring(anillo) b.change_ring(anillo) S=[max(abs(A[i,j]) for j in range(n)) for i in range(n)] show(html("$\\text{factores=}" + latex(S) +"$")) A=A.n(digits=digitos).augment(b.n(digits=digitos)) alpha=0 #indice de permutaciones de filas for i in range(n): #Si el pivote es nulo se realiza un cambio de fila pivote=abs(A[i,i]/S[i]) filapivote=i for s in range((i + 1),n): if abs(A[s, i])/S[s]> pivote: filapivote=s pivote=abs(A[s, i]/S[s]) if (pivote==0): #no se ha encontrado enu pivote distinto de cero print 'el sistema no tiene solución única o es incompatibe' return False if (filapivote<>i): A[i],A[filapivote]=A[filapivote],A[i] alpha=alpha+1 for j in range((i + 1),n): A[j] = A[j] - (A[i]*A[j, i]/A[i, i]) #calculamos determinate print "determinante: " +((-1.0)^alpha*prod(A[i,i] for i in range(0,n))).str() #recuperamos matices Aux=A.delete_columns([n]) #extraemos A^n b=A.delete_columns(range(0,n)) #extraemos b^n A=Aux # metodo de la subida x = matrix(anillo,n,1) x[n-1] = b[n-1]/A[n-1,n-1] for k in range(n-2,-1,-1): suma = sum(x[i]*A[k,i] for i in range(n-1,k,-1)) x[k] = (b[k]-suma)/A[k,k] return x 
       
# Método de Gauss con pivoteo parcial para Sistemas de Ecuaciones Lineales # con matriz de los coeficientes cuadrada def gaussparcial(A,b,anillo=QQ,digitos=53): m,n = A.ncols(),A.nrows() tam = b.nrows() A.change_ring(anillo) b.change_ring(anillo) if ((m<>n) or (tam <> n)): print 'Las dimensiones del sistema no son válidas' return False A=A.n(digits=digitos).augment(b.n(digits=digitos)) alpha=0 #indice de permutaciones de filas for i in range(n): #Si el pivote es nulo se realiza un cambio de fila pivote=abs(A[i, i]) filapivote=i for s in range((i + 1),n): if abs(A[s, i])> pivote: filapivote=s pivote= A[s, i] if (pivote==0): #no se ha encontrado un pivote distinto de cero print 'el sistema no tiene solución única o es incompatibe' return False if (filapivote<>i): A[i],A[filapivote]=A[filapivote],A[i] alpha=alpha+1 for j in range((i + 1),n): A[j] = A[j] - (A[i]*A[j, i]/A[i, i]) #calculamos determinate print "determinante: " +((-1.0)^alpha*prod(A[i,i] for i in range(0,n))).str() #recuperamos matices Aux=A.delete_columns([n]) #extraemos A^n b=A.delete_columns(range(0,n)) #extraemos b^n A=Aux # metodo de la subida x = matrix(QQ,n,1) x[n-1] = b[n-1]/A[n-1,n-1] for k in range(n-2,-1,-1): suma = sum(x[i]*A[k,i] for i in range(n-1,k,-1)) x[k] = (b[k]-suma)/A[k,k] return x 
       
# Método de Gauss para Sistemas de Ecuaciones Lineales # con matriz de los coeficientes cuadrada def gauss(A,b,anillo=QQ,digitos=53): m,n = A.ncols(),A.nrows() tam = b.nrows() if ((m<>n) or (tam <> n)): print 'Las dimensiones del sistema no son válidas' return False A.change_ring(anillo) b.change_ring(anillo) A=A.n(digits=digitos).augment(b.n(digits=digitos)) alpha=0 #indice de permutaciones de filas for i in range(n): #Si el pivote es nulo se realiza un cambio de fila if (A[i, i] == 0): encontrado = True for s in range((i + 1),n): if A[s, i]<> 0: encontrado = False aux1 = A[i] A[i] = A[s] A[s] = aux1 alpha=alpha+1 if encontrado: #no se ha encontrado un pivote distinto de cero print 'el sistema no tiene solución única o es incompatibe' return False for j in range((i + 1),n): A[j] = A[j] - (A[i]*A[j, i]/A[i, i]) #calculamos determinate print "determinante: " +((-1.0)^alpha*prod(A[i,i] for i in range(0,n))).str() #recuperamos matices Aux=A.delete_columns([n]) #extraemos A^n b=A.delete_columns(range(0,n)) #extraemos b^n A=Aux # metodo de la subida x = matrix(QQ,n,1) x[n-1] = b[n-1]/A[n-1,n-1] for k in range(n-2,-1,-1): suma = sum((x[i]*A[k,i]) for i in range(n-1,k,-1)) x[k] = (b[k]-suma)/A[k,k] return x 
       
# Metodo de sustitucion progresiva o bajada para # sistemas de ecuaciones lineales triangulares inferiores def bajada(A, b,anillo=QQ,digitos=53): A.change_ring(anillo) b.change_ring(anillo) # comprobamos validez del sistema m,n = A.nrows(),A.ncols() tam = b.nrows() if ((m <> n) or (tam <> m)): print 'Las dimensiones del sistema no son válidas' return True if A.det()==0: print 'El sistema no tiene solución única o es incompatible' return True # Vector de soluciones x = matrix(anillo,n,1) x[0] = b[0].n(digits=digitos)/A[0,0].n(digits=digitos) for k in range(1,n): suma = sum(x[i]*A[k,i].n(digits=digitos) for i in range(0,k)) x[k] = (b[k].n(digits=digitos)-suma)/A[k,k].n(digits=digitos) return x 
       
# Metodo de sustitucion regresiva o subida para # sistemas de ecuaciones lineales triangulares superiores def subida(A, b,anillo=QQ,digitos=53): A.change_ring(anillo) b.change_ring(anillo) # comprobamos validez del sistema m,n = A.nrows(),A.ncols() tam = b.nrows() if ((m <> n) or (tam <> m)): print 'Las dimensiones del sistema no son válidas' return True if A.det()==0: print 'El sistema no tiene solución única o es incompatible' return True # Vector de soluciones x = matrix(anillo,n,1) x[n-1] = b[n-1].n(digits=digitos)/A[n-1,n-1].n(digits=digitos) for k in range(n-2,-1,-1): suma = sum(x[i]*A[k,i].n(digits=digitos) for i in range(n-1,k,-1)) x[k] = (b[k].n(digits=digitos)-suma)/A[k,k].n(digits=digitos) return x