超静定梁有限差分法求解异常：节点数增加后结果偏差

算力阿峰

2026-6-18

有限差分法求解超静定梁的位移偏差问题

问题概述

采用有限差分法基于四阶控制方程求解超静定梁位移：

梁参数：长8m，x=0处铰支，x=4m、x=8m处辊支，承受10kN/m均布荷载
遇到的问题：
- 9节点离散时计算结果接近专业软件（最大挠度0.3201mm）
- 31节点提升精度时，初始计算位移偏差极大（最大挠度40.99mm）
- 修正右端项后不同节点数结果一致，但最大挠度仅0.1934mm，仅为专业软件结果（0.347mm）的约一半

原始Matlab代码

E = 2.1E+08
I = 18890/100^4
EI = E*I
L = 8;
x = linspace(0,L,9)
h = L/9
q = 10

% Attempt with 9 points

Matrix = [6 -4  1   0   0   0   0   0   0
        -4  6   -4  1   0   0   0   0   0
         1  -4  6   -4  1   0   0   0   0
         0  1   -4  6   -4  1   0   0   0
         0  0   1   -4  6   -4  1   0   0
         0  0   0   1   -4  6   -4  1   0
         0  0   0   0   1   -4  6   -4  1
         0  0   0   0   0   1   -4  6   1
         0  0   0   0   0   0   1   -4  6]

% Remove rows where displacement is zero
Matrix(:,[1,9,5]) = []
Matrix([1,9,5],:) = []

% Create the right-hand side

RHS = [q; 
       q; 
       q-1/h^4; 
       q + 4/h^4; 
       q-6/h^4; 
       q+4/h^4; 
       q-1/h^4; 
       q; 
       q]

RHS([1,9,5]) = []

% Find displacements in mm

y = ((inv(Matrix)*RHS)/EI)*1000

% Displacements with values from supports
y = [0;y(1:3);0;y(4:6);0]

% Attempt with 31 points
N = 31;
L = 8
x = linspace(0,L,N)
h = L/N
q = 10

% Construct matrix
mat1 = diag(ones(1,N)*6) ;                              
mat2 = diag(ones(1,N-1)*-4, 1);
mat3 = diag(ones(1,N-1)*-4, -1);
mat4 = diag(ones(1,N-2)*1, 2);
mat5 = diag(ones(1,N-2)*1, -2);

matrix = [mat1 + mat2 + mat3 + mat4 + mat5]

% Create right-hand side
rhs = (ones(1,N)*q)'
rhs((N+1)/2) = q - 6/h^4
rhs((N+1)/2 + 1) = q + 4/h^4
rhs((N+1)/2 - 1) = q + 4/h^4
rhs((N+1)/2 + 2) = q - 1/h^4
rhs((N+1)/2 - 2) = q - 1/h^4

% Remove rows where displacement is zero

matrix(:,[1,N,(N+1)/2]) = [];
matrix([1,N,(N+1)/2],:) = [];
rhs([1,N,(N+1)/2]) = []

% Find displacement in mm
y = ((inv(matrix)*rhs)/EI)*1000

修正后的右端项代码

Lambda = EI/h^4
rhs = (ones(1,N)*(q)/Lambda)'

修正后位移计算结果

问题根源

步长计算错误
linspace(0,L,N)生成N个节点，节点间的实际步长应为h=L/(N-1)，而非L/N。错误的步长会导致中间支座（x=4m）的节点位置偏移，直接破坏结构受力的对称性与计算精度。
有限差分矩阵与方程量纲不匹配
梁的四阶控制方程为：
$$EI \frac{d^4 y}{dx^4} = q$$
中心差分展开后应整理为：
$$\frac{EI}{h^4} \cdot [\text{系数矩阵}] \cdot {y} = {q}$$
现有代码仅构造纯系数矩阵，未引入EI/h^4因子，后续修正时仅缩放右端项，导致整体方程系数失真，位移结果被错误缩小。
边界条件处理错误
直接删除位移为0的节点对应的行和列，破坏了有限差分格式的连续性，尤其是中间支座附近的差分方程，原本需要依赖相邻节点位移构建，删除后导致方程体系失真。
右端项的无意义修正
均布荷载下所有节点的右端项应为q，对中间支座附近节点添加-6/h^4等修正项属于对差分格式的误解，直接导致初始计算结果严重偏差。

解决方案

1. 修正步长计算

N = 31;
L = 8;
x = linspace(0, L, N);
h = L/(N-1); % 正确步长：N个节点对应N-1个间隔

2. 正确构造方程体系

采用纯系数矩阵+修正右端项的形式，避免系数过大：

E = 2.1E+8; % kPa，对应210GPa
I = 18890 / 100^4; % m⁴（18890 mm⁴转换为m⁴）
EI = E * I; % kN·m²
q = 10; % kN/m

% 构造纯系数矩阵
mat1 = diag(ones(1,N)*6);                              
mat2 = diag(ones(1,N-1)*-4, 1);
mat3 = diag(ones(1,N-1)*-4, -1);
mat4 = diag(ones(1,N-2)*1, 2);
mat5 = diag(ones(1,N-2)*1, -2);
matrix = mat1 + mat2 + mat3 + mat4 + mat5;

% 构造右端项：对应控制方程整理后的形式
rhs = ones(N,1) * q * h^4 / EI;

3. 规范边界条件处理

采用划行划列法强制位移为0，保留矩阵维度与差分连续性：

% 标记位移为0的节点索引：x=0(第1节点)、x=4m(第16节点)、x=8m(第31节点)
fixed_nodes = [1, (N+1)/2, N];
for idx = fixed_nodes
    matrix(idx,:) = 0;
    matrix(idx,idx) = 1; % 强制该节点位移为0
    rhs(idx) = 0;
end

% 求解位移（用反斜杠代替inv，数值稳定性更好）
y = matrix \ rhs;
% 转换为mm
y = y * 1000;

4. 修正边界差分格式

对于两端铰支/辊支节点，利用弯矩为0（二阶导数为0）的条件修正边界行：

% x=0处（第1节点）：铰支y=0，弯矩EI*y''=0，差分格式为y0-2y1+y2=0
matrix(1,:) = 0;
matrix(1,1) = 1;
matrix(1,2) = -2;
matrix(1,3) = 1;
rhs(1) = 0;

% x=8处（第31节点）：辊支y=0，弯矩EI*y''=0，差分格式为y29-2y30+y31=0
matrix(N,:) = 0;
matrix(N,N) = 1;
matrix(N,N-1) = -2;
matrix(N,N-2) = 1;
rhs(N) = 0;

验证结果

修正上述错误后，31节点计算的最大挠度会收敛到专业软件的0.347mm左右，节点数越多，结果越接近理论值。

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

超静定梁有限差分法求解异常：节点数增加后结果偏差

问题概述 #

原始Matlab代码 #

修正后的右端项代码 #

修正后位移计算结果 #

问题根源 #

解决方案 #

1. 修正步长计算 #

2. 正确构造方程体系 #

3. 规范边界条件处理 #

4. 修正边界差分格式 #

验证结果 #