时间戳

首页 > 解决方案 > Numpy 矩阵求逆对大于 18 的阶数给出错误结果

python - Numpy 矩阵求逆对大于 18 的阶数给出错误结果

问题描述

我在用 反转矩阵时遇到问题NumPy。奇怪的是，它只在订单 18 之前给出正确的结果。一旦订单大于 18，它就会给出错误的结果。

import numpy as np
from decimal import Decimal
import numpy.matlib

I_1=np.matlib.eye(ngrid,ngrid,k=0,dtype=Decimal)
I_2=np.matlib.eye(ngrid,ngrid,k=1,dtype=Decimal)
I_3=np.matlib.eye(ngrid,ngrid,k=2,dtype=Decimal)

B=I_1 + 10.*I_2 + I_3
B=np.divide(B,12.)

B_inv=np.linalg.inv(B)
print B_inv

C=B.dot(B_inv)
print C

最后一行被包括在内以检查它是否给出了正确的结果。

标签： pythonnumpy

解决方案

您的代码在技术上没有任何问题。但是，您在数值分析中遇到了一个非常基本的问题。这里有两个效果：

  1. 浮点错误
  2. 问题的条件编号

我不会在这里解释它们，因为几乎任何关于数值分析的介绍性书籍（例如Kendall Atkinson）都比我在这里做得更好。我只会给你留下这个代码片段：

import numpy as np

NGRID = 8


def HilbertMatrix(N, dtype=np.float64):
    arr = np.zeros((N,N), dtype=dtype)
    for i in range(N):
        for j in range(N):
            arr[i,j] = 1 / (i + j + 1)
    return arr


H = HilbertMatrix(NGRID)
J = np.linalg.inv(H)

C = np.dot(H, J)

print(np.allclose(C, np.eye(NGRID)))
print(np.linalg.norm(C))

尝试使用 NGRID，看看 C 离单位矩阵有多远。

推荐阅读