零基础使用 MATLAB 求解偏微分方程（建议收藏） |
您所在的位置：网站首页 › 有限元法matlab代码 › 零基础使用 MATLAB 求解偏微分方程（建议收藏） |
零基础使用 MATLAB 求解偏微分方程（建议收藏）
|
零基础使用 MATLAB 求解偏微分方程（建议收藏）
文章目录
零基础使用 MATLAB 求解偏微分方程（建议收藏）偏微分开源工具介绍PDE 工具箱函数汇总介绍0 基础：GUI 界面操作示例问题工具箱求解导出为代码形式代码导出相关数据
0.1 基础：编程调用 PDE 工具箱PDE 工具箱的局限性
偏微分开源工具介绍
百分之九十以上的重要的工程和数学科学研究，和偏微分方程都脱不开关系。在所有的偏微分方程中，百分之九十九都是没有解析解的。没有解析解怎么办，我们只能通过有限元或者有限差分等方法，求解偏微分方程数值解。如果您有一些代码基础，建议参考我的几篇有限经典博文，简单问题可在此基础上进行修改。 有限元方法入门：有限元方法简单的一维算例 有限元方法入门：有限元方法简单的二维算例（三角形剖分） 有限元方法入门：有限元方法简单的二维算例（矩形剖分） 对于做工程的朋友，不会偏微分方程数值解，怎么办？没关系，我推荐一些求解各类偏微分方程的容易入门的开源的软件包和工具，它们是： Free FEM++（足够傻瓜又不失自由度，强烈建议做工程的朋友可以学习一下）FEniCS(C++/Python, 开始于芝加哥大学和查尔姆斯理工大学，升级版是 FEniCSX)PETSC (C/Python, 美国阿贡国家实验室)deal.II (C++, 开始于德国海德堡大学，如果你非要学一个 C++ 有限元工具，而又不知道选哪个的话，可以看看这个)MFEM (C++, 美国劳伦斯利弗莫尔国家实验室)PHG (C, 张林波, 中国科学院)AFEPACK (C++, 李若, 北京大学）FEALPy（Python，魏华祎，湘潭大学）IFEM (MATLAB, 陈龙, UCI)NGSolve（C++/Python，Christoph 等）PHOEBESolver（ Fortran and C/C++，宁夏大学，葛永斌）GCGE（C/MATLAB，中国科学院，谢和虎，特征值求解），这是代数特征值求解的工具包，不是 PDE 的。……当然，商业有限元软件 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 文件，并打开。
当前目录下，保存如下代码为 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 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 nargin1 % 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') PDE 工具箱的局限性当然，MATLAB 的 PDE 工具箱也有很多局限性，比如一些 CFD 的问题没办法弄，比如说柱坐标的问题，据我所知，他目前只支持直角坐标。再比如说经典的 NS 方程，可能没法搞，因为 PDE 工具箱没有为我们提供指定对流项的方法，所以我们可能无法使用它来模拟 Navier-Stokes 方程。 如果不将所有非线性项放在右侧，就无法将 Navier-Stokes 方程转换为工具箱接受的这种形式。 因此，一般来说，PDE 工具箱可能不足以进行流体力学模拟。 我们无法将对流项放尽到当中的。 那 MATLAB 应该拿流体问题怎么办？不要慌，你可以考虑一下 FEATool Multiphysics。 |
【本文地址】
今日新闻 |
推荐新闻 |