Untitled

446 days ago by pub

def jacobi(A,b,x=None,N=25,epsilon=5): """Resuleve la ecuación Ax=b por el métodod iterativo de Jacobi.""" n=A.nrows() # Creamos solución inicial si no está dada if x is None: x = matrix(RR,n,1) # Creamos las matrices diagonal, D, y las estrictamente # triangulares inferor, L, y superior, U. D = diagonal_matrix(RR,n,A.diagonal()) L = matrix(RR,n,n, lambda i,j: A[i,j] if i>j else 0) U = matrix(RR,n,n, lambda i,j: A[i,j] if i<j else 0) show(html("<b><u>Método de Jacobi</u></b>: $x^{k+1}=-D^{-1}(L+U)x^k +D^{-1}b$<br>")) # Iteramos N veces for i in range(N): sol =-D^(-1)*(L+U)*x+ D^(-1)*b error_absoluto=(sol-x).norm(1) error_relativo=error_absoluto/x.norm(1) "Si alcancamos el error deseado salimos" if error_relativo<=10^(-epsilon): show(html("iteración: " + str(i)+"<br>")) show(html("$\\text{solución aceptada:}"+latex(x)+"$")) show(html("error relativo: " + error_relativo.str())) return x x=sol show(html("iteración: " + str(i))) show(html("$\\text{última solución:}"+latex(x)+"$")) show(html("error relativo: " + error_relativo.str())) show(html("<br><B><U>superado el numero máximo de iteraciones</U></B>")) return x 
       
A=matrix(RR,4,4,[5,2,1,1,1,4,2,1,1,1,3,1,1,2,3,7]) b=matrix(RR,4,1,[1,1,1,1]) x=jacobi(A,b,N=20,x=None,epsilon=8) show(html("<br>==================<br>")) sol=A.change_ring(RR).solve_right(b.change_ring(RR)) show(html("$\\text{solucion exacta:}"+latex(sol)+"$")) show(html("$\\text{diferencia:}"+latex(sol-x)+"$")) show(html("error absoluto: "+(sol-x).norm(1).str()+"<br>")) show(html("error relativo: "+((sol-x).norm(1)/sol.norm(1)).str()+"<br>")) 
       
Método de Jacobi: 
iteración: 19
error relativo: 0.380606900861

superado el numero máximo de iteraciones

==================
error absoluto: 0.104566924158
error relativo: 0.213960014046

                                
                            
Método de Jacobi: 
iteración: 19
error relativo: 0.380606900861

superado el numero máximo de iteraciones

==================
error absoluto: 0.104566924158
error relativo: 0.213960014046

                                
def gauss_seidel(A,b,x=None,N=25,epsilon=5): """Solves the equation Ax=b via the Gauss-Seidel iterative method.""" # Create an initial guess if needed n=A.nrows() if x is None: x = matrix(RDF,n,1) # Create a vector of the diagonal elements of A # and subtract them from A D = diagonal_matrix(RDF,n,A.diagonal()) L = matrix(RDF,n,n, lambda i,j: A[i,j] if i>j else 0) U = matrix(RDF,n,n, lambda i,j: A[i,j] if i<j else 0) show(html("<b><u>Método de Gauss-Seidel</u></b>: $x^{k+1}=-(D+L)^{-1}Ux^k +(D+L)^{-1}b$<br>")) # Iteramos N veces for i in range(N): sol = -(D+L)^(-1)*U*x +(D+L)^(-1)*b error_absoluto=(sol-x).norm(1) error_relativo=error_absoluto/x.norm(1) # Si alcanzamos el error deseado salimos if error_relativo<=10^(-epsilon): show(html("iteración: " + str(i)+"<br>")) show(html("$\\text{solución aceptada: }"+latex(x)+"$")) show(html("error relativo: " + error_relativo.str())) return x x=sol show(html("iteración: " + str(i)+"<br>")) show(html("$\\text{última solución: }"+latex(x)+"$")) show(html("error relativo: " + error_relativo.str())) show(html("<br><B><U>superado el numero máximo de iteraciones</U></B>")) return x 
       
A=matrix(RR,4,4,[5,2,1,1,1,4,2,1,1,1,3,1,1,2,3,7]) b=matrix(RR,4,1,[1,1,1,1]) x=gauss_seidel(A,b,N=20,x=None,epsilon=8) show(html("<br>==================<br>")) sol=A.change_ring(RR).solve_right(b.change_ring(RR)) show(html("$\\text{solucion exacta:}"+latex(sol)+"$")) show(html("$\\text{diferencia:}"+latex(sol-x)+"$")) show(html("error absoluto: "+(sol-x).norm(1).str()+"<br>")) show(html("error relativo: "+((sol-x).norm(1)/sol.norm(1)).str()+"<br>")) 
       
Método de Gauss-Seidel: 
iteración: 16
error relativo: 6.03553715042e-09

==================
error absoluto: 4.13049155296e-09
error relativo: 8.45162117759e-09

                                
                            
Método de Gauss-Seidel: 
iteración: 16
error relativo: 6.03553715042e-09

==================
error absoluto: 4.13049155296e-09
error relativo: 8.45162117759e-09