偏微分方程的MATLAB数值解法

打造最优秀、专业和权威的Matlab技术交流平台! http://www.matlabsky.cn 偏微分方程的 MATLAB求解精讲© 作者：dynamic 时间：2008.12.10 版权：All Rights Reserved By www.matlabsky.cn ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★ Matlab Sky联盟----打造最优秀、专业和权威的Matlab技术交流平台! 网址：http://www.matlabsky.cn /com/org/net 邮箱：matlabsky@gmail.com QQ群：23830382 40510634 16233891(满了) 44851559(满了) 论坛拥有 40多个专业版块，内容涉及资料下载、视频教学、数学建模、数学运算、程序设计、GUI开发、simulink 仿真、统计概率、拟合优化、扩展编程、算法研究、控制系统、信号通信、图像处理、经济金融、生物化学、航 空航天、人工智能、汽车设计、机械自动化、毕业设计等几十个方面! 请相信我们：1.拥有绝对优秀的技术人员，热情的版主，严谨负责的管理团队 2.免费提供技术交流和在线解答 ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★ 打造最优秀、专业和权威的Matlab技术交流平台! http://www.matlabsky.cn MATLAB求解微分/偏微分方程，一直是一个头大的问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 ，两个字，“难过”，由于 MATLAB对 LaTeX的支持有 限，所有方程必须化成MATLAB可接受的标准形式，不支持像其他三个数学软件那样直接傻瓜式输入，这个真 把人给累坏了！ 不抱怨了，还是言归正传，回归我们今天的主体吧！ MATLAB提供了两种方法解决 PDE问题，一是 pdepe()函数，它可以求解一般的 PDEs，据用较大的通用性，但 只支持命令行形式调用。二是 PDE工具箱，可以求解特殊 PDE问题，PDEtool有较大的局限性，比如只能求解 二阶 PDE问题，并且不能解决偏微分方程组，但是它提供了 GUI界面，从繁杂的编程中解脱出来了，同时还可 以通过 File->Save As直接生成M代码 一、一般偏微分方程组(PDEs)的MATLAB求解.............................................................................................................3 1、pdepe函数 说明 关于失联党员情况说明岗位说明总经理岗位说明书会计岗位说明书行政主管岗位说明书 ......................................................................................................................................................3 2、实例讲解 ................................................................................................................................................................4 二、PDEtool求解特殊 PDE问题 ......................................................................................................................................6 1、典型偏微分方程的描述.........................................................................................................................................6 （1）椭圆型 ........................................................................................................................................................6 （2）抛物线型.....................................................................................................................................................6 （3）双曲线型.....................................................................................................................................................6 （4）特征值型.....................................................................................................................................................7 2、偏微分方程边界条件的描述.................................................................................................................................8 （1）Dirichlet条件 .............................................................................................................................................8 （2）Neumann条件 ............................................................................................................................................8 3、求解实例 ................................................................................................................................................................9 Ivan 高亮 Ivan 高亮 打造最优秀、专业和权威的Matlab技术交流平台! http://www.matlabsky.cn 一、一般偏微分方程组(PDEs)的MATLAB求解 1、pdepe函数说明 MATLAB语言提供了 pdepe()函数，可以直接求解一般偏微分方程(组)，它的调用格式为 sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t) 【输入参数】 @pdefun：是 PDE的问题描述函数，它必须换成下面的标准形式 ( , , ) [ ( , , , )] ( , , , ) ( )m mu u u uc x t x x f x t u s x t u x t x x x -¶ ¶ ¶ ¶ ¶= + ¶ ¶ ¶ ¶ ¶ 式1 这样，PDE就可以编写下面的入口函数 [c,f,s]=pdefun(x,t,u,du) m,x,t就是对应于(式 1)中相关参数，du是 u的一阶导数，由给定的输入变量即可表示出出 c,f,s这三个函数 @pdebc：是 PDE的边界条件描述函数，必须先化为下面的形式 ( , , ) ( , , ).* ( , , , ) 0up x t u q x t u f x t u x ¶ + = ¶ 于是边值条件可以编写下面函数描述为 [pa,qa,pb,qb]=pdebc(x,t,u,du) 其中 a表示下边界，b表示下边界 @pdeic：是 PDE的初值条件，必须化为下面的形式 0 0( , )u x t u= 股我们使用下面的简单的函数来描述为 u0=pdeic(x) m,x,t：就是对应于(式 1)中相关参数 【输出参数】 sol：是一个三维数组，sol(:,:,i)表示 ui的解，换句话说 uk对应 x(i)和 t(j)时的解为 sol(i,j,k) 通过 sol，我们可以使用 pdeval()直接计算某个点的函数值 打造最优秀、专业和权威的Matlab技术交流平台! http://www.matlabsky.cn 2、实例讲解 试求解下面的偏微分 2 1 1 1 22 2 2 2 1 22 0.024 ( ) 0.17 ( ) u u F u u t x u u F u u t x ì¶ ¶ = - -ïï ¶ ¶ í ¶ ¶ï = - -ï ¶ ¶î 其中， 5.73 11.46( ) x xF x e e-= - ，且满足初始条件 1 2( ,0) 1, ( ,0) 0u x u x= = 及边界条件 1 2 2 1(0, ) 0, (0, ) 0, (1, ) 1, (1, ) 0 u ut u t u t t x x ¶ ¶ = = = = ¶ ¶ 【解】 （1）对照给出的偏微分方程，根据标注形式，则原方程可以改写为 1 1 1 2 2 2 1 2 0.024 ( )1 .* 1 ( )0.17 u u F u ux u u F u ut t x ¶é ù ê ú - -é ù é ùé ù ¶ ¶ ¶= +ê úê ú ê úê ú ¶ -¶ ¶ë û ë û ê ú ë û ê ú¶ë û 可见 m=0，且 1 1 2 2 1 2 0.024 ( )1 , , 1 ( )0.17 u F u uxc f s u F u u x ¶é ù ê ú - -é ùé ù ¶= = =ê ú ê úê ú ¶ -ë û ê ú ë û ê ú¶ë û %% 目标 PDE函数 function [c,f,s]=pdefun (x,t,u,du) c=[1;1]; f=[0.024*du(1);0.17*du(2)]; temp=u(1)-u(2); s=[-1;1].*(exp(5.73*temp)-exp(-11.46*temp)); （2）边界条件改写为 1 2 0 11 0 1 0 .* .* 0 0 0 00 u f f u -é ù é ùé ù é ù é ù é ù + = + =ê ú ê úê ú ê ú ê ú ê ú ë û ë û ë û ë ûë û ë û 下边界 上边界 %% 边界条件函数 function [pa,qa,pb,qb]=pdebc(xa,ua,xb,ub,t) %a表示下边界，b表示上边界 pa=[0;ua(2)]; 打造最优秀、专业和权威的Matlab技术交流平台! http://www.matlabsky.cn qa=[1;0]; pb=[ub(1)-1;0]; qb=[0;1]; （3）初值条件改写为 1 2 1 0 u u é ù é ù =ê ú ê ú ë ûë û %% 初值条件函数 function u0=pdeic(x) u0=[1;0]; （4）最后编写主调函数 clc x=0:0.05:1; t=0:0.05:2; m=0; sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t); figure('numbertitle','off','name','PDE Demo——by Matlabsky') subplot(211) surf(x,t,sol(:,:,1)) title('The Solution of u_1') xlabel('X') ylabel('T') zlabel('U') subplot(212) surf(x,t,sol(:,:,2)) title('The Solution of u_2') xlabel('X') ylabel('T') zlabel('U') 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 00.2 0.40.6 0.8 11.2 1.41.6 1.8 2 0 0.2 0.4 0.6 0.8 1 X The Solution of u1 T U 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 00.2 0.40.6 0.81 1.21.4 1.6 1.82 0 0.2 0.4 0.6 0.8 X The Solution of u2 T U 打造最优秀、专业和权威的Matlab技术交流平台! http://www.matlabsky.cn 二、PDEtool求解特殊 PDE问题 MATLAB的偏微分工具箱(PDE toolbox)可以比较规范的求解各种常见的二阶偏微分方程，但是惋惜的是只能求 解特殊二阶的 PDE问题，并且不支持偏微分方程组！ PDE toolbox支持命令行形式求解 PDE问题，但是要记住那些命令以及调用形式真的很累人，还好MATLAB提 供了 GUI可视交互界面 pdetool，在 pdetool中可以很方便的求解一个 PDE问题，并且可以帮我们直接生成M代 码(File->Save As)。 下面我们先了解下三个典型的二阶 PDE，然后介绍 pdetool，至于命令行我们就免了，真的很累人，如果的确需 要的话，那就让Matalb直接生成就好了。 1、典型偏微分方程的描述 （1）椭圆型偏微分方程的一般形式为 2 2 2 2 2 2 1 2 ( ) * ( , ) * * ( , ) , , n div c u a u f x t c u a u f x t x x x c a f - Ñ + = ¶ ¶ ¶ - + + + + = ¶ ¶ ¶ L 即 ( ) 其中 为给定的函数或者常数 （2）抛物线型偏微分方程的一般形式 2 2 2 2 2 2 1 2 * ( ) * ( , ) * * * ( , ) , , , n ud div c u a u f x t t ud c u a u f x t t x x x d c a f ¶ - Ñ + = ¶ ¶ ¶ ¶ ¶ - + + + + = ¶ ¶ ¶ ¶ L 即 ( ) 其中 必须是常数 （3）双曲线型偏微分方程的一般形式 打造最优秀、专业和权威的Matlab技术交流平台! http://www.matlabsky.cn 2 2 2 2 2 2 2 2 2 2 1 2 * ( ) * ( , ) * * * ( , ) , , , n ud div c u a u f x t t ud c u a u f x t t x x x d c a f ¶ - Ñ + = ¶ ¶ ¶ ¶ ¶ - + + + + = ¶ ¶ ¶ ¶ L 即 ( ) 其中 必须是常数 （4）特征值型偏微分方程的一般形式，注意它是（1）的变形，不能算独立的一类 2 2 2 2 2 2 1 2 ( ) * * * * * * * n div c u a u d u c u a u d u x x x l l - Ñ + = ¶ ¶ ¶ - + + + + = ¶ ¶ ¶ L 即 ( ) 从上面可以看出，三类典型二阶偏微分方程的区别在于 u对 t的导数阶次。椭圆型 PDEs中，c、a、d和 f可以 是给定的函数或者常数，但是其它两类必须都是常数。 MATLAB是采用有限元的方法求解各种 PDE。MATLAB为我们提供一个 pdetool的交互界面，可以求解二元偏 微分 u(x1,x2)(注意只能求解二元)。 方程的参数由 a、c、d和 f确定，求解域由图形确定，求解域确定好后，需要对求解域进行栅格化(这个是自动)。 打造最优秀、专业和权威的Matlab技术交流平台! http://www.matlabsky.cn 2、偏微分方程边界条件的描述 一般在 PDE中边界条件包括 Dirichlet(狄利克莱)条件和 Neumann(纽曼)条件： （1）Dirichlet条件 一般描述为 ( , , , )* | ( , , , )u uh x t u u r x t u x x¶W ¶ ¶ = ¶W ¶ ¶ ，其中 表示求解域的边界 假设在边界上满足该方程，则只需给出 r和 h即可，它们可以是常数也可以是给定的函数 （2）Neumann条件 一般描述为 [ ( ) * ] | uc u q u g u n n¶W ¶ ¶ Ñ + = ¶ ¶ ，其中 表示 的法向偏导数 通过下面的操作调出边界条件设置，注意在这之前一定要使用【区域边界】按钮制定边界 打造最优秀、专业和权威的Matlab技术交流平台! http://www.matlabsky.cn 3、求解实例 试求解双曲线型偏微分方程 2 2 2 2 2 2 2 10 u u u u t x y ¶ ¶ ¶ - - + = ¶ ¶ ¶ 求解域 s为 2 2 2 2 : 9 1 4 16 ( ) x y x y + £ + £ = U I s1 s2: s s1 s2)-(s1 s2 边界条件为 构成求解域的边界值都为 5 【解】 由给定的 PDE，可以的 d=1,c=1,a=2,f=10，再说一次，对于抛物线和双曲线型偏微分方程 4个系数必须是常数， 否则MATLAB无能为力 打造最优秀、专业和权威的Matlab技术交流平台! http://www.matlabsky.cn step1:点击工具栏的【PDE】按钮，如下输入 PDE的参数，注意选择 Hyperbolic step2:绘制求解域 对坐标轴的操作可以在【Options】主菜单中操作，包括设置网格、坐标系范围等 (1)【Options】->Axis Limits设置如下 其它设置如下 (2)点击工具栏上的第三个按钮【绘制椭圆】，任意绘制一个椭圆，双击椭圆，设置如下 打造最优秀、专业和权威的Matlab技术交流平台! http://www.matlabsky.cn 重复上面的操作，参数如下 于是得到 打造最优秀、专业和权威的Matlab技术交流平台! http://www.matlabsky.cn (3)在 set formula中如下输入，“+”表示求并集，“-”表示求差集，注意没有直接求交接的操作符 step3:边界条件和初值条件 初值条件可以通过【Solve】->【Parameters…】设置 边值条件设置如下 (1)点击工具栏的第 6个按钮【区域边界】，显示如下 (2)【Boundary】->【Remove All Subdomain Borders】移除所有子域的边界，将得到所有子域合并成一个求解域 (3) 【Boundary】->【Secify Boundary Conditons…】设置边界如下，注意我们这里只有 Dirichlet条件 打造最优秀、专业和权威的Matlab技术交流平台! http://www.matlabsky.cn step4:生成使用有限元方法求解方程所需的栅格 点击工具栏的第 8/9个按钮，对求解域生成栅格，多次点击可以在原来基础上继续细化栅格，直到自己觉得满意 为止，当然可以通过【Mesh】主菜单进行精确控制 step5:求解方程 点解工具栏的第 10个按钮“=”【求解方程】 step6:求解结果绘图 点击第 11个按钮【绘制图形】，里面的选项很丰富，可以绘制等高线等好多，甚至播放动画，具体大家可以自 己慢慢摸索 打造最优秀、专业和权威的Matlab技术交流平台! http://www.matlabsky.cn 动画播放设置： (1)【Solve】->【Parameters】设置合适的时间向量 Time (2)【Plot】->【Parameters】选中【Animation】，点击后面的【Options】，设置播放速度和次数，比如 6fps表示 每秒 6帧 (3)【Plot】->【Export Movie…】输入动画保存的变量名，比如M (4)在 Command Windows中直接输入 movie(M)即可播放 (5)使用 movie2ve(M,’demo.avi’)命令可以将动画保存为 avi文件