• 2021-04-18 12:38:59

function [u,x,y] = poisson(f,g,bx0,bxf,by0,byf,D,Mx,My,tol,MaxIter)

%solve u_xx + u_yy + g(x,y)u = f(x,y)

% over the region D = [x0,xf,y0,yf] = {(x,y) |x0 <= x <= xf, y0 <= y <= yf}

% with the boundary Conditions:

% u(x0,y) = bx0(y), u(xf,y) = bxf(y)

% u(x,y0) = by0(x), u(x,yf) = byf(x)

% Mx = # of subintervals along x axis

% My = # of subintervals along y axis

% tol : error tolerance

% MaxIter: the maximum # of iterations

x0 = D(1); xf = D(2); y0 = D(3); yf = D(4);

dx = (xf - x0)/Mx; x = x0 + [0:Mx]*dx;

dy = (yf - y0)/My; y = y0 + [0:My]’*dy;

Mx1 = Mx + 1; My1 = My + 1;

%Boundary conditions

for m = 1:My1, u(m,[1 Mx1])=[bx0(y(m)) bxf(y(m))]; end %left/right side

for n = 1:Mx1, u([1 My1],n) = [by0(x(n)); byf(x(n))]; end %bottom/top

%initialize as the average of boundary values

sum_of_bv = sum(sum([u(2:My,[1 Mx1]) u([1 My1],2:Mx)’]));

u(2:My,2:Mx) = sum_of_bv/(2*(Mx + My - 2));

for i = 1:My

for j = 1:Mx

F(i,j) = f(x(j),y(i)); G(i,j) = g(x(j),y(i));

end

end

dx2 = dx*dx; dy2 = dy*dy; dxy2 = 2*(dx2 + dy2);

rx = dx2/dxy2; ry = dy2/dxy2; rxy = rx*dy2;

for itr = 1:MaxIter

for j = 2:Mx

for i = 2:My

u(i,j) = ry*(u(i,j + 1)+u(i,j - 1)) + rx*(u(i + 1,j)+u(i - 1,j))...

+ rxy*(G(i,j)*u(i,j)- F(i,j)); %Eq.(9.1.5a)

end

end

if itr > 1 & max(max(abs(u - u0))) < tol, break; end

u0 = u;

end

以上是possion.m文件，下面给个例子。

%solve_poisson in Example 9.1

f = inline(’0’,’x’,’y’); g = inline(’0’,’x’,’y’);

x0 = 0; xf = 4; Mx = 20; y0 = 0; yf = 4; My = 20;

bx0 = inline(’exp(y) - cos(y)’,’y’); %(E9.1.2a)

bxf = inline(’exp(y)*cos(4) - exp(4)*cos(y)’,’y’); %(E9.1.2b)

by0 = inline(’cos(x) - exp(x)’,’x’); %(E9.1.3a)

byf = inline(’exp(4)*cos(x) - exp(x)*cos(4)’,’x’); %(E9.1.3b)

D = [x0 xf y0 yf]; MaxIter = 500; tol = 1e-4;

[U,x,y] = poisson(f,g,bx0,bxf,by0,byf,D,Mx,My,tol,MaxIter);

clf, mesh(x,y,U), axis([0 4 0 4 -100 100])

更多相关内容
• 求解偏微分的代码偏微分方程求解器 此 MATLAB 代码用于可视化存在振动欧拉梁时流体域的压力和速度场。 求解器使用有限差分来求解梁的四阶微分方程。 流体是根据分析推导实现的，并与结构振动耦合。
• 不连续伽辽金法是一种求解各种偏微分方程的数值方法。 此代码的开发基于由Jan S. Hesthaven和Tim Warburton在其著作Nodal Discontinuous Galerkin Method 中开发的MATLAB代码版本。 此代码是使用 python 库开发的，...
• 描述这些过程的偏微分方程具有这样的性质；若初始时刻t=t0的已给定，则t>t0时刻的完全取决于初始条件和某些边界条件。利用MATLAB有限差分法这类问题，就是从初始值出发，通过差分格式沿时间增加的方向，逐步求...
• 有限元法的Matlab实现，以求解偏微分方程。 主程序 (FEM_Basico) 调用多个 Matlab 函数，以加载网格、执行计算和可视化结果。 该问题通过实现不同的拉格朗日元素来解决，但上传的代码特定于拉格朗日 P3。
• MATLAB求解偏微分方程（扩散方程）有限差分法，处理偏微分方程
• 是一系列已实现的神经网络，用于求解偏微分方程。 该项目属于作者在米兰理工大学航空工程中的硕士论文“求解偏微分方程的现代非确定性方法：机器学习应用于纳维-斯托克斯方程”。 为了完整起见，这里报告了硕士论文...
• 代码为2018年全国数学建模竞赛a题第一问matlab源程序，可动态生成三层隔热服距离与温度的关系图，以及三层隔热服的温度分布图。主要内容建立在一维非稳态热传导以及偏微分方程的解法程序
• 求解偏微分的代码薛定谔偏微分方程：数值和模拟。 作者：Alejandro Galván、Daniel Quiles 和 Bernat Ramis - UPC 数学和工程物理专业的学生。 这项工作的目的是研究数值方法对求解薛定谔偏微分方程的适用性。 ...
• 求解偏微分的代码数值偏微分方程 包含有限差分方法对偏微分方程的实际实现的存储库。 这个存储库包含我在约克完成的一个模块中的一些项目，以及我自己的一些工作。 在这个存储库中编写的代码应该能够在 MATLAB 和 ...
• 零基础使用 MATLAB 求解偏微分方程（建议收藏） 文章目录零基础使用 MATLAB 求解偏微分方程（建议收藏）偏微分开源工具介绍PDE 工具箱函数汇总介绍0 基础：GUI 界面操作示例问题工具箱求解导出为代码形式代码导出...

# 零基础使用 MATLAB 求解偏微分方程（建议收藏）

## 偏微分开源工具介绍

百分之九十以上的重要的工程和数学科学研究，和偏微分方程都脱不开关系。在所有的偏微分方程中，百分之九十九都是没有解析解的。没有解析解怎么办，我们只能通过有限元或者有限差分等方法，求解偏微分方程数值解。如果您有一些代码基础，建议参考我的几篇有限经典博文，简单问题可在此基础上进行修改。

对于做工程的朋友，不会偏微分方程数值解，怎么办？没关系，我推荐一些求解各类偏微分方程的容易入门的开源的软件包和工具，它们是：

