Untitled

509 days ago by pub

def m_QR(A): m,n= A.nrows(),A.ncols() A= copy(A) Q= matrix(RR,m,n) R= matrix(RR,n,n) for k in range(n): R[k,k]= A.matrix_from_columns([k]).norm(p='frob') Q[:,k:k+1]= A[:,k:k+1]/R[k,k] R[k:k+1,k+1:n+1]= Q[:,k:k+1].transpose()* A[:,k+1:n+1] A[:,k+1:n+1]= A[:, k+1:n+1] - Q[:,k:k+1]*R[k:k+1,k+1:n+1] return Q,R # Metodo de sustitucion regresiva o subida para # sistemas de ecuaciones lineales triangulares superiores def subida(A, b): # 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(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 
       
X=[0,0.1,0.2,0.3,0.4,0.5] Y=[1.5,1,-0.5,0,0.5,-1] grado=2 fig1= point(zip(X,Y),color='blue',size=20,legend_label="$(x,y)$") A=matrix(RDF,len(X),grado+1,lambda i,j: X[i]^j) b=matrix(RDF,len(Y),1,lambda i,j: Y[i]) Q,R=m_QR(A) x=subida(R,Q.transpose()*b) show(html("$\\text{Ax - b = r} \Rightarrow "+latex(A.n(digits=5))+latex(x.n(digits=5))+"-"+latex(b.n(digits=5))+"="+latex((A*x-b).n(digits=5))+"$")) show(html("$\\text{QRx = b } \Rightarrow " + latex(Q)+latex(R)+latex(x.n(digits=5))+"="+latex(b)+ "$")) show(html("$\\text{Rx=Q}^t\\text{b} \Rightarrow "+latex(R.n(digits=5))+latex(x.n(digits=5))+"="+latex((Q.transpose()*b).n(digits=5))+"$")) var('t') f=sum(x[i,0].n(digits=4)*t^i for i in range(grado+1)) leg="$y="+latex(f)+"$" show(html(leg)) show(html("<br>=================<br>")) fig2=plot(f,min(X),max(X),color="red",legend_label=leg) show(fig1+fig2,figsize=8) 
       

=================

                                
                            

=================