MATLAB求解偏微分方程

发表于 | 分类于 |

1. 偏微分方程工具箱介绍

本工具箱采用有限元分析求解偏微分方程

​ 偏微分方程工具箱提供有限元分析在 2-D，3-D和时间变量的空间下求解偏微分方程。通过它你可以选定和网格化 2-D 和 3-D 的几何空间并且规定边界条件和方程。也就是说在几何上，你可以用它求解 定常域（不含时间）、时域、频域和特征值问题。用于后处理的函数，如绘制结果的函数可以可视化研究结果。

​ 它能用于求解从标准问题的提炼出的偏微分方程。例如，扩散、热传导、结构力学、静电学、电磁学 和 交流电电磁学 （原作: AC power electromagnetics）以及定制耦合PDEs系统的问题。

1.1主要特点

  • 可以求解耦合PDEs系统：定常，时域，频域和特征值
  • PDE要符合椭圆型、抛物线型和双曲线型问题
  • 边界条件的选择： Dirichlet（第一类边界条件），generalized Neumann（第二类边界条件)，and mix(Robin or 第三类边界条件)
  • 由STL文件输入，可用函数创建2-D几何图形和3-D几何图像
  • 要自动网格化，请使用 tetrahedra and triangles
  • 同时可视化多方案的参量、网格覆盖和动画处理

2. 例子

2.1 在方域上的非均匀热方程

这个方程是为了展示怎么用 solvepde函数求解热方程

​ 单位热源的热方程的一般形式是

2.1.1 定义问题

​ 创建一个$x \in\left[ {-1,1} \right] ,y\in\left[ {-1,1} \right]​$的正方形域，里面包含一个半径为0.4与正方形同心的圆作为选定的初始条件。

1
2
3
4
5
6
7
R1 = [3;4;-1;1;1;-1;-1;-1;1;1];
C1 = [1;0;0;0.4];
C1 = [C1;zeros(length(R1)-length(C1),1)];
gd = [R1,C1];
sf = 'R1+C1';
ns = char('R1','C1')';
g = decsg(gd,sf,ns);

2.1.2 创建一维偏微分方程

1
2
number0fPDE = 1;
pdem = createpde(number0fPDE);

2.1.3 创建一个几何体附加到PDE模型

1
geometryFormEdges(pdem,g)

2.1.4 应用边界条件

​ 画出几何图形并且用边界条件定义的显示出边缘的标识。

1
2
3
4
5
6
7
8
figure
pdegplot(pdem,'edgeLabels','on','subdomainLabels','on')  %
axis([-1.1 1.1 -1.1 1.1]);                   %
axis equal                                   %
title 'Geometry with Edge and Subdomain Labels'

% Solution is zero at all four outer edges of the square
applyBoundaryCondition(pdem,'Edge',(1:4),'u',0);      %

2.1.5 选定PDE的系数

1
2
3
4
5
specifyCoefficients(pdem,'m',0,...
                         'd',1,...
                         'c',0,...
                         'a',0,...
                         'f',1):

2.1.6 选定初始条件

​ 在圆内，不连续初始条件值为$1$，在圆外则为$0$.选定每处初始条件都为0。

1
2
3
setInitialConditions(pdem,0);
% 对于圆内选择非零初始条件，face ID 2
setInitialConditions(pdem,1,'face',2)

2.1.7 生成网格

1
2
3
4
msh = generateMesh(pdem);
figure;
pdemesh(pdem)
axis equal

2.1.8 对时间离散化

​ 这里我们在$t \in \left[ {0,0.1} \right]​$求20个点的值。

1
2
nframes = 20;
tlist = linspace(0,0.1,nframes);

2.1.9 用pdem的SolverOptions.Reportstartistics设置为on求解

1
2
3
pdem.SolveOptions.Reportstatistics = 'on';
result = solvepde(pdem,tlist);
u1 = result.NodalSolution;

2.1.10 绘画结果

1
2
3
4
5
6
7
8
9
10
figure
umax = max(max(u1));
umin = min(min(u1));
for j = i:nframes,
    pdeplot(pdem,'xydata',u1(:,j),'zdata',u1(:,j));
    caxis([umin,umax]);
    axis([-1 1 -1 1 0 1]);
    Mv(j) = getframe;
end
movie(Mv,1);

函数附录

本附录主要对一些函数进行解释，且按问题学习分布。例如，2.1 在方域上的非均匀热方程中的陌生函数将出现在A部分。详情请看matlab帮助文档。

A

1. decsg

将构建的实体几何体分解为小部分。

1.1 语法
1
2
3
4
5
6
dl = decsg(gd,sf,ns)
dl = decsg(gd)
[dl,bt] = decsg(gd)
[dl,bt] = decsg(gd,sf,ns)
[dl,bt,dl1,dt1,msb] = decsg(gd)
[dl,bt,dl1,dt1,msb] = decsg(gd,sf,ns)
1.2 描述