有限元法

一、直观理解：分片简单函数拼出整体解

有限元法（Finite Element Method, FEM）把求解区域切成许多简单小单元（线段、三角形、四面体等），在每个单元上用多项式等“形状函数”近似未知场，再通过拼装得到整体近似解。

对比有限差分：

差分强调网格点上的差商；
有限元强调分片函数空间中的最佳逼近（在弱形式意义下）。

可以把温度场想象成一块不规则铁片上的连续曲面。有限元不是强迫曲面只在矩形格点上取值，而是用许多小三角形“帐篷函数”叠加，去贴合真实曲面。区域边界弯曲时，三角形网格比矩形差分网格自然得多。

有限元在结构力学、椭圆型方程、复杂几何热传导与电磁中极为常用。MATLAB PDE Toolbox 的 createpde、solvepde 即基于有限元思想。

二、适用场景与优缺点

类型	说明
适用场景	复杂几何、局部加密需求强；变系数、各向异性介质；以稳态椭圆或抛物问题为主；需要弱形式自然施加 Neumann/Robin 边界；工程场问题与建模中的二维区域方程。
核心优点	几何适应性强；数学框架统一（弱形式 + 有限维子空间）；高阶元可系统提升精度；局部加密（自适应）成熟。
局限与缺点	理论与实现门槛高于差分；网格质量影响很大；对流占优问题需稳定化；守恒形式与激波问题不如有限体积直观；工具箱依赖可能掩盖方程设置错误。

使用前应确认：

能写出或确认 PDE 的强形式与边界类型；
理解 Dirichlet / Neumann / Robin 的物理含义；
网格尺度与感兴趣梯度区域匹配；
有简化几何上的解析解或实验数据做校验。

三、核心原理：弱形式、基函数与刚度矩阵

以二维 Poisson 方程为例：

-\nabla\cdot(c\nabla u)+a u=f\quad\text{in }\Omega,

边界 $\partial\Omega$ 上给定 $u=g_D$ 或 $c\nabla u\cdot n=g_N$ 等。

1. 弱形式

取满足齐次 Dirichlet 条件的试验函数 $v$ ，分部积分得

\int_\Omega c\nabla u\cdot\nabla v\,dx +\int_\Omega a u v\,dx =\int_\Omega f v\,dx +\text{边界项}.

弱形式把二阶导数降为一阶，对近似空间更友好。

2. Galerkin 离散

选有限维空间 $V_h=\mathrm{span}\{\phi_1,\ldots,\phi_N\}$ ，令

u_h=\sum_{j=1}^N U_j\phi_j,

代入弱形式并取 $v=\phi_i$ ，得到线性系统

KU=F,

其中刚度矩阵

K_{ij}=\int_\Omega c\nabla\phi_j\cdot\nabla\phi_i\,dx+\int_\Omega a\phi_j\phi_i\,dx,

载荷

F_i=\int_\Omega f\phi_i\,dx+\text{边界贡献}.

3. 单元装配

实际计算在单元上积分，再按自由度编号累加到全局 $K,F$ 。线性三角元是最常见入门选择。

4. 误差直观

在正则问题与合适范数下，线性元常有 $\|u-u_h\|_{H^1}=O(h)$ ， $\|u-u_h\|_{L^2}=O(h^2)$ （ $h$ 为网格尺寸）。网格越细、单元阶越高，通常越准——前提是解本身足够光滑且网格不畸形。

四、关键决策

方程系数与 toolbox 约定
MATLAB 中 specifyCoefficients 的 m,d,c,a,f 对应标准型，必须与教材方程逐项对齐。
边界条件
绝热、固定温度、对流换热对应不同 BC；写错边界比网格粗往往更致命。
网格策略
几何拐角、材料界面、载荷集中处加密；避免极扁三角形。
验证
- 网格加密后关键输出（最大温度、通量）收敛；
- 与一维简化或有限差分对照；
- 检查对称性、极值原理（无源时极值在边界等）。
瞬态问题
空间 FEM + 时间差分/ode15s 半离散是常见路线；时间步与空间误差应匹配。

五、MATLAB 示例

%% Poisson：-Δu = 1，圆域，边界 u=0
model = createpde();
geometryFromEdges(model, @circleg);
applyBoundaryCondition(model, "dirichlet", ...
    "Edge", 1:model.Geometry.NumEdges, "u", 0);
specifyCoefficients(model, "m", 0, "d", 0, "c", 1, "a", 0, "f", 1);
generateMesh(model, "Hmax", 0.15);
result = solvepde(model);

figure;
pdeplot(model, "XYData", result.NodalSolution, "ZData", result.NodalSolution);
title("圆域 Poisson 有限元解");
colormap turbo; colorbar;

%% 网格加密收敛：监测最大解
HmaxList = [0.3, 0.2, 0.1, 0.05];
umax = zeros(size(HmaxList));
ndof = zeros(size(HmaxList));
for i = 1:numel(HmaxList)
    generateMesh(model, "Hmax", HmaxList(i));
    R = solvepde(model);
    umax(i) = max(R.NodalSolution);
    ndof(i) = numel(R.NodalSolution);