• Free FEM++（足够傻瓜又不失自由度，强烈建议做工程的朋友可以学习一下）
• FEniCS(C++/Python, 开始于芝加哥大学和查尔姆斯理工大学)
• PETSC (C/Python, 美国阿贡国家实验室)
• deal.II (C++, 开始于德国海德堡大学)
• MFEM (C++, 美国劳伦斯利弗莫尔国家实验室)
• PHG (C, 张林波, 中国科学院)
• AFEPACK (C++, 李若, 北京大学）
• FEALPy（Python，魏华祎，湘潭大学）
• IFEM (MATLAB, 陈龙, UCI)
• NGSolve（C++/Python，Christoph 等）
• PHOEBESolver（ Fortran and C/C++，宁夏大学，葛永斌）
[1]: GCGE（C/MATLAB，中国科学院，谢和虎，特征值求解）
• ……

当然，商业有限元软件 ansys 等，也非常推荐学工程的朋友去学习，如果要深挖算法的，建议还是用开源的。

好，有的同学说，这些对你们来说还是太难了。没关系，我可以祭出大招：MATLAB PDE工具箱。为什么它比上面的简单呢？主要是因为，它有可视化的 GUI 工具，你实在不会写代码，你用鼠标点点，也能 “写” 出像模像样的代码。

## PDE 工具箱函数汇总介绍

PDE 工具箱包含比较多的工具，典型的几个函数如下所示。

% 偏微分方程工具箱
%
% 使用结构化的工作流定义和求解 PDEs
%   createpde               - 创建 PDE 分析模型
%   geometryFromEdges       - 从 DECSG or PDEGEOM 创建 2D 几何图形
%   importGeometry          - 从 STL 文件创建 3D 几何图形
%   geometryFromMesh        - 从三角网格创建几何图形
%   multicuboid             - 组合立方体胞单元创建 3D 几何图形
%   multicylinder           - 组合若干柱状胞单元创建 3D 几何图形
%   multisphere             - 组合若干球单元创建 3D 几何图形
%   specifyCoefficients     - 指定区域和或者子区域的 PDE 系数
%   applyBoundaryCondition  - 给几何区域施加边界条件
%   setInitialConditions    - 设定 PDE 初始条件
%   generateMesh            - 从几何生成一个网格
%   solvepde                - 求解 PDE
%   solvepdeeig             - 求解 PDE 特征值问题
%   assembleFEMatrices      - 组装中间的有限元矩阵
%   createPDEResults        - 创建一个用于后处理的结果对象
%   pdegplot                - 绘制 PDE 几何表示
%   pdemesh                 - 绘制 PDE 网格
%   pdeplot                 - 绘制二维 PDE 网格和结果
%   pdeplot3D               - 绘制 3D PDE 网格和结果
%   interpolateSolution     - 在指定的空间位置插入解
%   evaluateCGradient       - 评估 PDE 解的通量
%
% 使用热模型解决以传导为主的传热问题
%   thermalProperties       - 为热模型指定材料的热属性
%   internalHeatSource      - 指定热模型的内部热源
%   thermalBC               - 指定热模型的边界条件
%   thermalIC               - 设置热模型的初始条件或初始猜测
%   solve                   - 求解热模型中指定的传热问题
%   interpolateTemperature	- 在任意空间位置的热结果中插入温度
%   evaluateHeatFlux        - 在节点或任意空间位置评估热解的热通量
%   evaluateHeatRate        - 评估垂直于指定边界的综合热流率
%
% 使用结构模型解决静态、模态和瞬态线性弹性问题
%   structuralProperties    - 为模型分配结构材料属性
%   structuralBC            - 将边界条件应用于结构模型
%   structuralIC            - 设置初始位移和速度
%   structuralDamping       - 为结构模型指定比例阻尼参数
%   solve                   - 求解 StructuralModel 中指定的结构模型
%   evaluateStress          - 评估节点位置的应力
%   evaluateStrain          - E评估节点位置的应变
%   evaluateVonMisesStress  - 评估节点位置的 von Mises 应力
%   evaluatePrincipalStrain - 计算节点位置的主要应变
%   evaluatePrincipalStress - 计算节点位置的主应力
%   evaluateReaction        - 评估边界上的反作用力
%   interpolateDisplacement - 在指定的空间位置插入位移
%   interpolateVelocity     - 在指定的空间位置插入速度
%   interpolateAcceleration - 在指定的空间位置插入加速度
%   interpolateStress       - 在指定的空间位置插入应力
%   interpolateStrain       - 在指定的空间位置插入应变
%   interpolateVonMisesStress - 在指定空间位置内插 von Mises 应力
%
% 使用非结构化工作流程求解 PDE
%   adaptmesh   - 自适应网格生成和 PDE 解
%   assema      - 组装面积积分贡献
%   assemb      - 组装边界条件贡献
%   assempde    - 组装 PDE 问题
%   hyperbolic  - 解决双曲线问题
%   parabolic   - 解决抛物线问题
%   pdeeig      - 解决特征值 PDE 问题
%   pdenonlin   - 解决非​​线性 PDE 问题
%   poisolv     - 矩形网格上泊松方程的快速解
%
% 用户界面算法和实用程序
%   pdecirc     - 画圆
%   pdeellip    - 绘制椭圆
%   pdemdlcv    - 转换 MATLAB 4.2c 模型 MATLAB 文件以与 MATLAB 5 一起使用
%   pdepoly     - 绘制多边形
%   pderect     - 绘制矩形.
%   pdeModeler  - PDE Modeler 图形用户界面 (GUI)
%
% 几何算法
%   csgchk      - 检查几何描述矩阵的有效性
%   csgdel      - 删除最小区域之间的边界
%   decsg       - 将构造实体几何分解为最小区域
%   initmesh    - 构建初始三角形网格
%   jigglemesh  - 抖动三角形网格的内部点
%   pdearcl     - 参数化表示和弧长之间的插值
%   poimesh     - 在矩形几何体上制作规则网格
%   refinemesh  - 细化三角形网格
%   wbound      - 写入边界条件规范数据文件
%   wgeom       - W写入几何规格数据文件
%
% 绘图函数
%   pdecont     - 等高线图的速记命令
%   pdegplot    - 绘制 PDE 几何
%   pdemesh     - 绘制 PDE 三角形网格
%   pdeplot     - 通用 PDE 工具箱绘图函数
%   pdesurf     - 曲面图的速记命令
%
% 实用算法
%   dst         - 离散正弦变换
%   idst        - 逆离散正弦变换
%   pdecgrad    - 计算 PDE 解的通量
%   pdeent      - 与给定三角形集相邻的三角形的索引
%   pdegrad     - 计算 PDE 解的梯度
%   pdeintrp    - 将函数值插入到三角形中点
%   pdejmps     - 适应的误差估计
%   pdeprtni    - 将函数值内插到网格节点
%   pdesde      - 与一组子域相邻的边的索引
%   pdesdp      - 一组子域中的点索引
%   pdesdt      - 一组子域中的三角形索引
%   pdesmech    - 计算结构力学张量函数
%   pdetrg      - 三角形几何数据
%   pdetriq     - 测量网格三角形的质量
%   poiasma     - 泊松方程的边界点矩阵贡献
%   poicalc     - 矩形网格上泊松方程的快速解
%   poiindex    - 矩形网格的规范排序点的索引
%   sptarn      - 求解决广义稀疏特征值问题
%   tri2grid    - 从 PDE 三角形网格插值到矩形网格
%
% 用户定义的算法
%   pdebound    - 边界 MATLAB 文件
%   pdegeom     - 几何 MATLAB 文件
%
% 对象创建函数。这些函数不是直接调用的。
%   PDEModel          - 表示 PDE 模型的容器
%   GeometricModel    - 模型边界的几何表示
%   AnalyticGeometry  - 来自 PDEGEOM 或 DECSG 几何矩阵的 2D 几何对象
%   DiscreteGeometry  - 分面边界的几何表示
%   BoundaryCondition - 定义 PDE 的边界条件
%   CoefficientAssignmentRecords  - 方程系数的分配
%   CoefficientAssignment         - 指定区域或子域上的所有 PDE 系数
%   InitialConditionsRecords   - 记录初始条件的分配
%   GeometricInitialConditions - 区域或区域边界上的初始条件
%   NodalInitialConditions  - 在网格节点指定的初始条件
%   PDEResults           - PDE 解及其派生量
%   StationaryResults    - PDE 解及其派生量
%   TimeDependentResults - PDE 解及其派生量
%   EigenResults         - PDE 解表示
%   StructuralModel      - 表示结构分析模型的容器
%   ThermalModel         - 表示热分析模型的容器
%   ThermalMaterialAssignment - 指定区域或子区域的材料属性
%   HeatSourceAssignment - 指定域或子域上的热源
%   ThermalBC            - 定义热模型的边界条件 (BC)
%   GeometricThermalICs  - 区域或区域边界上的初始温度
%   NodalThermalICs      - 在网格节点指定的初始温度
%   ThermalResults       - 热解及其派生量
%   TransientThermalResults   - 瞬态热模型解及其派生量
%   StructuralMaterialAssignment - 区域或子域上的结构材料属性分配
%   StructuralBC         - 定义结构模型的边界载荷或边界条件 (BC)
%   StructuralResults    - 结构解及其派生量
%   StaticStructuralResults - 静态结构模型解及其派生量
%   StructuralDampingAssignment - 结构分析模型的阻尼分配
%   GeometricStructuralICs - 区域上的初始位移和速度
%   NodalStructuralICs - 在网格节点指定的初始位移和速度
%   ModalStructuralResults - 结构模态分析结果
%   TransientStructuralResults - 瞬态结构模型解及其派生量
%

% 未记录的类和函数
%   pdeCalcFullU
%   pdeODEInfo
%   pdeParabolicInfo
%   pdeHyperbolicInfo



## 0 基础：GUI 界面操作

#### 示例问题

没有什么编程基础，但是又想快速写出有限元程序的同学，建议使用图形界面进行编程，然后导出代码。做个简单的示例操作。比如要求解：
− Δ u = λ u u ∣ ∂ Ω = 0 Ω 是 一 个 L 型 区 域 ， 如 下 图 所 示 -\Delta u = \lambda u\\ u|_{\partial \Omega}=0\\ \Omega 是一个L 型区域，如下图所示

#### 工具箱求解

• 打开 MATLAB
• 命令行窗口口输入 pdetool 回车
• 依次点击菜单栏如下按钮，其中点击 PDE 的时候，改成特征值模式

基础通过 GUI 界面生成的代码此由 pdetool 编写和读取，不应编辑。 有两个推荐的替代方案：

• 从 pdetool 导出所需的变量并创建一个 MATLAB 脚本，对这些变量执行操作。
• 使用 MATLAB 脚本完全定义问题。

#### 导出为代码形式

得到求解结果后，保存为 main.m 文件，并打开。

function pdemodel
[pde_fig,ax]=pdeinit;
pdetool('appl_cb',1);
set(ax,'DataAspectRatio',[1.5 1 1]);
set(ax,'PlotBoxAspectRatio',[1 0.74375917767988253 0.74375917767988253]);
set(ax,'XLim',[-1.5 1.5]);
set(ax,'YLim',[-1 1]);
set(ax,'XTickMode','auto');
set(ax,'YTickMode','auto');

% Geometry description:
pdepoly([ -1,...
1,...
1,...
0,...
0,...
-1,...
],...
[ -1,...
-1,...
1,...
1,...
0,...
0,...
],...
'P1');
set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','P1')

% Boundary conditions:
pdetool('changemode',0)
pdesetbd(6,...
'dir',...
1,...
'1',...
'0')
pdesetbd(5,...
'dir',...
1,...
'1',...
'0')
pdesetbd(4,...
'dir',...
1,...
'1',...
'0')
pdesetbd(3,...
'dir',...
1,...
'1',...
'0')
pdesetbd(2,...
'dir',...
1,...
'1',...
'0')
pdesetbd(1,...
'dir',...
1,...
'1',...
'0')

% Mesh generation:
setappdata(pde_fig,'refinemethod','regular');
setappdata(pde_fig,'jiggle',char('on','mean',''));
setappdata(pde_fig,'MesherVersion','preR2013a');
pdetool('initmesh')
pdetool('refine')

% PDE coefficients:
pdeseteq(4,...
'1.0',...
'0.0',...
'10.0',...
'1.0',...
'0:10',...
'0.0',...
'0.0',...
'[0 100]')
setappdata(pde_fig,'currparam',...
['1.0 ';...
'0.0 ';...
'10.0';...
'1.0 '])

% Solve parameters:
setappdata(pde_fig,'solveparam',...
'0.5','longest','0','1E-4','','fixed','Inf'))

