一维微分方程有限元刚度矩阵推导求助
嘿，我来帮你理清楚这个问题！首先你卡在分部积分那一步其实是没注意到原方程可以整理成更友好的散度形式，咱们一步步拆解：
第一步：简化弱形式的分部积分
你的原方程是：
$$-e^xu'' -e^xu' + u = f(x)$$
仔细看前两项，其实可以用乘积法则重新组合：
$$-e^xu'' -e^xu' = -\frac{d}{dx}\left(e^x u'\right)$$
因为对$e^x u'$求导的话，$\frac{d}{dx}(e^x u') = e^x u' + e^x u''$，刚好和前两项的相反数对应！
那原方程就可以写成：
$$-\frac{d}{dx}\left(e^x u'\right) + u = f(x)$$
现在再写弱形式就简单多了：对任意满足边界条件$\phi_i(0)=\phi_i(1)=0$的试函数$\phi_i$，有：
$$\int_0^1 -\frac{d}{dx}\left(e^x u'\right)\phi_i dx + \int_0^1 u\phi_i dx = \int_0^1 f\phi_i dx$$
对第一个积分做分部积分：
$$\int_0^1 -\frac{d}{dx}\left(e^x u'\right)\phi_i dx = -\left.e^x u'\phi_i\right|_0^1 + \int_0^1 e^x u'\phi_i' dx$$
因为$\phi_i$在端点0和1处为0，所以边界项直接消失，最终弱形式简化为：
$$\int_0^1 e^x u'\phi_i' dx + \int_0^1 u\phi_i dx = \int_0^1 f\phi_i dx$$
是不是一下子清爽多了？之前的交叉项直接抵消了，不用纠结多函数的分部积分~
第二步：手动构造N×N刚度矩阵
假设咱们用最常用的分段线性基函数，节点为$x_0=0, x_1, x_2, ..., x_N=1$（这里的N×N矩阵对应内部N个节点，端点0和1因齐次边界条件被排除）。每个基函数$\phi_j(x)$只在相邻两个节点的区间内非零，所以刚度矩阵的元素$K_{ij}$只有当$i$和$j$相邻（或相等）时非零：
1. 对角元（i=j）
$$K_{ii} = \int_{x_{i-1}}^{x_i} e^x (\phi_i')^2 dx + \int_{x_i}^{x_{i+1}} e^x (\phi_i')^2 dx + \int_{x_{i-1}}^{x_i} \phi_i^2 dx + \int_{x_i}^{x_{i+1}} \phi_i^2 dx$$
对于分段线性基函数：
- 在$[x_{i-1},x_i]$上，$\phi_i(x)=\frac{x-x_{i-1}}{h_i}$，$\phi_i'=\frac{1}{h_i}$（$h_i=x_i-x_{i-1}$）
- 在$[x_i,x_{i+1}]$上，$\phi_i(x)=\frac{x_{i+1}-x}{h_{i+1}}$，$\phi_i'=-\frac{1}{h_{i+1}}$（$h_{i+1}=x_{i+1}-x_i$）
代入计算后：
$$K_{ii} = \frac{e{x_i}-e{x_{i-1}}}{h_i^2} + \frac{e{x_{i+1}}-e{x_i}}{h_{i+1}^2} + \frac{h_i + h_{i+1}}{3}$$
（内部节点都有左右两个区间，端点处的对角元只需保留单侧区间的项）
2. 非对角元（|i-j|=1）
比如i=j+1（下对角），重叠区间是$[x_j,x_{j+1}]$，此时：
$$K_{i,j} = \int_{x_j}^{x_{j+1}} e^x \phi_i'\phi_j' dx + \int_{x_j}^{x_{j+1}} \phi_i\phi_j dx$$
代入基函数表达式计算得：
$$K_{i,j} = -\frac{e{x_{j+1}}-e{x_j}}{h_{j+1}^2} + \frac{h_{j+1}}{6}$$
因为矩阵是对称的，$K_{j,i}=K_{i,j}$，上对角元直接复用这个结果就行。
第三步：用程序自动生成矩阵
如果不想手动算积分，推荐用有限元库或者符号计算工具：
用FEniCS（DOLFINx）快速组装
FEniCS是专门的有限元库，能直接定义弱形式并自动组装矩阵，伪代码如下：
from dolfin import * # 创建单位区间网格，N是内部节点数 mesh = UnitIntervalMesh(N) # 定义分段线性函数空间 V = FunctionSpace(mesh, 'P', 1) # 定义齐次边界条件 def boundary(x, on_boundary): return on_boundary bc = DirichletBC(V, Constant(0), boundary) # 定义试函数和测试函数 u = TrialFunction(V) v = TestFunction(V) # 定义弱形式（和咱们推导的一致） a = exp(x)*dot(grad(u), grad(v))*dx + u*v*dx # 定义右端项，这里假设f是已知函数，比如f=1的话就写Constant(1) f = Expression("你的函数表达式", degree=1) L = f*v*dx # 组装刚度矩阵A和荷载向量b A = assemble(a) b = assemble(L) # 应用边界条件 bc.apply(A, b) # 现在A就是你要的N×N刚度矩阵
用SymPy做符号推导
如果想验证手动计算的结果，可以用SymPy做符号积分，比如对单个区间的积分进行符号计算，确保公式正确：
import sympy as sp x = sp.symbols('x') h = sp.symbols('h') x0 = sp.symbols('x0') # 定义基函数 phi_i = (x - x0)/h phi_j = (x0 + h - x)/h # 计算积分 int1 = sp.integrate(sp.exp(x)*sp.diff(phi_i, x)*sp.diff(phi_j, x), (x, x0, x0+h)) int2 = sp.integrate(phi_i*phi_j, (x, x0, x0+h)) print(int1 + int2)
这样就能快速得到符号形式的矩阵元素，和手动推导的结果对比验证。
备注：内容来源于stack exchange，提问作者Makogan
