python 实现matlab左除
import numpy as np
import scipy
def mldivide_new(A,b):
[mm,nn]=np.shape(A)
ss=np.shape(b)[1]
x=np.zeros((nn,ss))
if np.shape(A)[0]==np.shape(A)[1]:
if A==np.tril(A):
for ii in range(ss):
x[:,ii]=backward_sub(A,b[:,ii])
return
elif A==np.triu(A):
for ii in range(ss):
x[:,ii]=backward_sub(A,b[:,ii])
return
else:
if A==A.T:
[R,p]=np.linalg.cholesky(A)
if (p==0):
for ii in range(ss):
x[:,ii]=backward_sub(R,forward_sub(R.T,b[:,ii]))
return
[L,U,P]=scipy.linalg.lu(A)
for ii in range(ss):
x[:,ii]=backward_sub(U,forward_sub(L,P.dot(b[:,ii])))
return
else:
print(np.shape(A))
[Q,R]=np.linalg.qr(A)
Q=np.hstack((Q,np.array([[ii] for ii in Q[:,-1]])))
R=np.vstack((R,R[-1,:]))
Q[np.abs(Q)<=2.261143e-14]=0
R[np.abs(R)<=2.261143e-14]=0
print(np.shape(Q))
print(np.shape(R))
print(Q)
print("R=", R)
for jj in range(ss):
print(np.shape(x))
x1 = backward_sub(R, Q.T.dot(b[:, jj]))
x1=[i[0] for i in x1]
x[:,jj]=x1
return x
return x
def forward_sub(A,b):
[m,n]=np.shape(A)
x=np.zeros((n,1))
for i in range(0,n):
if np.abs(A[i,i])>=2.261143e-14:
x[i]=(b[i]-A[i,0:i-1].dot(x[0,i-1]))/A[i,i]
elif b[i]-A[i,0:i-1].dot(x[0:i-1])==0:
x[i]=0
return x
def backward_sub(A,b):
[m, n] = np.shape(A)
x = np.zeros((n, 1))
for i in range(n-1,-1,-1):
print(int(A[i,i]))
print(type(A[i,i]))
if int(A[i,i])!=0:
x[i]=(b[i]-A[i,i+1:-1].dot(x[i+1:-1]))/A[i,i]
elif b[i]-A[i,i+1:-1].dot(x[i+1:-1])==0:
x[i]=0
return x
if __name__=='__main__':
A=np.ones((10,9))*10
B=np.ones((10,3))*3
x=mldivide_new(A,B)
print(x)
这里面应该有bug,请大神们,指点一二。
点赞
收藏
分享
文章举报
qestion_yz_10086
发布了69 篇原创文章 · 获赞 2 · 访问量 1662
私信
关注