% Plotflags and user data strings:
setappdata(pde_fig,'plotflags',[1 1 1 1 1 1 1 1 0 0 0 1 1 0 0 0 0 1]);
setappdata(pde_fig,'colstring','');
setappdata(pde_fig,'arrowstring','');
setappdata(pde_fig,'deformstring','');
setappdata(pde_fig,'heightstring','');

% Solve PDE:
pdetool('solve')


#### 代码导出相关数据

当前目录下，保存如下代码为 matqueque。

function y = matqueue(p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15,p16,p17,p18)
%MATQUEUE Creates and manipulates a figure-based matrix queue.
%   FIG = MATQUEUE('create');
%   Create a queue figure and return its number.
%
%   FIG = MATQUEUE('find');
%   Searches the root window's children to find the queue
%   figure.  Returns 0 if no queue exists.
%
%   MATQUEUE('put', X1, X2, ..., X18);
%   Insert up to 18 matrices into the queue.  Create the
%   queue if none exists.
%
%   X = MATQUEUE('get');
%   Get a matrix out the queue.  Return [] if the queue is
%   empty.
%
%   NUM_ITEMS = MATQUEUE('length');
%   Return the number of matrices in the queue.  Return -1 if
%   no buffer exists.
%
%   MATQUEUE('clear')
%   Empty the queue.
%
%   MATQUEUE('close')
%   Close the queue figure.

buffer_name = 'FIFO Buffer';

if (nargin < 1)
action = 'create';
else
action = lower(p0);
end

if (strcmp(action, 'create'))

%==================================================================
% Create a new queue.
%
% matqueue('create');
%==================================================================

narginchk(1,1);
nargoutchk(0,1);

oldFig = findobj(allchild(0), 'flat', 'Visible', 'on');

buffer_fig = matqueue('find');
if (buffer_fig ~= 0)
% Buffer already exists; do nothing
return;
end

buffer_fig = figure('Name', buffer_name, ...
'Visible', 'off',...
'HandleVisibility', 'callback', ...
'IntegerHandle', 'off', ...
'NumberTitle', 'off', ...
'Tag', buffer_name);
if (~isempty(oldFig))
figure(oldFig(1));
end

queue_holder = uicontrol(buffer_fig, 'Style', 'text', ...
'Visible', 'off', 'Tag', 'QueueHolder');

y = buffer_fig;

return;

elseif (strcmp(action, 'find'))

%==================================================================
% Find the queue figure.  If no queue figure exists, return 0.
%
% matqueue('find');
%==================================================================

narginchk(1,1);
nargoutchk(0,1);

% Search the root's children for a figure with the right name
buffer_number = findobj(allchild(0), 'flat', 'Tag', buffer_name);
if (isempty(buffer_number))
y = 0;
else
y = buffer_number(1);
end

return;

elseif (strcmp(action, 'put'))

%==================================================================
% Put matrices into the queue.  Queue figure is created if none
% exists.
%
% matqueue('put', X1, X2, ..., X18);
%==================================================================

narginchk(2,19);
nargoutchk(0,0);

buffer_fig = matqueue('find');
if (buffer_fig == 0)
buffer_fig = matqueue('create');
end

queue_holder = findobj(get(buffer_fig, 'Children'), 'flat', 'Tag', 'QueueHolder');
if (isempty(queue_holder))
error(message('pde:matqueue:corruptMatrixQueue'));
end

handles = get(queue_holder, 'UserData');

num_inputs = nargin-1;
new_handles = zeros(1, num_inputs);
for i = 1:num_inputs
arg_name = ['p', num2str(i)];
try_string = ['new_handles(num_inputs+1-i)=uicontrol(buffer_fig,', ...
' ''Style'',''text'',''Visible'','...
' ''off'',''UserData'',', arg_name, ');'];
eval(try_string);
end

set(queue_holder, 'UserData', [new_handles handles]);

return;

elseif (strcmp(action, 'get'))

%==================================================================
% Return earliest matrix item in the queue.  Errors out if there's
% no queue.  Returns empty matrix if queue is empty.
%
% X = matqueue('get');
%==================================================================

narginchk(1,1);
nargoutchk(0,1);