end
disp(table(HmaxList(:), ndof(:), umax(:), ...
    'VariableNames', {'Hmax', 'nDoF', 'u_max'}));

%% 矩形导热稳态：-div(k grad T)=Q，三边绝热一边定温
% 使用 decsg 构造矩形
model2 = createpde();
R1 = [3,4, 0,1,1,0, 0,0,0.5,0.5]';  % 矩形 [0,1]x[0,0.5]
g = decsg(R1, 'R1', ('R1')');
geometryFromEdges(model2, g);
% 边编号需用 pdegplot(model2,'EdgeLabels','on') 查看后设置
pdegplot(model2, "EdgeLabels", "on");
% 假设 Edge 1 为底边定温，其余 Neumann=0（默认）
applyBoundaryCondition(model2, "dirichlet", "Edge", 1, "u", 0);
% 其余边不设 Dirichlet 即齐次 Neumann（绝热）
k = 1; Q = 2;
specifyCoefficients(model2, "m", 0, "d", 0, "c", k, "a", 0, "f", Q);
generateMesh(model2, "Hmax", 0.05);
R2 = solvepde(model2);
figure;
pdeplot(model2, "XYData", R2.NodalSolution, "Contour", "on");
title("矩形区域稳态热传导（示意）"); colorbar;

%% 一维“有限元味道”的手写对照：-u''=pi^2 sin(pi x), u(0)=u(1)=0
% 线性元等价于某些情形下的特殊差分；用于理解 KU=F
n = 20;               % 单元数
h = 1/n;
x = linspace(0,1,n+1);
% 组装内点 1..n-1
K = zeros(n-1); F = zeros(n-1,1);
for e = 1:n
    % 单元 [x_{e-1}, x_e]，局部基
    nodes = [e-1, e];  % 0-based 全局节点号
    Ke = (1/h) * [1 -1; -1 1];
    % 载荷用节点插值近似
    fe = [pi^2*sin(pi*x(e)); pi^2*sin(pi*x(e+1))];
    Fe = h/2 * fe;     % 简单梯形/节点规则示意
    for a = 1:2
        A = nodes(a);
        if A==0 || A==n, continue; end  % Dirichlet 固定
        F(A) = F(A) + Fe(a);
        for b = 1:2
            B = nodes(b);
            if B==0 || B==n, continue; end
            K(A,B) = K(A,B) + Ke(a,b);
        end
    end
end
U = K\F;
u = [0; U; 0];
uTrue = sin(pi*x)';
figure;
plot(x, u, "o-", x, uTrue, "-"); grid on;
legend("1D 线性元组装", "解析解");
title(sprintf("max error = %.3e", max(abs(u-uTrue))));

六、如何解释结果

有限元结果解释应覆盖：

强形式 PDE、系数含义与单位；
弱形式或 toolbox 系数如何对应；
几何、边界条件、网格规模（节点/单元数、Hmax）；
场分布图与关键标量（最大/最小、某点值、总通量）；
网格加密表，证明量已进入收敛区；
物理解读是否符合极值原理、对称性与守恒直觉；
不确定性：材料参数、边界理想化带来的模型误差（往往大于离散误差）。

避免只贴彩色云图而不给网格与边界说明——那会使结果不可复核。

七、常见误区

1. 边界条件与物理场景错位

把对流换热写成定温、或漏设边界，都会得到“漂亮但错误”的场。

2. 网格极粗却讨论细结构

若关心局部热点，必须局部加密并做收敛，不能凭默认网格下结论。

3. 把 toolbox 默认系数模型当成自己的方程

c,a,f,m,d 填错一个，方程就变了。应回到文档中的标准形式逐项核对。

4. 畸形单元

大长宽比三角形导致条件数恶化与局部误差大，需要检查网格质量。

5. 对流占优仍用朴素 Galerkin

可能出现非物理振荡，需迎风/稳定化或改用更合适的方法。

八、小结

有限元法的核心思想可以概括为：

将 PDE 转化为弱形式，在分片多项式空间中求 Galerkin 近似，通过单元积分与全局装配得到稀疏代数系统。

它特别擅长复杂几何与工程场问题。建模中使用 FEM 时，几何与边界的正确性优先于“单元是否足够高级”；务必用网格收敛和简化情形校验，再采信定量场分布与导出指标。

一、直观理解：分片简单函数拼出整体解​

二、适用场景与优缺点​

三、核心原理：弱形式、基函数与刚度矩阵​

1. 弱形式​

2. Galerkin 离散​

3. 单元装配​

4. 误差直观​

四、关键决策​

五、MATLAB 示例​

六、如何解释结果​

七、常见误区​

1. 边界条件与物理场景错位​

2. 网格极粗却讨论细结构​

3. 把 toolbox 默认系数模型当成自己的方程​

4. 畸形单元​

5. 对流占优仍用朴素 Galerkin​

八、小结​