3. 一维杆问题弱形式的程序实现

上一节我们得到了以试探函数和检验函数表述的弱形式，这一节把它真正落到代码上。一维杆问题的弱形式为：求 $u\in H^1$ 且满足本质边界条件，使对任意检验函数 $v$（在本质边界上为零）有

$$\int_0^L EA\,\frac{du}{dx}\frac{dv}{dx}\,dx=\int_0^L f v\,dx+\left[EA\frac{du}{dx}v\right]_0^L$$

程序实现的核心任务，就是把连续形式转成线性方程组 $\mathbf{K}\mathbf{u}=\mathbf{F}$。整个过程可拆成四件事：Galerkin 离散、单元矩阵与载荷向量推导、全局组装、边界条件处理。下图给出这四步之间的数据流与关键量。

Galerkin 离散：把积分变成求和

Galerkin 法的要点是让检验函数与试探函数取自同一组基。把杆划分为 $n$ 个节点、$n-1$ 个单元，取全局形函数 $\varphi_i(x)$（在节点 $i$ 处为 1、其余节点为 0 的分片线性函数），令

$$u_h=\sum_{j=1}^{n}u_j\varphi_j(x),\qquad v=\varphi_i(x)$$

代入弱形式，得到第 $i$ 个方程：$\sum_j u_j\int EA\,\varphi_i'\varphi_j'\,dx=\int f\varphi_i\,dx+\text{边界项}$。写成矩阵即 $\mathbf{K}\mathbf{u}=\mathbf{F}$，其中 $K_{ij}=\int EA\,\varphi_i'\varphi_j'\,dx$，$F_i=\int f\varphi_i\,dx$。由于 $\varphi_i$ 只在相邻单元非零，$\mathbf{K}$ 是稀疏三对角矩阵——这正是程序里按单元循环组装、而不是逐元素积分的原因。

单元刚度矩阵与载荷向量推导

在单元 $e$ 上引入局部坐标 $\xi\in[0,1]$，$x=x_e+\xi L^e$，形函数 $N_1=1-\xi$、$N_2=\xi$，导数 $dN_1/dx=-1/L^e$、$dN_2/dx=1/L^e$。于是

$$K^e_{ab}=\int_0^{L^e}EA\,\frac{dN_a}{dx}\frac{dN_b}{dx}\,dx=\frac{EA}{L^e}\begin{bmatrix}1&-1\\-1&1\end{bmatrix}$$

对均布轴力 $f$（单位长度载荷），载荷向量为 $f^e_a=\int_0^{L^e}fN_a\,dx=\frac{fL^e}{2}\begin{bmatrix}1\\1\end{bmatrix}$，即把总轴力平均分给两个节点。若载荷为集中力 $P$ 作用于某节点，则直接累加到该节点的 $F_i$ 上，无需积分。这一步的推导结果直接决定了代码里两个常数矩阵的写法。

组装过程

组装是把每个单元的 $2\times2$ 矩阵按节点编号“对号入座”叠加到全局矩阵。设单元 $e$ 连接节点 $e$ 与 $e+1$，则其四个元素分别累加到 $K(e,e)$、$K(e,e+1)$、$K(e+1,e)$、$K(e+1,e+1)$。用伪代码表示：

K = zeros(n, n);  F = zeros(n, 1)
for e in 1..n-1:
    Le = x[e+1] - x[e]
    ke = EA/Le * [[1, -1], [-1, 1]]
    fe = f*Le/2 * [1, 1]
    i, j = e, e+1
    K[i,i] += ke[0][0];  K[i,j] += ke[0][1]
    K[j,i] += ke[1][0];  K[j,j] += ke[1][1]
    F[i] += fe[0];       F[j] += fe[1]

注意相邻单元在共享节点处会叠加，这正是刚度矩阵对角元随单元数增加而增大的来源。组装完成后，$\mathbf{K}$ 对称、半正定，尚未施加边界条件时是奇异的（存在刚体平移），必须处理边界才能求解。

本质边界条件处理

设节点 1 固定，$u_1=0$；若为指定位移 $u_1=\bar u$，常用“置零行列法”：把第 1 行、第 1 列除对角元外全部置零，对角元置 1，右端项第 1 个分量置 $\bar u$；同时对所有其他行 $i$，执行 $F_i\leftarrow F_i-K_{i1}\bar u$，再把 $K_{i1}$ 置零。这样既保留了矩阵对称性，又把已知位移的影响正确移到右端。若 $\bar u=0$，则只需置零行列并对角置 1、右端置 0，其余行无需修正。

程序编写与验证

下面给出一个可直接运行的最小实现（Python，仅用列表，便于对照公式）：

def solve_bar(x, EA, f, fixed, ubar):
    n = len(x)
    K = [[0.0]*n for _ in range(n)]
    F = [0.0]*n
    for e in range(n-1):
        Le = x[e+1] - x[e]
        k = EA/Le
        K[e][e]   += k;  K[e][e+1] -= k
        K[e+1][e] -= k;  K[e+1][e+1] += k
        F[e]   += f*Le/2
        F[e+1] += f*Le/2
    # 置零行列法施加 u[fixed] = ubar
    for i in range(n):
        if i != fixed:
            F[i] -= K[i][fixed]*ubar
            K[i][fixed] = 0.0
        K[fixed][i] = 0.0
    K[fixed][fixed] = 1.0
    F[fixed] = ubar
    # 高斯消元（略）
    return gauss(K, F)

验证算例：取长 $L=2\,\text{m}$、$EA=200\,\text{kN}$、均布轴力 $f=10\,\text{kN/m}$、左端固定、右端自由。解析解为 $u(x)=\dfrac{f}{2EA}x(2L-x)$，右端位移 $u(L)=\dfrac{fL^2}{2EA}=\dfrac{10\times4}{2\times200}=0.1\,\text{m}$。用 4 个等长单元（$L^e=0.5$）计算，右端节点位移应精确等于 0.1 m——因为分片线性形函数对二次位移场虽不能精确表示，但对本问题（常 $EA$、均布 $f$）的节点位移恰好精确。若把单元数改为 2 或 8，右端位移仍为 0.1 m，可作为程序正确性的判据。

进一步验证可检查：① 反力 $R=K_{11}u_1+\dots$ 应等于 $-fL=-20\,\text{kN}$（与外力平衡）；② 单元应力 $\sigma^e=EA(u_{e+1}-u_e)/L^e$ 沿杆线性变化，左端最大、右端为零；③ 将 $f$ 改为集中力 $P=20\,\text{kN}$ 作用于右端，右端位移应为 $PL/EA=0.2\,\text{m}$。三项检查全部通过，说明离散、组装、边界处理与求解链路均正确。

至此，一维杆弱形式已完整转化为可执行程序。这套“单元矩阵—组装—边界处理—求解—验证”的骨架，在二维、三维问题中只是把标量 $K^e$ 换成矩阵、把节点编号换成拓扑连接表，思想完全一致。