buffer_fig = matqueue('find');
if (buffer_fig == 0)
% No buffer; return empty matrix.
y = [];
return;
end

queue_holder = findobj(get(buffer_fig, 'Children'), 'flat', 'Tag', 'QueueHolder');
if (isempty(queue_holder))
error(message('pde:matqueue:corruptMatrixQueue'));
end

handles = get(queue_holder, 'UserData');
N = length(handles);
if (N > 0)
y = get(handles(N), 'UserData');
delete(handles(N));
handles(N) = [];
set(queue_holder, 'UserData', handles);
else
% Nothing in the buffer; return empty matrix
y = [];
end

return;

elseif (strcmp(action, 'length'))

%==================================================================
% Returns the length of the queue.  Returns -1 if no queue
% figure exists.
%
% num_items = matqueue('length');
%==================================================================

narginchk(1,1);
nargoutchk(0,1);

buffer_fig = matqueue('find');
if (buffer_fig == 0)
% No buffer!  Return 0.
y = 0;
return;
end

queue_holder = findobj(get(buffer_fig, 'Children'), 'flat', 'Tag', 'QueueHolder');
if (isempty(queue_holder))
error(message('pde:matqueue:corruptMatrixQueue'));
end

y = length(get(queue_holder, 'UserData'));

return;

elseif (strcmp(action, 'clear'))

%==================================================================
% Clear the queue.
%
% matqueue('clear');
%==================================================================

narginchk(1,1);
nargoutchk(0,1);

buffer_fig = matqueue('find');
if (buffer_fig == 0)
% No buffer; nothing to do.
return;
end

queue_holder = findobj(get(buffer_fig, 'Children'), 'flat', 'Tag', 'QueueHolder');
if (~isempty(queue_holder))
delete(queue_holder);
end

queue_holder = uicontrol(buffer_fig, 'Style', 'text', ...
'Visible', 'off', 'Tag', 'QueueHolder');

return;

elseif (strcmp(action, 'close'))

%==================================================================
% Close the queue figure.
%
% matqueue('close');
%==================================================================

narginchk(1,1);
nargoutchk(0,1);

buffer_fig = matqueue('find');
if (buffer_fig ~= 0)
close(buffer_fig);
end

return;

else

error(message('pde:matqueue:unknownAction'));

end



同样地，将如下代码存为 matqdlg。

function y = matqdlg(P0,P1,V1,P2,V2,P3,V3,P4,V4,P5,V5,P6,V6,P7,V7,P8,V8,P9,V9)
%MATQDLG Workspace transfer dialog box.
%   MATQDLG('ws2buffer', {prop/value pairs})
%   Put up a dialog box that invites the user to enter a comma-separated
%   list of expressions.  When user clicks OK,  eval the expressions
%   one at a time, putting the results into the buffer.  Allowable
%   properties include 'PromptString', which may be a string matrix,
%   'OKCallback', which will be eval'ed when the user finishes
%   with the dialog box by typing <Return> in the entry field or
%   clicking OK, 'CancelCallback', which will be eval'ed when the user
%   clicks Cancel, and 'EntryString', the default user entry.
%   Figure properties are also allowed in this list, such as 'Name'.
%
%   MATQDLG('buffer2ws', {prop/value pairs})
%   Put up a dialog box that invites the user to enter N comma-separated
%   variable names, where N is the number of items in the buffer.  Get
%   items out of the buffer one at a time, storing the results in
%   indicated workspace variables.  Allowable properties include
%   'PromptString', 'OKCallback', 'CancelCallback', and 'EntryString'.
%   These work the same way as in the 'ws2buffer' action.
%
%   Y = MATQDLG('get_entry');
%   Return the user-entered string in the entry field of the workspace
%   transfer dialog box.
%
%   H = MATQDLG('find')
%   Return the handle of the dialog box figure.
%
%   H = MATQDLG('create')
%   Create the dialog box figure, returning its handle.

% MATLAB-files required:  matqparse.m, ws2matq.m, matq2ws.m.

buffer_tag = 'Workspace Transfer';
buffer_name = '';

if (nargin < 1)
action = 'create';
else
action = lower(P0);
end

if (strcmp(action, 'create'))

%==================================================================
% Create a new queue.
%
% matqdlg('create');
%==================================================================

narginchk(1,1);
nargoutchk(0,1);

buffer_fig = matqdlg('find');
if (buffer_fig ~= 0)
% Queue already exists; do nothing
return;
end

screenSize = get(0,'ScreenSize');
width = 415;
height = 136;
left = (screenSize(3) - width) / 2;
bottom = (screenSize(4) - height) / 2;

buffer_fig = figure('Name', buffer_name, 'Visible', 'off',...
'HandleVisibility', 'callback', ...
'IntegerHandle', 'off', ...
'Units', 'pixels', ...
'Position', [left bottom width height], ...
'Colormap', [], ...
'Color', get(0, 'DefaultUIControlBackgroundColor'), ...
'DefaultUIControlInterruptible','on', ...
'Tag', buffer_tag, ...
'NumberTitle', 'off');

axes('Visible', 'off', 'Parent', buffer_fig);

ok = uicontrol(buffer_fig, ...
'Style', 'push', ...
'Units', 'normalized', ...
'Position', [.63 .06 .15 .24], ...
'Tag', 'OK', ...
'String', 'OK');

cancel = uicontrol(buffer_fig, ...
'Style', 'push', ...
'Units', 'normalized', ...
'Position', [.80 .06 .15 .24], ...
'Tag', 'Cancel', ...
'String', 'Cancel');

prompt = uicontrol(buffer_fig, ...
'Style', 'text', ...
'Units', 'normalized', ...
'Position', [.05 .76 .90 .15], ...
'String', '', ...
'Min', 1, ...
'Max', 3, ...
'Tag', 'Prompt', ...
'Horizontal', 'left');

entry = uicontrol(buffer_fig, ...
'Style', 'edit', ...
'Units', 'normalized', ...
'Position', [.05 .46 .90 .31], ...
'BackgroundColor', 'w', ...
'ForegroundColor', 'k', ...
'Tag', 'Entry', ...
'Horizontal', 'left');

y = buffer_fig;

return;

elseif (strcmp(action, 'find'))

%==================================================================
% Find the queue figure.  If no queue figure exists, return 0.
%
% matqdlg('find');
%==================================================================

narginchk(1,1);
nargoutchk(0,1);

% Search the root's children for a figure with the right tag
buffer_number = findobj(allchild(0), 'flat', 'Type', 'figure', ...
'Tag', buffer_tag);
if (isempty(buffer_number))
y = 0;
else
y = buffer_number(1);
end

return;

elseif (strcmp(action, 'ws2buffer'))

%==================================================================
% Invoke the dialog box in workspace-to-buffer mode.
%
% matqdlg('ws2buffer')
%==================================================================

narginchk(1,19);
if (rem(nargin,2) ~= 1)
error(message('pde:matqdl:invalidNumberOfArgs'));
end

% Remember the current visible figure.
figHandles = findobj(allchild(0), 'flat', 'Visible', 'on');

% Set up default properties.
ok_string = '';
cancel_string = '';
entry_string = '';
prompt_string = 'Enter workspace variable names or expressions:';

buffer_fig = matqdlg('find');
if (buffer_fig == 0)
buffer_fig = matqdlg('create');
end
if (~isempty(figHandles))
set(buffer_fig, 'UserData', figHandles(1));
end

% Process param/value pairs.
num_properties = (nargin - 1)/2;
for k = 1:num_properties
prop_arg_name = ['P' num2str(k)];
val_arg_name = ['V' num2str(k)];
prop_arg = lower(eval(prop_arg_name));
val_arg = eval(val_arg_name);
if (strcmp(prop_arg, 'promptstring'))
prompt_string = val_arg;
elseif (strcmp(prop_arg, 'okcallback'))
ok_string = val_arg;
elseif (strcmp(prop_arg, 'cancelcallback'))
cancel_string = val_arg;
elseif (strcmp(prop_arg, 'entrystring'))
entry_string = val_arg;
else
set(buffer_fig, prop_arg, val_arg);
end
end

