用MATLAB函数解方程

一、求解非特征值问题

工具箱用solvepde函数求解偏微分方程非特征值问题。该函数的语法格式为：

• result = solvepde(model) 返回model表示的偏微分方程的解。对于平稳偏微分方程，有model.IsTimeDependent = false，即model.EquationCoefficients中的时间导数系数m和d的值必须为0。
• result = solvepde(model,tlist) 返回tlist表示的各时间点处的解。要求model.EquationCoefficients中的时间导数系数m和d至少有一个非零。tlist必须是一个单调递增或递减的实数向量，如1:10。

返回值result为StationaryResults对象或TimeDependentResults对象。

StationaryResults对象是与时间无关的对象，它有Mesh, NodalSolution, XGradients, YGradients和ZGradients等5个属性，分别表示有限元网格、节点处的解、以及节点处梯度的x分量、y分量和z分量。

TimeDependentResults对象则是与时间有关的对象，除了有5个跟StationaryResults对象相同的属性，它还比StationaryResults对象多一个SolutionTimes属性，该属性的值为各时间点处的解组成的向量。向量的大小跟tlist参数的大小相同。

【例19】下面首先用polyshape函数创建PDE模型的几何模型。

>> t = pi/12:pi/12:2*pi;
>> pg = polyshape({[-0.9638 0.8717 0.8717 -0.9638], 0.4424*cos(t)-0.0082}, ...
                {[0.5312 0.5312 -0.4589 -0.4589], 0.1859*sin(t)+0.0296})
>> tr = triangulation(pg)
>> model = createpde;
>> tns = tr.Points';
>> tes = tr.ConnectivityList';
>> geometryFromMesh(model,tns,tes);

用applyBoundaryCondition函数设置边界条件。

>> applyBoundaryCondition(model,'dirichlet', ...
                             'Edge',1:6,'u',0);

用setInitialConditions函数设置初始条件。

>> setInitialConditions(model,0);

用specifyCoefficients函数设置模型的系数。

>> specifyCoefficients(model,'m',0,...
                          'd',0,...
                          'c',1,...
                          'a',0,...
                          'f',1);

用generateMesh函数生成几何模型的网格。

>> m=generateMesh(model,’Hmax’,0.06);

用solvepde函数进行求解，计算结果返回到results对象中。

>> results = solvepde(model)
results =
  StationaryResults - 属性:
    NodalSolution: [2108×1 double]
       XGradients: [2108×1 double]
       YGradients: [2108×1 double]
       ZGradients: []
             Mesh: [1×1 FEMesh]

results是一个StationaryResults对象，引用其属性可以获得网格节点处的解、X方向、Y方向和Z方向的梯度，以及网格对象。

二、求解特征值问题

工具箱用solvepdeeig函数求解偏微分方程特征值问题。该函数的语法格式为：

• result = solvepdeeig(model,evr) 求解model表示的偏微分方程特征值问题，要求特征值落在evr指定的范围内。如果该范围内不包含任何特征值，solvepdeeig函数返回一个EigenResults对象，对象的EigenVectors, EigenValues和Mesh属性的值都为空。

对方程的解进行可视化

常用图形来表现偏微分方程的解，这些图形包括色谱图、等值线图、向量场图、网格图和曲面图等。工具箱提供了pdeplot函数和pdeplot3D函数绘制二维问题和三维问题的解。

一、绘制二维问题的解

工具箱使用pdeplot函数绘制二维问题的解。该函数的语法格式为：

• pdeplot(model,"XYData",results.NodalSolution) 用彩色曲面图绘制节点处模型的解。默认时使用"jet"颜色查找表进行着色。
• pdeplot(model,"XYData",results.Temperature,"ColorMap","hot") 对于二维热力学模型，绘制节点处的温度。使用"hot"颜色查找表定义的颜色创建彩色曲面图。
• pdeplot(model,"XYData",results.VonMisesStress,"Deformation",results.Displacement) 对于二维结构力学模型，绘制米色斯应力和结构变形图。
• pdeplot(model,"XYData",results.ModeShapes.ux) 对于二维结构模态分析模型，绘制模态位移的x分量。
• pdeplot(model,"XYData",results.ElectricPotential) 对于二维静电学分析模型，绘制节点处的电势。
• pdeplot(model) 绘制模型的网格。
• pdeplot(mesh) 绘制二维模型对象的Mesh属性定义的网格。
• pdeplot(nodes,elements) 绘制节点和单元定义的网格。
• pdeplot(p,e,t) 绘制p, e和t数据定义的网格。
• pdeplot(___,Name,Value) 根据Name和Value匹配对指定的设置进行绘图。Name和Value匹配对的设置如表1-5所示。
• h = pdeplot(___) 返回图形句柄。

表1-5 Name和Value匹配对设置

pdeplot函数可以在pdetool GUI中和命令行使用。它可以同时显示3种图形。XYData可以用曲面图显示。曲面图可以使用刻面着色和插值着色两种方式，默认时使用插值着色。在曲面图上可以叠加等值线图（用黑色表示）或通过设置Contour为on来独立绘制等值线图（用彩色表示）。ZData可以用高度图表示。流数据可以用箭头图表示。所有类型的数据都可以是节点数据或三角形数据（流数据只能是三角形数据）。节点数据用长度为size(p,2)的列矢量表示，三角形数据用长度为size(t,2)的行矢量表示。若没有提供XYData、ZData或者FlowData，则pdeplot函数根据p，e和t绘制网格图。

