数值线性代数求解器

Question

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

我有4个n×n的系数矩阵，即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)

请注意，在这种情况下，1D列表与具有nx1 (或1xn)向量的列表相同。矢量的方向性对于np.linalg.solve()函数并不重要。