修改Python结构分析代码遇numpy.linalg.LinAlgError错误求助

问题概述 #

修改Oliver Natt所著《Physik mit Python》中的结构分析代码以计算复杂结构时，添加新节点后出现两类错误：一是求解方程组时触发numpy.linalg.LinAlgError，提示系数矩阵非方阵；二是直接修改原书代码时出现索引越界错误。常规语法检查和结构调整未解决问题，推测为数学层面的方程组适定性问题。

报错代码 #

import numpy as np

#Set the scaling factor for the force vectors [m/N].
force_scale = 0.005

#Number of dimensions for the problem (2 or 3).
dim = 3

#Set the positions of the points [m].
points = np.array([
[0, 0, 0], [1.5, 0, 0], [0, 1.5, 0], [0, 0, 2], [1, 0, 2],
[1, 1, 2], [0, 1, 2], [1.5, 1.5, 0]
])

#Create a list of the indices of the support points.
support_indices = [0, 1, 2, 7]

#Each bar connects exactly two points. We store the indices
#of the corresponding points in an array.
bars = np.array([
[0, 3], [1, 4], [2, 6], [3, 4], [4, 5], [5, 6], [6, 3],
[3, 5], [0, 4], [0, 6], [0, 5], [1, 5], [7, 5]
])

#Set the external force acting on each point [N].
#For the support points, we initially set this force to 0.
#This will be calculated later.
external_force = np.array([
[0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, -98.1], [0, 0, -98.1],
[0, 0, -98.1], [0, 0, -98.1], [0, 0, 0]
])

#Define the number of points, bars, etc.
n_points = points.shape[0]
n_bars = bars.shape[0]
n_support = len(support_indices)
n_nodes = n_points - n_support
n_equations = n_nodes * dim

#Create a list of the indices of the nodes.
node_indices = list(set(range(n_points)) - set(support_indices))

def unit_vector(point_index, bar_index):
    """Returns the unit vector pointing from the point with index
    point_index along the bar with index bar_index. """
    i1, i2 = bars[bar_index]
    if point_index == i1:
        vec = points[i2] - points[i1]
    else:
        vec = points[i1] - points[i2]
    return vec / np.linalg.norm(vec)

# Set up the equation system for the forces.
A = np.zeros((n_equations, n_bars))
for i, bar in enumerate(bars):
    for k in np.intersect1d(bar, node_indices):
        n = node_indices.index(k)
        A[n * dim:(n + 1) * dim, i] = unit_vector(k, i)

#Solve the equation system A @ F = -external_force for the forces F.
b = -external_force[node_indices].reshape(-1)
forces = np.linalg.solve(A, b)

#Calculate the external forces.
for i, bar in enumerate(bars):
    for k in np.intersect1d(bar, support_indices):
        external_force[k] -= forces[i] * unit_vector(k, i)

可正常运行的代码 #

points = np.array([
[0, 0, 0], [1.5, 0, 0], [0, 1.5, 0], [0, 0, 2], [1, 0, 2],
[1, 1, 2], [0, 1, 2]
])

support_indices = [0, 1, 2]

bars = np.array([
[0, 3], [1, 4], [2, 6], [3, 4], [4, 5], [5, 6], [6, 3],
[3, 5], [0, 4], [0, 6], [0, 5], [1, 5]
])

external_force = np.array([
[0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, -98.1], [0, 0, -98.1],
[0, 0, -98.1], [0, 0, -98.1]
])

错误回溯信息 #

线性代数求解错误 #

Traceback (most recent call last):
  File "C:\...\main.py", line 62, in <module>
    forces = np.linalg.solve(A, b)
  File "<__array_function__ internals>", line 200, in solve
  File "C:\...\venv\lib\site-packages\numpy\linalg\linalg.py", line 373, in solve
    _assert_stacked_square(a)
  File "C:\...\venv\lib\site-packages\numpy\linalg\linalg.py", line 190, 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

索引越界错误（原书代码修改后） #

Traceback (most recent call last):
  File "C:\...\original.py", line 66, in <module>
    A[n * dim:(n + 1) * dim, i] = einheitsvektor(k, i)
IndexError: index 12 is out of bounds for axis 1 with size 12

错误原因分析 #

1. 线性代数求解错误核心原因：
   np.linalg.solve仅支持求解适定方程组（方程数=未知数个数，即系数矩阵为方阵）。修改后的参数计算：

   • 自由节点数n_nodes = 8 - 4 = 4，3维问题下方程数n_equations = 4*3 = 12
   • 杆件数n_bars = 13，因此系数矩阵A的形状为(12,13)，非方阵，不满足solve的要求。
     而正常运行的代码中，n_bars=12，A为(12,12)方阵，符合要求。

2. 索引越界错误原因：
   原书代码中einheitsvektor（对应修改后的unit_vector）返回的向量维度与矩阵切片不匹配，或节点索引计算逻辑在新增节点后出现偏移，导致访问了超出矩阵列数的索引。

解决办法 #

针对线性代数求解错误 #

有两种方案可选：

• 方案1：调整结构，使方程组适定
  减少一根杆件，让n_bars = n_equations = 12（例如删除杆件[7,5]），此时A为方阵，可继续使用np.linalg.solve。
• 方案2：使用最小二乘法求解超定方程组
  若需保留所有杆件，将求解语句替换为最小二乘法，该方法可处理超定（方程数<未知数个数）或欠定方程组：

  # 替换原solve语句
  forces, residuals, rank, singular_values = np.linalg.lstsq(A, b, rcond=None)
  

  此方法会返回使误差平方和最小的杆件力解。

针对索引越界错误 #

• 确保node_indices的生成逻辑正确，可打印node_indices确认自由节点索引无误；
• 检查unit_vector函数返回的向量维度是否为dim=3，确保与矩阵切片A[n*dim:(n+1)*dim, i]的行数匹配；
• 原书代码中的einheitsvektor函数若存在硬编码的维度或索引，需同步修改为动态计算（如使用dim变量而非固定值）。

内容的提问来源于stack exchange，提问作者klmo6000