Mesh选项控制图中三角形网格的显示或隐藏（默认选项）。选项XYGrid首先将数据转换为xy数据（使用tri2grid函数），然后使用标准MATLAB绘图算法。GridParam 属性将tri2grid数据传到pdeplot函数中，这会加速动画显示。ColorBar属性控制图中色条的添加。输入Title 参数，可以给当前图形输入标题。Levels参数仅适用于等值线图。

【例20】1.5.6小节结合一个实例进行了模型求解，下面在该例的基础上对问题的解进行可视化。利用results对象的NodalSolution属性可以获取网格各节点处的解，然后使用pdeplot函数利用解数据绘色谱图。

>> u = results.NodalSolution;
>> pdeplot(model,'XYData',u)
>> axis equal

生成图5-7。

图5-7 利用解数据绘制色谱图

下面利用解数据绘制三维曲面图。

>> pdeplot(model,"XYData",u,"ZData",u)

生成图5-8。

图5-8 利用解数据绘制三维曲面图

下面仍然用解数据绘图，但是将XYStyle参数的值设置为Flat，此时网格中每个单元用同一种颜色进行绘制。这样生成的图称为刻面图。

>> pdeplot(model,"XYData",u,"ZData",u,"XYStyle","Flat")

生成图5-9。

图5-9 利用解数据绘制三维刻面图

下面利用results对象的XGradients属性和YGradients属性获取网格节点处X方向和Y方向的梯度，然后利用梯度数据绘制向量场图。

>> ux = results.Xgradients;
>> uy = results.Ygradients;
>> pdeplot(model,”FlowData”,[ux,uy])
>> axis equal

生成图5-10。

图5-10 利用梯度数据绘制向量场图

下面利用解数据绘图时，设置Contour参数的值为on，在色谱图上叠加等值线图。

>> pdeplot(model,'XYData',u,'Contour','on')
>> axis equal

生成图5-11。

图5-11 在色谱图上叠加等值线图

下面利用解数据绘制曲面图，同时设置了多个参数的值。用ColorMap参数设置颜色查找表的名字为hot，用FaceAlpha参数设置面的透明度为0.6，用FlowData参数叠加向量场图，用Contour参数叠加等值线图，用Mesh参数叠加网格线。

>> pdeplot(model,'XYData',u,'ZData',u,...
             'ColorMap','hot',...
             'FaceAlpha',0.6,...
             'FlowData',[ux,uy],...
             'Contour','on',...
             'Mesh','on')

生成图5-12。

图5-12 在曲面图上叠加多个参数的设置效果

二、绘制三维问题的解

工具箱使用pdeplot3D函数绘制三维问题的解。该函数的语法格式为：

• pdeplot3D(model,"ColorMapData",results.NodalSolution) 绘制三维模型表面节点处解的色谱图。
• pdeplot3D(model,"ColorMapData",results.Temperature) 对于三维热力学分析模型，绘制节点处的温度。
• pdeplot3D(model,"ColorMapData",results.VonMisesStress,"Deformation",results.Displacement) 对于三维结构力学分析模型，绘制米色斯应力和网格变形图。
• pdeplot3D(model,"ColorMapData",results.ElectricPotential) 对于三维静电学分析模型，绘制节点处的电势。
• pdeplot3D(model) 绘制模型的曲面网格。
• pdeplot3D(mesh) 用三维模型对象的Mesh属性定义的数据绘制网格。
• pdeplot3D(nodes,elements) 绘制节点和单元定义的网格。
• pdeplot3D(___,Name,Value) 绘制Name和Value匹配对设置的曲面网格、节点处的数据，或二者兼而有之。Name和Value匹配对的设置请参见pdeplot函数的对应设置。
• h = pdeplot3D(___) 返回图形的句柄。

【例21】下面创建一个二维几何模型，然后将它垂向拉伸1个单位长度生成三维几何模型。

>> t = pi/12:pi/12:2*pi;
>> pg = polyshape({[-0.9638 0.8717 0.8717 -0.9638], 0.4424*cos(t)-0.0082}, ...
{[0.5312 0.5312 -0.4589 -0.4589], 0.1859*sin(t)+0.0296});
>> tr = triangulation(pg);
>> model = createpde;
>> tns = tr.Points';
>> tes = tr.ConnectivityList';
>> [g,mesh]=geometryFromMesh(model,tns,tes);
>> extrude(g,1);    %垂向拉伸1个单位长度
>> pdegplot(g,'FaceLabels','on','FaceAlpha',0.5)
>> axis equal

生成的三维几何模型如图5-13所示。

图5-13 创建三维几何模型

给模型添加边界条件并设置模型系数，生成网格，然后用solvepde函数进行求解。

>> applyBoundaryCondition(model,"dirichlet","Face",[1:6],"u",0);
>> specifyCoefficients(model,"m",0,"d",0,"c",1,"a",0,"f",1);
>> generateMesh(model);
>> results = solvepde(model)
results =
  StationaryResults - 属性:
    NodalSolution: [32738×1 double]
       XGradients: [32738×1 double]
       YGradients: [32738×1 double]
       ZGradients: [32738×1 double]
             Mesh: [1×1 FEMesh]

用pdeplot3D函数绘制三维模型的网格。

>> pdeplot3D(model,'Mesh','on')

生成图5-14。

图5-14 三维模型的网格

下面利用解数据会三维模型的色谱图，设置面的透明度为0.5。

>> pdeplot3D(model,"ColorMapData",u,'FaceAlpha',0.5)

生成图5-15。

图5-15 利用解数据绘图