[已解决]求助大神！做矩阵PALU分解时，遇见了不懂的事情！

wuxianke

def exc_row(B,Q,j):
    n = B.shape[0]
    A = B.copy()
    P = Q.copy()
    for i in range(j,n):
        if A[i,j] >= A[j,j]:
            maxIndex = i   
    tempA = A[maxIndex,:].copy()
    A[maxIndex,:] = A[j,:]
    A[j,:] = tempA
    tempP = P[maxIndex,:].copy()
    P[maxIndex,:] = P[j,:]
    P[j,:] = tempP
    return A,P
   
def palu(A):
    n = A.shape[0]
    U = A.copy()
    L = eye(n)
    P = eye(n)
    for j in range(n-1):
        [U,P] = exc_row(U,P,j)
        for i in range(j+1,n):
            mult = U[i,j]/U[j,j]
            U[i,j] = mult
            print(mult)
            print(U[i,j])
            for k in range(j+1,n):
                U[i,k] = U[i,k] - mult*U[j,k]

    return P,L,U

就是在红色那一段，最后U[i,j]居然等于0，我print之后明明mult都可以显示数值，但是一赋值给U[i,j]就等于0了。

有没有大神能帮忙解决一下这个问题啊

最佳答案

月排行榜 / 总排行榜

BngThea

2018-1-23 13:48:13

wuxianke 发表于 2018-1-23 13:29
已在七楼附上完整代码。

测试了一下，因为你计算的mult都是小数，而你初始化A的时候都是用的整数，所以其被取整，变成了0
请初始化的时候修改A中的某一个值，比如
A = array([[2.0,1,5],[4,4,-4],[1,3,1]])

跳转到最佳答案楼层

BngThea

你的代码是全的？

wuxianke

我又试了下，直接用某个矩阵B，mult = B[2,1]/B[1,1],B[2,1] = mult，B[2,1]还是返回0啊。
这是什么道理啊

DarkmasterSugar

请给出你的完整代码
如果你使用了诸如numpy之类的第三方库，你要在前面使用import导入，不然根本不知道你的shape, eye这样的东西哪里来的
其次，你只给给出了2个函数定义，具体运行时你传入的参数又是什么样的，矩阵A是什么样的，起码要有个测试用的样例
还有，代码请用代码格式发

wuxianke

BngThea 发表于 2018-1-23 12:22
你的代码是全的？

没有，是部分的，稍等一下，我传全部的。

wuxianke

DarkmasterSugar 发表于 2018-1-23 12:56
请给出你的完整代码
如果你使用了诸如numpy之类的第三方库，你要在前面使用import导入，不然根本不知道你 ...

恩恩，好的谢谢提醒。稍等

wuxianke

1. from numpy import dot,array,eye,where,zeros
   
2. 
   
3. def exc_row(B,Q,j):
   
4.     n = B.shape[0]
   
5.     A = B.copy()
   
6.     P = Q.copy()
   
7.     for i in range(j,n):
   
8.         if A[i,j] >= A[j,j]:
   
9.             maxIndex = i   
   
10.     tempA = A[maxIndex,:].copy()
    
11.     A[maxIndex,:] = A[j,:]
    
12.     A[j,:] = tempA
    
13.     tempP = P[maxIndex,:].copy()
    
14.     P[maxIndex,:] = P[j,:]
    
15.     P[j,:] = tempP
    
16.     return A,P
    
17.    
    
18. def palu(A):
    
19.     n = A.shape[0]
    
20.     U = A.copy()
    
21.     L = eye(n)
    
22.     P = eye(n)
    
23.     for j in range(n-1):
    
24.         [U,P] = exc_row(U,P,j)
    
25.         for i in range(j+1,n):
    
26.             mult = U[i,j]/U[j,j]
    
27.             U[i,j] = mult
    
28.             print(mult)
    
29.             print(U[i,j])
    
30.             for k in range(j+1,n):
    
31.                 U[i,k] = U[i,k] - mult*U[j,k]
    
32.     return P,L,U    #这里的L还没有写，还在测试刚刚提出的问题
    
33.         
    
34. def backsolve(U,b):
    
35.     n = U.shape[0]
    
36.     x = zeros((n,1))
    
37.     for i in range(n-1,-1,-1):
    
38.         if i != n-1:
    
39.             for k in range(i+1,n-1):
    
40.                 b[i] = b[i] - x[i]*U[i,k]
    
41.         x[i] = b[i]/U[i,i]
    
42.         
    
43. def l2u(L):
    
44.     n = L.shape[0]
    
45.     Lu = eye(n)
    
46.     for j in range(n):
    
47.         for i in range(j,n):
    
48.             Lu[n-1-i,n-1-j] = L[i,j]
    
49.     return Lu
    
50. 
    
51. def one_reverse(b):
    
52.     n = b.shape[0]
    
53.     b_rev = zeros((n,1))
    
54.     for i in range(n):
    
55.         b_rev[i,0] = b[n-1-i,0]
    
56.     return b_rev   
    
57. 
    
58. def palusolve(A,b):
    
59.     [P,L,U] = palu(A)
    
60.     Lu = l2u(L)   
    
61.     c = backsolve(Lu,one_reverse(dot(P,b)))
    
62.     x = backsolve(U,one_reverse(c))
    
63.     return x

复制代码

测试的矩阵是：A = array([[2,1,5],[4,4,-4],[1,3,1]])，测试的函数是[P,L,U] = palu(A)
输出的U是：array([[ 4,  4, -4],[ 0,  2,  2],[ 0,  0,  8]])；
就是那几个0的地方，我是进行了U[i,j] = mult操作的，而mult是不等零的。不知道是哪里的问题，所以向各位请教了

wuxianke

BngThea 发表于 2018-1-23 12:22
你的代码是全的？

已在七楼附上完整代码。

wuxianke

DarkmasterSugar 发表于 2018-1-23 12:56
请给出你的完整代码
如果你使用了诸如numpy之类的第三方库，你要在前面使用import导入，不然根本不知道你 ...

已在七楼附上完整代码。

BngThea

wuxianke 发表于 2018-1-23 13:29
已在七楼附上完整代码。

测试了一下，因为你计算的mult都是小数，而你初始化A的时候都是用的整数，所以其被取整，变成了0
请初始化的时候修改A中的某一个值，比如
A = array([[2.0,1,5],[4,4,-4],[1,3,1]])

wuxianke

BngThea 发表于 2018-1-23 13:48
测试了一下，因为你计算的mult都是小数，而你初始化A的时候都是用的整数，所以其被取整，变成了0
请初始 ...

啊，原来如此，谢谢谢谢！