用numpy计算零空间的有理基

2024-10-05 14:31:25 发布

您现在位置:Python中文网/ 问答频道 /正文

我试图计算一个矩阵的零空间的有理基。关于如何使用Python/numpy计算nullspace的文章有很多,但是他们计算的是正交基,而不是有理基。以下是如何在MATLAB中实现的:

ns = null(A,'r')

当我查看源代码时,我看到它是这样计算的:

function Z = null(A,how)
[m,n] = size(A)
%...
[R,pivcol] = rref(A);
r = length(pivcol);
nopiv = 1:n;
nopiv(pivcol) = [];
Z = zeros(n,n-r,class(A));
if n > r
  Z(nopiv,:) = eye(n-r,n-r,class(A));
  if r > 0
     Z(pivcol,:) = -R(1:r,nopiv);
  end
end
%...

function [A,jb] = rref(A,tol)
%...
[m,n] = size(A);

[num, den] = rat(A);
rats = isequal(A,num./den);

if (nargin < 2), tol = max(m,n)*eps(class(A))*norm(A,'inf'); end

i = 1;
j = 1;
jb = [];
while (i <= m) && (j <= n)

   [p,k] = max(abs(A(i:m,j))); k = k+i-1;
   if (p <= tol)

      A(i:m,j) = zeros(m-i+1,1);
      j = j + 1;
   else

      jb = [jb j];

      A([i k],j:n) = A([k i],j:n);

      A(i,j:n) = A(i,j:n)/A(i,j);

      for k = [1:i-1 i+1:m]
         A(k,j:n) = A(k,j:n) - A(k,j)*A(i,j:n);
      end
      i = i + 1;
      j = j + 1;
   end
end

if rats
    [num,den] = rat(A);
    A=num./den;
end

这里rref是简化的行梯队形式。因此,通过查看此源代码,我尝试用以下代码重新创建它:

^{pr2}$

有两件事我不知道。首先,我不知道如何将比率部分添加到rref函数中。其次,我不确定我的索引是否正确,因为MATLAB的索引是从1开始的,当你选择一个切片时,索引包含最后一个元素(即1:5同时包含1和5)。在


Tags: sizeif源代码functionnullnumclassend
1条回答
网友
1楼 · 发布于 2024-10-05 14:31:25

SymPy开箱即用,尽管(作为符号,在Python中)不如NumPy或Scipy快。浮点输入示例:

from sympy import Matrix, S, nsimplify
M = Matrix([[2.75, -1.2, 0, 3.2], [8.29, -4.8, 7, 0.01]])
print(nsimplify(M, rational=True).nullspace())

打印两列向量的列表,用一列矩阵表示。在

^{pr2}$

必须使用nsimplify将浮点数转换为它们要表示的有理数。如果矩阵是作为整数/有理数项的矩阵创建的,那就没有必要了。在

M = Matrix([[1, 2, 3, 5, 9], [9, -3, 0, 2, 4], [S(3)/2, 0, -1, 2, 0]])
print(M.nullspace())

[Matrix([
[ -74/69],
[-176/69],
[   9/23],
[      1],
[      0]]), Matrix([
[ -70/69],
[-118/69],
[ -35/23],
[      0],
[      1]])]

这里,S(3)/2代替了`3/2,以强制创建SymPy对象,而不是浮点求值。在

相关问题 更多 >