numpy 线性代数求解器

Question

让我直截了当地描述一下我试图在 python 中解决的简单数学问题：

我有 4 个 nxn 刚度系数矩阵，即 s11、s12、s21、s22，因此我试图求解的代数方程组如下所示：

[s11]*{u}+[s12]*{v} = {fx}
[s22]*{v}+[s21]*{u} = {fy}

其中 u 和 v - 是 x 和 y 方向上未知向量的 nx1 个分量。 IE。总向量将为 g = sqrt(u^2+v^2)。 fx 和 fy - 是给定“力”向量的 nx1 x 和 y 分量。

所以任务相对容易，但我坚持使用 python 中的 numpy 实现高斯消除。我尝试使用的第一个也是最明显的方法之一是以以下格式重写所有内容：

s = [ [[s11], [s12]],
      [[s21], [s22]] ]
f = [ [fx],
      [fy] ]

使得 s 是 2nx2n，f 是 2nx1。组装方式如下：

s1 = np.concatenate((s11, s12), axis=1)
s2 = np.concatenate((s21, s22), axis=1)
s = np.concatenate((s1, s2), axis=0)
f = np.concatenate((fx, fy), axis=0)

但是，当我尝试用

解决系统问题时

u = np.linalg.solve(s, f)

我希望得到 2nx1 未知向量 u（其中第一组 n 个元素是 x 分量，第二组是 y），但我收到以下错误：

Traceback (most recent call last):
File "/home/art/PycharmProjects/Cavern/main.py", line 204, in <module>
u = np.linalg.solve(s, f)
File "<__array_function__ internals>", line 5, in solve
File "/home/art/.pyenv/versions/test380/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 399, in solve
r = gufunc(a, b, signature=signature, extobj=extobj)
File "/home/art/.pyenv/versions/test380/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 97, in _raise_linalgerror_singular
raise LinAlgError("Singular matrix")
numpy.linalg.LinAlgError: Singular matrix

然后我又尝试了另一种有争议的方式，可以总结如下：

s1 = s11 + s21
s2 = s12 + s22
s = np.concatenate((s1, s2), axis=1)
f = fx + fy

使得 s 现在是 nx2n，f 是 nx1。再说一次，我尝试使用高斯消除如下：

u = np.linalg.solve(s, f)

我再次希望得到 2nx1 的未知向量，但我得到了以下错误：

Traceback (most recent call last):
File "/home/art/PycharmProjects/Cavern/main.py", line 204, in <module>
u = np.linalg.solve(s, f)
File "<__array_function__ internals>", line 5, in solve
File "/home/art/.pyenv/versions/test380/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 386, in solve
_assert_stacked_square(a)
File "/home/art/.pyenv/versions/test380/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 213, in _assert_stacked_square
raise LinAlgError('Last 2 dimensions of the array must be square')
numpy.linalg.LinAlgError: Last 2 dimensions of the array must be square

因此，我渴望寻求您的帮助。任何帮助/提示将不胜感激。是的，mb 这并不重要，但整件事就是将我的代码从 Matlab 转移到 python。在 Matlab 中的过程如下：

S1 = [S11, S12]; S2 = [S21, S22];
U1 = S1 \ fx; U2 = S2 \ fy;
U = U1 + U2;
u(:,1) = U(1:nnodes);
v(:,1) = U(nnodes+1:end);

反斜杠在这里代表高斯消除。

Answer 1

如果 s 和 f 系数是标量，为什么要使用 np.concatenate？您可以直接将它们制作成 numpy 数组（或任何有序序列）并在它们上使用numpy.linalg.solve()。

s1 = [s11, s12]
s2 = [s21, s22]
s = np.column_stack((s1, s2))
f = [s11, s12]

u = np.linalg.solve(s, f)

注意在这种情况下，一维列表与具有 nx1（或 1xn）向量相同。向量的方向性与np.linalg.solve() 函数无关。