MATLAB中借助pdetool 实现有限元求解Possion方程

二维泊松方程求解

问题描述

考虑二维矩形区域 $\Omega =[0,1]\times [0,1]$，求解泊松方程：

$$-{\nabla }^{2}u=f(x,y)\text{in}\Omega$$

边界条件：

$$u=0\text{on}\partial \Omega (\text{Dirichlet边界})$$

假设源项 $f(x,y)=2{\pi }^2\sin (\pi x)\sin (\pi y)$，其解析解为 $u(x,y)=\sin (\pi x)\sin (\pi y)$。

解析解验证

将解析解 $u(x,y)=\sin (\pi x)\sin (\pi y)$ 代入泊松方程：

$$\begin{array}{rl}\frac{{\partial }^{2}u}{\partial x^2}& =-{\pi }^2\sin (\pi x)\sin (\pi y)\\ \frac{{\partial }^{2}u}{\partial y^2}& =-{\pi }^2\sin (\pi x)\sin (\pi y)\\ {\nabla }^{2}u& =-2{\pi }^2\sin (\pi x)\sin (\pi y)\\ -{\nabla }^{2}u& =2{\pi }^2\sin (\pi x)\sin (\pi y)=f(x,y)\end{array}$$

验证了该解析解的正确性。

边界条件验证

在边界 $\partial \Omega$ 上：

• 当 $x=0$ 或 $x=1$ 时，$\sin (\pi x)=0$
• 当 $y=0$ 或 $y=1$ 时，$\sin (\pi y)=0$

因此 $u(x,y)=\sin (\pi x)\sin (\pi y)=0$，满足Dirichlet边界条件。
使用 MATLAB 的 generateMesh 函数生成三角形网格（需安装 PDE Toolbox）：

刚度矩阵组装

% 创建几何模型并生成网格
model = createpde();
geometryFromEdges(model, @squareg);
mesh = generateMesh(model, 'Hmax', 0.1);
p = mesh.Nodes;
t = mesh.Elements;% 绘制网格
figure;
pdeplot(model);
title('有限元网格');% 初始化参数
nNodes = size(p, 2);
nElements = size(t, 2);
K = sparse(nNodes, nNodes);
F = zeros(nNodes, 1);% 主循环 - 组装全局矩阵
for e = 1:nElementsnodes = t(1:3, e);coords = p(:, nodes);% 计算单元矩阵[Ke, fe] = computeElementMatrices(coords);% 组装到全局矩阵K(nodes, nodes) = K(nodes, nodes) + Ke;F(nodes) = F(nodes) + fe;
end% 施加边界条件
tol = 1e-6;
boundaryNodes = find(p(1,:) < tol | p(1,:) > 1-tol | p(2,:) < tol | p(2,:) > 1-tol);% 修正刚度矩阵和载荷向量（置一法）
K(boundaryNodes, :) = 0;
K(:, boundaryNodes) = 0;
K(boundaryNodes, boundaryNodes) = speye(length(boundaryNodes));
F(boundaryNodes) = 0;% 求解线性系统
u = K \ F;% 可视化有限元解
figure;
pdeplot(p, [], t, 'XYData', u, 'Contour', 'on', 'ColorMap', 'jet');
title('有限元解');
xlabel('x'); ylabel('z');% 计算并绘制解析解
[X, Y] = meshgrid(0:0.05:1, 0:0.05:1);
U_analytic = sin(pi*X) .* sin(pi*Y);
figure;
surf(X, Y, U_analytic);
title('解析解');
xlabel('x'); ylabel('y'); zlabel('z');% 单元矩阵计算函数
function [Ke, fe] = computeElementMatrices(coords)% 计算三角形面积x = coords(1, :);y = coords(2, :);A = 0.5 * abs((x(2)-x(1))*(y(3)-y(1)) - (x(3)-x(1))*(y(2)-y(1)));% 计算B矩阵（应变-位移矩阵）B = [y(2)-y(3), y(3)-y(1), y(1)-y(2);x(3)-x(2), x(1)-x(3), x(2)-x(1)] / (2*A);% 计算单元刚度矩阵Ke = A * (B' * B);% 计算单元载荷向量（假设f=2π²sin(πx)sin(πy)）xc = mean(x);yc = mean(y);f_val = 2 * pi^2 * sin(pi*xc) * sin(pi*yc);fe = (A/3) * f_val * ones(3, 1); % 平均分配到三个节点
end