promptButton = findobj(buffer_fig, 'Tag', 'Prompt');
entry = findobj(buffer_fig, 'Tag', 'Entry');
okButton = findobj(buffer_fig, 'Tag', 'OK');
cancelButton = findobj(buffer_fig, 'Tag', 'Cancel');
set(okButton, 'UserData', ok_string, 'Callback', {@okButtonCallback,'ws'});
set(cancelButton, 'UserData', cancel_string, 'Callback', @cancelButtonCallback);
set(promptButton, 'String', prompt_string);
set(entry, 'String', entry_string);

% Adjust height of figure and prompt box.
prompt_lines = size(prompt_string,1);
if (prompt_lines > 1)
set([promptButton entry okButton cancelButton], 'Units', 'pixels');
prompt_position = get(promptButton, 'Position');
prompt_height = prompt_position(4);
increase = (prompt_lines-1) * prompt_height;
fig_position = get(buffer_fig, 'Position');
set(promptButton, 'Position', prompt_position + [0 0 0 increase]);
set(buffer_fig, 'Position', fig_position + [0 0 0 increase]);
set([promptButton entry okButton cancelButton], 'Units', 'normalized');
end

drawnow;
figure(buffer_fig);

return

elseif (strcmp(action, 'buffer2ws'))

%==================================================================
% Invoke the dialog box in workspace-to-buffer mode.
%
% matqdlg('ws2buffer')
%==================================================================

narginchk(1,19);
if (rem(nargin,2) ~= 1)
error(message('pde:matqdl:invalidNumberOfArgs'));
end

% Remember the current visible figure.
figHandles = findobj(allchild(0), 'flat', 'Visible', 'on');

num_items = matqueue('length');
if (num_items == 0)
% If the queue is empty, there's nothing to do!
error(message('pde:matqdl:emptyQueue'));
end

buffer_fig = matqdlg('find');
if (buffer_fig == 0)
buffer_fig = matqdlg('create');
end
if (~isempty(figHandles))
set(buffer_fig, 'UserData', figHandles(1));
end

% Set up default properties.
ok_string = '';
entry_string = '';
cancel_string = '';
if (num_items == 1)
prompt_string = 'Enter a variable name:';
else
prompt_string = sprintf('Enter %d variable names:', num_items);
end

% Process param/value pairs.
num_properties = (nargin - 1)/2;
for k = 1:num_properties
prop_arg_name = ['P' num2str(k)];
val_arg_name = ['V' num2str(k)];
prop_arg = lower(eval(prop_arg_name));
val_arg = eval(val_arg_name);
if (strcmp(prop_arg, 'promptstring'))
prompt_string = val_arg;
elseif (strcmp(prop_arg, 'okcallback'))
ok_string = val_arg;
elseif (strcmp(prop_arg, 'cancelcallback'))
cancel_string = val_arg;
elseif (strcmp(prop_arg, 'entrystring'))
entry_string = val_arg;
else
set(buffer_fig, prop_arg, val_arg);
end
end

promptButton = findobj(buffer_fig, 'Tag', 'Prompt');
set(findobj(buffer_fig, 'Tag', 'OK'), 'UserData', ok_string, ...
'Callback', {@okButtonCallback,'mat'});
set(findobj(buffer_fig, 'Tag', 'Cancel'), 'UserData', cancel_string, ...
'Callback', @cancelButtonCallback);
set(promptButton, 'String', prompt_string);
set(findobj(buffer_fig, 'Tag', 'Entry'), 'String', entry_string);

% Adjust height of figure and prompt box.
prompt_lines = size(prompt_string,1);
if (prompt_lines > 1)
set([promptButton entry okButton cancelButton], 'Units', 'pixels');
prompt_position = get(promptButton, 'Position');
prompt_height = prompt_position(4);
increase = (prompt_lines-1) * prompt_height;
fig_position = get(buffer_fig, 'Position');
set(promptButton, 'Position', prompt_position + [0 0 0 increase]);
set(buffer_fig, 'Position', fig_position + [0 0 0 increase]);
set([promptButton entry okButton cancelButton], 'Units', 'normalized');
end

drawnow;
figure(buffer_fig);

return;

elseif (strcmp(action, 'get_entry'))

%==================================================================
% Get the user-entered string in the entry field.
%
% string = matqdlg('get_entry');
%==================================================================

buffer_fig = matqdlg('find');
if (buffer_fig < 1)
y = [];
return;
end

y = get(findobj(buffer_fig, 'Tag', 'Entry'), 'String');

return;

else

error(message('pde:matqdl:invalidAction'));

end

%--------------------------------------------
function   okButtonCallback(obj,evd,wsormat)

switch(wsormat)
case 'ws'
ws2matq

case 'mat'
matq2ws

otherwise
error(message('pde:matqdl:okButtonCallback:unknownOption'));

end
%-------------------------------------------

function cancelButtonCallback(obj,evd)

matqueue('clear');
eval(get(findobj(gcbf,'Tag','Cancel'),'UserData'));
close(matqdlg('find'));



将如下代码存为文件 matq2ws 。

%MATQ2WS Helper script for matqdlg.
%  MATQ2WS gets the user-entered comma-delimited string,
%  parses it, and then tries to put the queue contents one
%  at a time into the resulting variable names.  For
%  recoverable errors, the prompt is reset and the user can
%  try again.  Recoverable errors include empty input
%  string, string containing "#", too few variable names,
%  and too many variable names.  If the user types something
%  that cannot be a workspace variable name, that's a
%  nonrecoverable error.  The queue is cleared and made
%  invisible, and an error message is printed to the command
%  window.

% Variable names (these need to be cleared before returning):
% var_string_ err_string_ new_prompt_ pound_ N_
% fatal_error_flag_ i_ expr_ try_string_ catch_string_ error_message_

var_string_ = get(findobj(gcbf,'Tag','Entry'), 'String');
[var_string_, err_string_] = matqparse(var_string_);
if (~isempty(err_string_))
errordlg(char('Could not parse your expression.', err_string_), ...
'Workspace Transfer Error', 'on');
clear var_string_ err_string_ new_prompt_ ...
pound_ N_ fatal_error_flag_ i_ expr_ try_string_ catch_string_ ...
error_message_
return;
end
N_ = size(var_string_, 1);
if (N_ < matqueue('length'))
errordlg(char('You did not enter enough variable names.', ...
'Please try again.'), 'Workspace Transfer Error', 'on');
clear var_string_ err_string_ new_prompt_ ...
pound_ N_ fatal_error_flag_ i_ expr_ try_string_ catch_string_ ...
error_message_
return;
elseif (N_ > matqueue('length'))
errordlg(char('You entered too many variable names.', ...
'Please try again.'), 'Workspace Transfer Error', 'on');
clear var_string_ err_string_ new_prompt_ ...
pound_ N_ fatal_error_flag_ i_ expr_ try_string_ catch_string_ ...
error_message_
return;
end
fatal_error_flag_ = 0;
for i_ = 1:N_
expr_ = deblank(var_string_(i_, :));
try
assignin('base', expr_, matqueue('get'));
catch
fatal_error_flag_ = 1;
end

