NUMPY 线性代数求解器

Question

让我直接进入交易，并描述我试图用python解决的简单数学问题：

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

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

其中 u 和 v - 是 x 和 y 方向上未知向量的 nx1 分量。 即总向量将为 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)

我希望得到未知数 u 的 2nx1 向量(其中第一组 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()函数无关紧要。