精华内容
下载资源
问答
  • 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求解偏微分方程(扩散方程)有限差分法,处理偏微分方程
  • 是一系列已实现的神经网络,用于求解偏微分方程。 该项目属于作者在米兰理工大学航空工程中的硕士论文“求解偏微分方程的现代非确定性方法:机器学习应用于纳维-斯托克斯方程”。 为了完整起见,这里报告了硕士论文...
  • matlab解偏微分方程

    2018-09-17 21:51:49
    代码为2018年全国数学建模竞赛a题第一问matlab源程序,可动态生成三层隔热服距离与温度的关系图,以及三层隔热服的温度分布图。主要内容建立在一维非稳态热传导以及偏微分方程的解法程序
  • 求解偏微分的代码薛定谔偏微分方程:数值和模拟。 作者:Alejandro Galván、Daniel Quiles 和 Bernat Ramis - UPC 数学和工程物理专业的学生。 这项工作的目的是研究数值方法对求解薛定谔偏微分方程的适用性。 ...
  • 求解偏微分的代码数值偏微分方程 包含有限差分方法对偏微分方程的实际实现的存储库。 这个存储库包含我在约克完成的一个模块中的一些项目,以及我自己的一些工作。 在这个存储库中编写的代码应该能够在 MATLAB 和 ...
  • 零基础使用 MATLAB 求解偏微分方程(建议收藏)

    万次阅读 多人点赞 2021-10-02 01:26:20
    零基础使用 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 几何图形
    %   addVertex               - 在几何区域边界上添加一个顶点
    %   specifyCoefficients     - 指定区域和或者子区域的 PDE 系数
    %   applyBoundaryCondition  - 给几何区域施加边界条件
    %   setInitialConditions    - 设定 PDE 初始条件
    %   generateMesh            - 从几何生成一个网格
    %   solvepde                - 求解 PDE
    %   solvepdeeig             - 求解 PDE 特征值问题
    %   assembleFEMatrices      - 组装中间的有限元矩阵
    %   createPDEResults        - 创建一个用于后处理的结果对象
    %   pdegplot                - 绘制 PDE 几何表示
    %   pdemesh                 - 绘制 PDE 网格
    %   pdeplot                 - 绘制二维 PDE 网格和结果
    %   pdeplot3D               - 绘制 3D PDE 网格和结果
    %   interpolateSolution     - 在指定的空间位置插入解
    %   evaluateGradient        - 在指定的空间位置评估解的梯度
    %   evaluateCGradient       - 评估 PDE 解的通量
    % 
    % 使用热模型解决以传导为主的传热问题
    %   thermalProperties       - 为热模型指定材料的热属性
    %   internalHeatSource      - 指定热模型的内部热源
    %   thermalBC               - 指定热模型的边界条件
    %   thermalIC               - 设置热模型的初始条件或初始猜测
    %   solve                   - 求解热模型中指定的传热问题
    %   interpolateTemperature	- 在任意空间位置的热结果中插入温度
    %   evaluateTemperatureGradient - 评估热解在任意空间位置的温度梯度
    %   evaluateHeatFlux        - 在节点或任意空间位置评估热解的热通量
    %   evaluateHeatRate        - 评估垂直于指定边界的综合热流率
    % 
    % 使用结构模型解决静态、模态和瞬态线性弹性问题
    %   structuralProperties    - 为模型分配结构材料属性
    %   structuralBodyLoad      - 将体载荷应用于结构模型
    %   structuralBoundaryLoad  - 在几何边界上施加结构载荷
    %   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        - 逆离散正弦变换
    %   pdeadgsc    - 使用相对容差标准挑选坏三角形
    %   pdeadworst  - 选择相对于最差值的坏三角形
    %   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       - 热解及其派生量
    %   SteadyStateThermalResults - 稳态热模型解及其派生量
    %   TransientThermalResults   - 瞬态热模型解及其派生量
    %   StructuralMaterialAssignment - 区域或子域上的结构材料属性分配
    %   BodyLoadAssignment   - 结构分析模型的体载荷分配
    %   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 型区域,如下图所示 Δu=λuuΩ=0Ω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,'Hgrad',1.3);
    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',...
    char('0','1548','10','pdeadworst',...
    '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', [], ...
          'MenuBar', 'none', ...
          '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'.
    %
    % See also MATQDLG, MATQUEUE.
    
    % 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,'Hgrad',1.3);
    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',...
        char('0','1548','10','pdeadworst',...
        '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',...
                'Tag','PDEMeshMenu'),'UserData');
            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');
            h=findobj(get(pde_fig,'Children'),'flat','Tag','PDEBoundMenu');
            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:
            h=findobj(get(pde_fig,'Children'),'flat','Tag','PDEMeshMenu');
            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:
            params=get(findobj(get(pde_fig,'Children'),'Tag','PDEPDEMenu'),...
                '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:
            u=get(findobj(get(pde_fig,'Children'),'flat','Tag','PDEPlotMenu'),...
                'UserData');
            l=get(findobj(get(pde_fig,'Children'),'flat','Tag','winmenu'),...
                '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 Δu=λ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;
    radius = 2;
    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解偏微分方程

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

    千次阅读 2021-04-15 09:33:55
    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
    使用方法与非刚性类似

    1.3高阶微分方程的解法

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    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偏微分方程快速上手:使用pde有限元工具箱求解二维偏微分方程
    偏微分方程的数值解(六): 偏微分方程的 pdetool 解法

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

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 3,106
精华内容 1,242
关键字:

matlab求解偏微分方程

matlab 订阅