Untitled

895 days ago by pub

def gauss_seidel(A,b,x=None,N=15,epsilon=7): """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,[4,1,-1,1,1,4,-1,-1,-1,-1,5,1,1,-1,1,3]) b=matrix(RR,4,1,[-2,-1,0,1]) x=gauss_seidel(A,b,N=15,x=None,epsilon=4) 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: 9
error relativo: 7.7958455154e-05

==================
error absoluto: 0.00021364504874
error relativo: 0.000120899911303

                                
                            
Método de Gauss-Seidel: 
iteración: 9
error relativo: 7.7958455154e-05

==================
error absoluto: 0.00021364504874
error relativo: 0.000120899911303

                                
def gauss_seidel_rel(A,b,x=None,N=15,epsilon=7,omega=None): """Solves the equation Ax=b via the Jacobi 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) if omega is None: G=-(D+L)^(-1)*U #calulo de omega a partir del radio espectral radio_espectral=max(G.eigenvalues()) omega=2/(1+sqrt(1-radio_espectral^2)) show(html("<b><u>Método de Gauss-Seidel con relajación</u></b>: $x^{k+1}=(D+\\omega L)^{-1}[(1-\\omega )D-\\omega U]x^k +\\omega (D+\\omega L)^{-1}b$<br>")) show(html("<b><u>omega</u></b>: "+str(omega)+"<br>")) # Iteramos N veces for i in range(N): sol = (D+omega*L)^(-1)*((1-omega)*D-omega*U)*x +omega*(D+omega*L)^(-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)+"<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,[4,1,-1,1,1,4,-1,-1,-1,-1,5,1,1,-1,1,3]) b=matrix(RR,4,1,[-2,-1,0,1]) x=gauss_seidel(A,b,N=15,x=None,epsilon=4) show(html("<br>==================<br>")) x=gauss_seidel_rel(A,b,N=15,x=None,epsilon=4) show(html("<br>==================<br>")) sol=A.solve_right(b) 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: 9
error relativo: 7.7958455154e-05

==================
Método de Gauss-Seidel con relajación: 
omega: 1.03371531225
iteración: 8
error relativo: 7.82231903152e-05

==================
error absoluto: 0.00019846787169
error relativo: 0.000112311276228

                                
                            
Método de Gauss-Seidel: 
iteración: 9
error relativo: 7.7958455154e-05

==================
Método de Gauss-Seidel con relajación: 
omega: 1.03371531225
iteración: 8
error relativo: 7.82231903152e-05

==================
error absoluto: 0.00019846787169
error relativo: 0.000112311276228