if (fatal_error_flag_)
errordlg(char(sprintf('Error using "%s" as a workspace variable.', ...
expr_), 'You will need to start over.'), ...
'Workspace Transfer Error', 'on');
set(matqdlg('find'), 'Visible', 'off');
if (~isempty(get(matqdlg('find'),'UserData')))
if (any(get(0,'Children') == get(matqdlg('find'),'UserData')))
figure(get(matqdlg('find'),'UserData'));
end
end
matqueue('clear');
clear var_string_ err_string_ new_prompt_ ...
pound_ N_ fatal_error_flag_ i_ expr_ try_string_ catch_string_ ...
error_message_
return;
end
end
set(matqdlg('find'), 'Visible', 'off');
if (~isempty(get(matqdlg('find'),'UserData')))
if (any(get(0,'Children') == get(matqdlg('find'),'UserData')))
figure(get(matqdlg('find'),'UserData'));
end
end

try
char(get(get(matqdlg('find'),'CurrentObject'),'UserData'));
catch
disp('Error evaluating button callback.')
end
close(matqdlg('find'));
clear var_string_ err_string_ new_prompt_ ...
pound_ N_ fatal_error_flag_ i_ expr_ try_string_ catch_string_ error_message_



如下代码为 matqparse 。

function [m,error_str] = matqparse(str,flag)
%MATQPARSE Dialog entry parser for MATQDLG.
%   [M,ERROR_STR] = MATQPARSE(STR,FLAG) is a miniparser
%   for MATQDLG.
%   eg: 'abc de  f ghij' becomes [abc ]
%                                [de  ]
%                                [f   ]
%                                [ghij]
%   Uses either spaces, commas, semi-colons, or brackets
%   as separators.  Thus 'a 10*[b c] d' will crash. User
%   must instead say 'a [10*[b c]] d'.
%

% Error checks
error_str = '';
if nargin==0
error_str = getString(message('pde:matqparse:StringReqd'));
return
elseif size(str,1)>1 | ~ischar(str)
error_str = getString(message('pde:matqparse:SingleRowStringReqd'));
return
end

if nargin<2
flag = 1;
end

l = length(str);
m = '';
i = 1;
j = 1;
k = 1;
while k<=l
% Check for missing [
if str(k)==']'
error_str = getString(message('pde:matqparse:UnmatchedRightBracket'));
return
elseif str(k)=='['
% Check for missing ]
index = find(str(k+1:l)==']');
if isempty(index)
error_str = getString(message('pde:matqparse:UnmatchedLeftBracket'));
return
else
% Check for mismatched brackets between k+1 and last element
index1 = find(str(k+1:l)=='[');
l_index = length(index);
l_index1 = length(index1);
if l_index~=l_index1+1
error_str = getString(message('pde:matqparse:BracketMismatch'));
return
else
% Everything OK so far
di = find([index1 index(l_index)+1]>index);
end_ind = index(di(1));
m_middle = ['[' matqparse(str(k+1:k+end_ind-1),2) ']'];
if flag==1
% m and m_end may be multiline matrices
m_end = matqparse(str(k+end_ind+1:l),1);
m = char(m,m_middle,m_end);
else
% m and m_end will be single line
m_end = matqparse(str(k+end_ind+1:l),2);
m = [m m_middle m_end];
end
k = l+1;
end
end
elseif any(str(k)==' ;,') & (flag==1)
if j>1
% Only reset to beginning of next row if
% NOT already at beginning of a row
j=1;
i = i+1;
end
k = k+1; % Increment index into str
else
m(i,j) = str(k);
j = j+1; % Increment column of resultant matrix, m
k = k+1; % Increment index into str
end
end

% Since char of zero is end-of-string flag, change to blanks
if ~isempty(m)
EndOfString = find(abs(m)==0);
m(EndOfString) = char(' '*ones(size(EndOfString)));

% Eliminate any empty rows
if size(m,2)>1
m  = m(find(any(m'~=' ')),:);
end

end

% end matqparse



在生成的主文件求解程序末尾添加和修改为如下代码，即可导出数据。

clc
clear
close all
[pde_fig,ax]=pdeinit;
pdetool('appl_cb',1);
set(ax,'DataAspectRatio',[1.5 1 1]);
set(ax,'PlotBoxAspectRatio',[1 0.74375917767988253 0.74375917767988253]);
set(ax,'XLim',[-1.5 1.5]);
set(ax,'YLim',[-1 1]);
set(ax,'XTickMode','auto');
set(ax,'YTickMode','auto');

% Geometry description:
pdepoly([ -1,...
1,...
1,...
0,...
0,...
-1,...
],...
[ -1,...
-1,...
1,...
1,...
0,...
0,...
],...
'P1');
set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','P1')

% Boundary conditions:
pdetool('changemode',0)
pdesetbd(6,...
'dir',...
1,...
'1',...
'0')
pdesetbd(5,...
'dir',...
1,...
'1',...
'0')
pdesetbd(4,...
'dir',...
1,...
'1',...
'0')
pdesetbd(3,...
'dir',...
1,...
'1',...
'0')
pdesetbd(2,...
'dir',...
1,...
'1',...
'0')
pdesetbd(1,...
'dir',...
1,...
'1',...
'0')

% Mesh generation:
setappdata(pde_fig,'refinemethod','regular');
setappdata(pde_fig,'jiggle',char('on','mean',''));
setappdata(pde_fig,'MesherVersion','preR2013a');
pdetool('initmesh')
pdetool('refine')

% PDE coefficients:
pdeseteq(4,...
'1.0',...
'0.0',...
'10.0',...
'1.0',...
'0:10',...
'0.0',...
'0.0',...
'[0 100]')
setappdata(pde_fig,'currparam',...
['1.0 ';...
'0.0 ';...
'10.0';...
'1.0 '])

% Solve parameters:
setappdata(pde_fig,'solveparam',...
'0.5','longest','0','1E-4','','fixed','Inf'))

% Plotflags and user data strings:
setappdata(pde_fig,'plotflags',[1 1 1 1 1 1 1 1 0 0 0 1 1 0 0 0 0 1]);
setappdata(pde_fig,'colstring','');
setappdata(pde_fig,'arrowstring','');
setappdata(pde_fig,'deformstring','');
setappdata(pde_fig,'heightstring','');

% Solve PDE:
pdetool('solve')
for flag=1:6
% case: export variables to workspace
pde_fig=findobj(allchild(0),'flat','Tag','PDETool');

if flag==1
% export geometry data:
gd=get(findobj(get(pde_fig,'Children'),'flat',...
ns=getappdata(pde_fig,'objnames');
evalhndl=findobj(get(pde_fig,'Children'),'flat','Tag','PDEEval');
sf=get(evalhndl,'String');
matqueue('put',gd,sf,ns)
pstr='Variable names for geometry data, set formula, labels:';
estr='gd sf ns';
elseif flag==2
% export decomposed list, boundary conditions:
dl1=getappdata(pde_fig,'dl1');
bl=get(findobj(get(h,'Children'),'flat',...
'Tag','PDEBoundMode'),'UserData');
matqueue('put',dl1,bl)
pstr='Variable names for decomposed geometry, boundary cond''s:';
estr='g b';
elseif flag==3
% export mesh:
p=get(findobj(get(h,'Children'),'flat','Tag','PDEInitMesh'),...
'UserData');
e=get(findobj(get(h,'Children'),'flat','Tag','PDERefine'),...
'UserData');
t=get(findobj(get(h,'Children'),'flat','Tag','PDEMeshParam'),...
'UserData');
matqueue('put',p,e,t)
pstr='Variable names for mesh data (points, edges, triangles):';
estr='p e t';
elseif flag==4
% export PDE coefficients:
'UserData');
ns=getappdata(pde_fig,'ncafd');
nc=ns(1); na=ns(2); nf=ns(3); nd=ns(4);
c=params(1:nc,:);
a=params(nc+1:nc+na,:);
f=params(nc+na+1:nc+na+nf,:);
d=params(nc+na+nf+1:nc+na+nf+nd,:);
matqueue('put',c,a,f,d)
pstr='Variable names for PDE coefficients:';
estr='c a f d';
elseif flag==5
% export solution:
'UserData');
'UserData');
if isempty(l)
pstr='Variable name for solution:';
estr='u';
matqueue('put',u)
else
pstr='Variable names for solution and eigenvalues:';
estr='u l';
matqueue('put',u,l)
end
elseif flag==6
% export movie:
M=getappdata(pde_fig,'movie');
matqueue('put',M)
pstr='Variable name for PDE solution movie:';
estr='M';
end
pdeinfo('Change the variable name(s) if desired. OK when done.',0);
%matqdlg('buffer2ws','Name','Export','PromptString',pstr,...
%    'OKCallback','pdeinfo;','CancelCallback','pdeinfo;','EntryString',estr);

end



出现 You did not enter enough variable names. Please try again. 如何解决？

不要调用matqdlg('buffer2ws','Name','Export','PromptString',pstr,'OKCallback','pdeinfo;','CancelCallback','pdeinfo;','EntryString',estr);即可。

## 0.1 基础：编程调用 PDE 工具箱

有的朋友说，PDE 工具箱求解 PDE 生成的代码，运行的时候会跳出界面，一看就显得很没有 B 格，有没有办法纯代码操作呢？当然有，界面也只不过是一些代码的集合而已。不想要 PDE 工具箱的界面，又想快速地写出有限元代码，需要一点点代码和有限元基础。

由工具箱界面生成的代码，一定是和图形界面高度耦合的，我们想把图形界面去掉不显示，并不容易。所以我们考虑使用 MATLAB 脚本完全定义问题。举一个简单的例子来说明它的操作。

我们还是以拉普拉斯特征值为例：
− Δ u = λ u -\Delta u=\lambda u

代码如下：

clc
clear
close all
model = createpde();%创建PDE模型
geometryFromEdges(model,@squareg);%从边界生成几何
pdegplot(model,'EdgeLabels','on')%可视化
ylim([-1.5,1.5])
axis equal
applyBoundaryCondition(model,'dirichlet','Edge',4,'u',0);%左边界 0 狄利克雷边界条件
applyBoundaryCondition(model,'neumann','Edge',[1,3],'g',0,'q',0);%上下边界 0 纽曼边界条件
applyBoundaryCondition(model,'neumann','Edge',2,'g',0,'q',-3/4);%由边界混合边界条件，\frac{\partial u}{\partial n}-\frac{3}{4} u=0
specifyCoefficients(model,'m',0,'d',1,'c',1,'a',0,'f',0);%指定系数，表示特征值问题
r = [-Inf,10];%找小于10的特征值和特征向量
generateMesh(model,'Hmax',0.05);%生成网格
results = solvepdeeig(model,r);%在指定范围求解特征值问题，找六个特征值
l = results.Eigenvalues;%获得特征值
u = results.Eigenvectors;%获得特征值对应的特征向量
pdeplot(model,'XYData',u(:,1));%画第一个特征函数
pdeplot(model,'XYData',u(:,length(l)));%画最后一个特征函数
%l(2) - l(1) - pi^2/4
%l(5) - l(1) - pi^2


上面用的是矩形区域。当然，更复杂的区域我们可以使用 decsg。区域的矩阵表示可以参考这个链接

decsg 使用的简单示例如下：

clc
clear
close all
rect1 = [3
4
-1
1
1
-1
0
0
-0.5
-0.5];
C1 = [1
1
-0.25
0.25];
C2 = [1
-1
-0.25
0.25];
C1 = [C1;zeros(length(rect1) - length(C1),1)];
C2 = [C2;zeros(length(rect1) - length(C2),1)];
gd = [rect1,C1,C2];
ns = char('rect1','C1','C2');
ns = ns';
sf = '(rect1+C1)-C2';
[dl,bt] = decsg(gd,sf,ns);


期间的刚度矩阵和质量矩阵也可以通过 assembleFEMatrices 获得。最后给一个区域稍微复杂一点的，通俗易懂的例子吧。顺便和直接把刚度矩阵和质量矩阵导出来调用 eigs 求解做个比较。

clc
clear
close all
%% 创建模型
model = createpde;
g = decsg([1 0 0 radius]','C1',('C1')');%通过简单集合图形生成区域
geometryFromEdges(model,g);
pdegplot(model,'EdgeLabels','on')%可视化
axis equal
title 'Geometry with Edge Labels'
c = 1;a = 0;f = 0;d = 1;
specifyCoefficients(model,'m',0,'d',d,'c',c,'a',a,'f',f);
applyBoundaryCondition(model,'dirichlet','Edge',(1:4),'u',0);
generateMesh(model,'Hmax',0.2);
%% 通过导出的刚度矩阵求特征值
FEMatrices = assembleFEMatrices(model,'nullspace');
K = FEMatrices.Kc;
B = FEMatrices.B;
M = FEMatrices.M;
sigma = 10;
numberEigenvalues = 5;
[eigenvectorsEigs,eigenvaluesEigs] = eigs(K,M,numberEigenvalues,sigma);
eigenvaluesEigs = diag(eigenvaluesEigs);
[maxEigenvaluesEigs,maxIndex] = max(eigenvaluesEigs);
eigenvectorsEigs = B*eigenvectorsEigs;
%% 通过工具箱直接求特征值
r = [min(eigenvaluesEigs)*0.99 max(eigenvaluesEigs)*1.01];
result = solvepdeeig(model,r);
eigenvectorsPde = result.Eigenvectors;
eigenvaluesPde = result.Eigenvalues;
%% 对比两种方法求出的特征值和特征向量的差距
eigenValueDiff = sort(eigenvaluesPde) - sort(eigenvaluesEigs);
fprintf('Maximum difference in eigenvalues from solvepdeeig and eigs: %e\n', ...
norm(eigenValueDiff,inf));
%% 画所选范围最大特征值对应的特征函数
h = figure;
h.Position = [1 1 2 1].*h.Position;
subplot(1,2,1)
axis equal
pdeplot(model,'XYData', eigenvectorsEigs(:,maxIndex),'Contour','on')
title(sprintf('eigs eigenvector, eigenvalue: %12.4e', eigenvaluesEigs(maxIndex)))
xlabel('x')
ylabel('y')
subplot(1,2,2)
axis equal
pdeplot(model,'XYData',eigenvectorsPde(:,end),'Contour','on')
title(sprintf('solvepdeeig eigenvector, eigenvalue: %12.4e',eigenvaluesPde(end)))
xlabel('x')
ylabel('y')


展开全文
• MATLAB解偏微分方程，附代码以及原理图。含拉普拉斯方程。绝热细杆求解方程等
• Matlab 求解偏微分的代码 :rocket: 建模技术-2 ...偏微分方程求解器。 它还包括用于 2D 和 3D MATLAB 绘图的通用动画图形回放 UI，以及自动从折线图图像中读取数据的函数。 使用的MATLAB版本是R2016b
• 脚本和函数文件，用于求解单位平方上的线性二维偏微分方程，使用具有双线性元素的标准 Galerkin FEM 和两种稀疏网格方法。 此代码用于在 Russell, S. 和 Madden, N. 中生成结果。稀疏网格有限元方法的分析和实现。 ...
• 非稳态的偏微分方程组是一个比较难解决的问题，也是在热质交换等方面的常常遇到的问题，因此需要一套程序来解决非稳态偏微分方程组的数值
• 资源名：MATLAB求解偏微分方程（扩散方程）有限差分法 源程序代码_偏微分方程_有限差分法_matlab 资源类型：matlab项目全套源码 源码说明： 全部项目源码都是经过测试校正后百分百成功运行的，如果您下载后不能运行...
• 构造偏微分方程差分格式，对其弱形式进行网格剖分，用matlab对其进行求解
• set(0,'defaultfigurecolor','w') h=0.1; N=30; dt=0.0001; M=10000; A=dt/(h^2); U=zeros([N+1,M+1]); Space=0:h:(N)*h; for k=1:M+1 U(1,k) = 0.0; U(N+1,k) = 0.0; end for i =1:N U(i,1)=4*(i-1)*h*(3-... U
set(0,'defaultfigurecolor','w')
h=0.1;
N=30;
dt=0.0001;
M=10000;
A=dt/(h^2);
U=zeros([N+1,M+1]);
Space=0:h:(N)*h;

for k=1:M+1
U(1,k) = 0.0;
U(N+1,k) = 0.0;
end

for i =1:N
U(i,1)=4*(i-1)*h*(3-(i-1)*h);
end

for k=1:M
for i=2:N
U(i,k+1)=A*U(i+1,k)+(1-2*A)*U(i,k)+A*U(i-1,k);
end
end
figure;
plot(Space,U(:,1), 'g-')
hold on;
plot(Space,U(:,3000), 'b-')
hold on;
plot(Space,U(:,6000), 'k-')
hold on;
plot(Space,U(:,9000),'r-')
hold on;
plot(Space,U(:,10000),'y-')
ylabel('u(x,t)')
xlabel('x')
xlim([0,3])
ylim([-2,10])
legend('t=0','t=0.3','t=0.6','t=0.9','t=1')
figure;
contourf(U)
colorbar



展开全文
• pde1dM在单个空间变量和时间中求解偏微分方程 (PDE) 系统。 输入主要与 MATLAB 函数pdepe兼容。 许多pdepe示例只需少量更改即可与pde1dM一起使用。 与pde1dM相比， pdepe的主要增强是pde1dM允许任意数量的常微分方程...
• MATLAB程序分享求解偏微分方程扩散方程有限差分法-MATLAB求解偏微分方程（扩散方程）有限差分法 源程序代码.rar 程序代码见附件，拿资料请顺便顶个贴~~ 如果下载有问题，请加我 qq 1530497909，给你在线传
• matlab使用有限元方法求解偏微分方程
• Matlab求解微分方程(组)及偏微分方程(组)
• 资源名：MATLAB求解偏微分方程（扩散方程）有限差分法 源程序代码.rar 资源类型：matlab项目全套源码 源码说明： 全部项目源码都是经过测试校正后百分百成功运行的，如果您下载后不能运行可联系我进行指导或者更换。...
• 介绍了应用最为广泛的椭圆型、双曲型、抛物型偏微分方程的数值解法，而且还详细编程实现了每种方程的多种常见数值解法。 附件使用MATLAB编程来实现这些算法。
• 资源名：MATLAB求解偏微分方程（扩散方程）有限差分法 源程序代码.zip 资源类型：matlab项目全套源码 源码说明： 全部项目源码都是经过测试校正后百分百成功运行的，如果您下载后不能运行可联系我进行指导或者更换。...
• 本篇将介绍用matlab求解微分方程的数值和解析，并非是一种完整的模型，仅仅是一些算法。由于数学原理过于复杂，故不探究背后的数学原理，仅将matlab求解的相关函数加以记录。所有代码均可跑通。 1.Matlab求常...

本篇将介绍用matlab求解常微分方程的数值解和解析解，并非是一种完整的模型，仅仅是一些算法。由于数学原理过于复杂，故不探究背后的数学原理，仅将matlab求解的相关函数加以记录。所有代码均可跑通。

# 1.Matlab求常微分方程的数值解

## 1.1非刚性常微分方程的数值解法：

功能函数：ode45，ode23，ode113
例：用RK方法（四阶龙格—库塔方法）求解方程
f=-2y+2x^2+2*x

matlab程序：

//doty.m
function f=doty(x,y)
f=-2*y+2*x^2+2*x;
end

//main.m
[x,y]=ode45('doty',[0,0.5],1)


注：[0,0.5]表示求解区间；1为初值列向量

## 1.2刚性常微分方程的数值解法

功能函数：如ode15s，ode23s，ode23t， ode23tb
使用方法与非刚性类似

# 2.Matlab求常微分方程的解析解

## 2.1求常微分方程的通解

syms x y
diff_equ='x^2+y+(x-2*y)*Dy=0'
dsolve(diff_equ,'x')


注：'x’代表x为自变量，D代表求导

## 2.2求常微分方程的初边值问题

syms x y
diff_equ='D3y-D2y=x'
dsolve(diff_equ,'y(1)=8,Dy(1)=7,D2y(2)=4','x')


## 2.3求常微分方程组

equ1='D2f+3*g=sin(x)';
equ2='Dg+Df=cos(x)';
[general_f,general_g]=dsolve(equ1,equ2,'x')
[f,g]=dsolve(equ1,equ2,'Df(2)=0,f(3)=3,g(5)=1','x')


# 3.Matlab求解偏微分方程

%(1)问题定义
g='circleg'; %单位圆
b='circleb1'; %边界上为零条件
c=1;a=0;f=1; %（2）产生初始的三角形网格
[p,e,t]=initmesh(g); %（3）迭代直至得到误差允许范围内的合格解
error=[]; err=1;
while err > 0.01,
[p,e,t]=refinemesh(g,p,e,t);
u=assempde(b,p,e,t,c,a,f); %求得数值解
exact=(1-p(1,:).^2-p(2,:).^2)/4;
err=norm(u-exact',inf);
error=[error err];
end %结果显示
subplot(2,2,1),pdemesh(p,e,t);
subplot(2,2,2),pdesurf(p,t,u)
subplot(2,2,3),pdesurf(p,t,u-exact')


# 4.Matlab pdetool工具箱求解偏微分方程

对于一般的区域，任意边界条件的偏微分方程，我们可以利用Matlab中pdetool提供的偏微分方程用户图形界面解法。pdetool提供的用户图形界面解法的使用步骤如下：
（i）在Matlab命令窗口运行pdetool，出现PDE Toolbox界面。
（ii）用鼠标点一下工具栏上的“PDE"按钮，在弹出的对话框中定义偏微分方程。
（iii）用鼠标点一下工具栏上的区域按钮，在下面的坐标系中画出偏微分方程的大致定解区域。
（iv）双击（iii）中画出的大致区域，在弹出的对话框中精确定位定解区域。
（v）用鼠标点一下工具栏上的边界按钮“ ”，画出区域的边界。
（vi）双击坐标系中的区域边界，定义偏微分方程的边界条件。
（vii）用鼠标点工具栏上的剖分按钮，对求解区域进行剖分。
（viii）如果求抛物型或双曲型方程的数值解，还需要通过“solve”菜单下的“parameters…”选项设置初值条件。
（ix）用鼠标点一下工具栏上的“＝”按钮，就画出偏微分方程数值解的图形。通过“solve”菜单下的“Export Solution…”选项可以把数值解u输出到Matlab的工作间。
（x）如要画出数值解的三维图形，需要设置“plot"菜单下的“parameters…”选项。

展开全文
• MATLAB算法-求解微分方程数值和解析.ppt
• (ODE)、偏微分方程 (PDE) 和涉及微分方程的特征值问题。 与传统的 ODE 求解方法相比，在目标函数足够平滑的情况下，谱方法自然具有收敛速度超快的优势。 有关光谱方法的更多详细信息，请查看 。 它列出了用于理解谱...
• 偏微分方程的有限差分法的 MATLAB 和 Python 实现 这些代码实现了有限差分法的数值方法来求解 Heat PDE 和 Black-Scholes PDE。 具体而言，Black-Scholes PDE 的代码旨在为普通期权定价，例如欧洲和美国的看涨和看跌...
• 我们开发了一个新的变分框架来解决偏微分方程 (PDE)，该偏微分方程控制连续时间随机系统的联合概率密度函数的流动。 这是 Fokker-Planck-Kolmogorov PDE 的原型求解器，用于模拟专门的一维和二维系统的密度传播。 ...

...

matlab 订阅