偏微分方程式之求解.doc
《偏微分方程式之求解.doc》由会员分享,可在线阅读,更多相关《偏微分方程式之求解.doc(37页珍藏版)》请在咨信网上搜索。
1、个人收集整理 勿做商业用途第六章 偏微分方程式之求解在化工的领域中,有不少程序之动态是由以偏微分方程式(Partial differential equation;PDE)所描述的,例如热与质量在空间中的传递等。这些用以描述实际问题的PDE,除非具有某些特定的方程式型态及条件,否则甚难以手算的方式找出解析解.而在数值求解方面,最常被采用的方法为有限差分法(finite difference)何有限元素法(finite element)。然对于某些不熟悉数值分析及程序编写的化工人而言,欲充分了解以偏微分方程式所描述之系统动态是相当不容易的,更遑论进一步的设计与分析了. 值得庆幸的是,MATLAB
2、的环境中提供了一个求解PDE问题的工具箱,让使用者得以利用简单的指令或图形接口工具输入欲解的PDE,并求解。使得PDE之数值解在弹指之间完成,使用者不在为数值法所苦恼,轻松掌握偏微分方程式系统的动态,并可进一步进行后续之设计工作。 本章将以循渐进的方式,介绍PDE工具箱及其用法,并以数个典型的化工范例进行 示范,期能使初学者很快熟悉PDE工具箱,并使用它来设计与分析以偏微方方程式所描述的程序系统。6.1 偏微分方程式之分类偏微分方程式可根据其阶数(order),线性或非线性型态,以及边界条件进行分类。6。1.1依阶数的分类 偏微分方程式是以偏微分项中之最高次偏微分来定义其阶数,例如: 一阶偏微
3、分方程式: 二阶偏微分方程式: 三阶偏微分方程式: 6.1.2 依非线性程度分类 偏微分方程式亦可以其线性或非线性情况,区分为线性(linear),似线性(quasilinear),以及非线性三类。例如,以下之二阶偏微分方程式(Constantinides and Mostoufi,1999) 可依系数之情况,进行如下表之归类类别 情况线性 系数为定值,或仅为(x,y)函数似线性 系数为依变数(dependent variable)u或其比方程式中之偏微分低阶之偏微分项的函数,如非线性 系数中,具有与原方程式之偏微分同阶数之变数,如 另外,对于线性二阶偏微分方程式,可进一步将其分类为椭圆型(e
4、lliptic),拋物线型(parabolic),以及双曲线型(hyperbolic)。具体上来说,此类偏微分方程式二阶线性之一般式为 系数是定值或为u的函数。若g=0,则上式为其次是偏微分方程式。式子( )之分类及代表性例子,请见下表方程式类别 判断式 代表性范例椭圆型 Laplace方程式,拋物线型 Poisson方程式, 热传导或扩散方程式 双曲线型 波动方程式 注:二维系统之运操作数之定义为6.1。3 起始条件和边界条件的分类 为了能获得偏微分方程式之解答,其起始条件和边界条件可依其特性区分为三类.现以一维之动态热传递方程式(拋物线型偏微分方程式) 为例,进一步说明如何区分这些边界条件
5、及起始条件(Constantinides and Mostoufi ,1999). (i) 第一类:Dirichlet Condiction 若依变量(T)本身,在某个独立变量值时,被指定,则此条件称为Dirichlet Condiction,亦称为essential边界条件。下图为一典型的Dirichlet条件示意图图6.1 平板Dirichlet Condiction示意图 由图中很清楚的显示,该平板之边界条件为 此边界条件依定义,即为Dirichlet Condiction。同时,若再起始时,各处之温度分布可以位置之函数表示,即 此亦属Dirichlet型之边界条件. (ii) 第二类:
6、Neumann condition Neumann condition系指依变量之变化率之边界条件为定值,抑或独立变量之函数之情况.例如 或 Neumann型边界条件,亦称为natural boundary condition。 (iii) 第三类:Robbins condition若依变量之变化率之边界条件,为自身之函数(非独立变量之函数)时,被称为Robbins condition。例如, 上式之边界条件,当发生在固液相间之传递上,亦即热流通量(heat flux)正比于固液两端之温差,其示意图如下: (iv) Cauchy condition Cauchy condition系指问题中同
7、时存在Dirichlet和Neumann边界条件.例如下图 即 “Dirichlet condition “Neumann condition” 6.2 MATLAB PDE工具箱6。2.1 MATLAB PDE解答器 MATLAB提供了一个指令pdepe,用以解以下的PDE方程式 其中时间介于之间,而位置则介于有限区域之间.m值表示问题之对称性,其可为0,1或2,以表示平板(slab),圆柱(cylindrical)或球体(spherical)之对秤情况.因而,如果如果m0,则a必等于b,也就是说其具有圆柱或球体之对称关系。同时,式中一项为流通量(flux),而为来源(source)项.为偏
8、微分方程式之对角线系数矩阵.若某一对角线元素为0,则表示该相应偏微分方程式为椭圆型偏微分方程式,若为正值(不为0),则为拋物线型偏微分方程式。请注意c之对角线元素必不全为0。偏微分方程式之起始值可表示为 而边界条件之形式为 其中x为两端点位置,即a或b。 用以解含上述起始值及边界值条件的偏微分方程式( )之指令pdepe的用法如下: 其中 为问题之对称参数 为独立变量x之网取点(mesh)位置向量,即,其中(起点),(终点)。 为独立变量t(时间)之向量,即,其中为起始时间,为终点时间。 为使用者提供之pde函数文件。其档案之格式如下: 亦即,使用者仅需提供偏微分方程式中之系数向量.,和均为行
9、(column)向量,而向量即为矩阵之对角线元素. 提供解u之起始值,其格式为值得注意的是u为行向量。 使用者提供之边界条件档,格式如下:其中ul和ur分别表示左边界和右边界之u 的近似解。输出变量中,pl和ql分别表示左边界p和q之行向量,而pr和qr则为右边界p和q之行向量. 为解答输出。为多维的输出向量,为的输出,即。另,元素表示在和时之答案. 为解答器之相关解法参数。详细请见odeset之使用法。 注:1. MATLAB PDE解答器pdepe之算法,主要事将原一组的椭圆型和拋物线型偏微分方程式转化为一组常微分方程式。此转换的过程系基于使用者所指定之mesh点,以二阶空间离散化(spa
10、tial discretization)技术为之(Keel and Berzins,1990),然后以ode15s的指令求解。采用ode15s的ode解法,主要是因为在离散化的过程中,椭圆型偏微分方程式为被转为一组代数式,而拋物线型的偏微分方程式则被转为一组联立的微分方程式。因而,原偏微分方程式被离散化后,遂变成一组同时伴有微分方程式与代数式之微分代数方程式系统,故以ode15s使可顺利求解。 2。 x之取点(mesh)位置影响解的精确度甚鉅,若pdepe解答器回应出has difficulty finding consistent initial considition之讯息时,使用者可进一
11、步将mesh点取密一点,即增加mesh点数。另外,若状态u再某些特定点上有较快速之变动时,亦需将此处之点取密集些,以增加精确度。值得注意的是pdepe并不会自动做并不会自动做xmesh的自动取点,使用者必须观察解的特性,自行作曲点的动作.一般而言,所取之点数至少需大于3以上。3. tspan之选取主要是基于使用者对那些特定时间之状态有兴趣而选定.而间距(step size)之控制由程序自动完成。4。 若要获得特定位置及时间下之解,可配合以pdeval的指令。跟pdepe指令之格式如下: 其中 代表问题之对称性。=0,代表平板;=1,圆柱体;=2表示球体。其意义同pdepe中之自变量。 使用者在
12、pdepe中所指定之输出点之位置向量。 即。也就是说其为pdepe输出中第i个输出变数在各点位置处,时间固定为下之解。 为所欲内插输出点位置向量。此为使用者重新指定之位置向量。 为基于所指定位置,固定时间下之相对应输出. 为相对应之输出值。 ref. Keel,R.D. and M。 Berzins,”A Method for the Spatial Discritization of Parabolic Equations in One Space Variable”,SIAM J。 Sci。 and Sat. Comput。,Vol。11,pp.132,1990。 以下吾人将以数个范例,详
13、细说明pdepe之用法。范例6.1 试解以下之偏微分方程式 其中,且满足以下之条件限制式(i)起始值条件IC:(ii)边界条件 BC1: BC2:注:本问题之解析解为 解答: 以下吾人将分述求解之步骤与过程。当以下之各步骤完成后,可进一步将其汇整为一主程序ex6_1.m,然后求解之。 步骤1 将欲解之偏微分方程式改写成如式子( )之标准式。 此即 和 步骤2 编写偏微分方程式之系数向量档function c,f,s=ex6_1pdefun(x,t,u,DuDx)c=pi2;f=dudx;s=0; 步骤3 编写起始值条件function u0=ex6_1ic(x)u0=sin(pix); 步骤4
14、 编写边界条件。在编写之前,先将边界条件改写成标准形式,如式( ),找出相对之和的函数,然后写成边界条件档,例如,原边界条件可写成 BC1: BC2: 即 故而,边界条件档可编写成 function pl,ql,pr,qr=ex6_1bc(xl,ul,xr,ur,t)pl=ul;ql=0;pr=piexp(t);qr=1; 步骤5 取点。例如 x=linspace(0,1,20); x取20点t=linspace(0,2,5); %时间取5点输出 步骤6 利用pdepe求解 m=0; 依步骤1之结果sol=pdepe(m,ex6_1pdefun,ex6_1ic,ex6_1bc,x,t); 步骤
15、7 显示结果 u=sol(:,:,1);surf(x,t,u)title(pde数值解)xlabel(位置)ylabel(时间 )zlabel(u)若要显示特定点上之解,可进一步指定x 或t的位置,以便绘图。例如,吾人欲了解时间为2(终点)时,各位置下之解,可输入以下指令(利用pdeval指令): figure(2); 绘成图2M=length(t); %取终点时间xout=linspace(0,1,100); 输出点位置uout,dudx=pdeval(m,x,u(M,:),xout);plot(xout,uout); %绘图title(时间为2时,各处之解)xlabel(x)ylabel(
16、u) 综合以上各步骤,吾人可写成一程序求解范例6.1。其参考程序如下:ex6_1.mfunction ex6_1范例6_1m=0;x=linspace(0,1,20); %xmesht=linspace(0,2,20); tspan%以pde求解sol=pdepe(m,ex6_1pdefun,ex6_1ic,ex6_1bc,x,t);u=sol(:,:,1); %取出答案%绘图输出figure(1)surf(x,t,u)title(pde数值解)xlabel(位置x)ylabel(时间t )zlabel(数值解u)%与解析解做比较%figure(2)surf(x,t,exp(-t)*sin(p
17、i*x);title(解析解)xlabel(位置x)ylabel(时间t )zlabel(数值解u)t=tf=2时各位置之解figure(3)M=length(t); 取终点时间xout=linspace(0,1,100); %输出点位置uout,dudx=pdeval(m,x,u(M,:),xout);plot(xout,uout); 绘图title(时间为2时,各处之解)xlabel(x)ylabel(u)%pde函数文件function c,f,s=ex6_1pdefun(x,t,u,DuDx)c=pi2;f=DuDx;s=0;%起始条件档%function u0=ex6_1ic(x)u
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分 方程式 求解
1、咨信平台为文档C2C交易模式,即用户上传的文档直接被用户下载,收益归上传人(含作者)所有;本站仅是提供信息存储空间和展示预览,仅对用户上传内容的表现方式做保护处理,对上载内容不做任何修改或编辑。所展示的作品文档包括内容和图片全部来源于网络用户和作者上传投稿,我们不确定上传用户享有完全著作权,根据《信息网络传播权保护条例》,如果侵犯了您的版权、权益或隐私,请联系我们,核实后会尽快下架及时删除,并可随时和客服了解处理情况,尊重保护知识产权我们共同努力。
2、文档的总页数、文档格式和文档大小以系统显示为准(内容中显示的页数不一定正确),网站客服只以系统显示的页数、文件格式、文档大小作为仲裁依据,平台无法对文档的真实性、完整性、权威性、准确性、专业性及其观点立场做任何保证或承诺,下载前须认真查看,确认无误后再购买,务必慎重购买;若有违法违纪将进行移交司法处理,若涉侵权平台将进行基本处罚并下架。
3、本站所有内容均由用户上传,付费前请自行鉴别,如您付费,意味着您已接受本站规则且自行承担风险,本站不进行额外附加服务,虚拟产品一经售出概不退款(未进行购买下载可退充值款),文档一经付费(服务费)、不意味着购买了该文档的版权,仅供个人/单位学习、研究之用,不得用于商业用途,未经授权,严禁复制、发行、汇编、翻译或者网络传播等,侵权必究。
4、如你看到网页展示的文档有www.zixin.com.cn水印,是因预览和防盗链等技术需要对页面进行转换压缩成图而已,我们并不对上传的文档进行任何编辑或修改,文档下载后都不会有水印标识(原文档上传前个别存留的除外),下载后原文更清晰;试题试卷类文档,如果标题没有明确说明有答案则都视为没有答案,请知晓;PPT和DOC文档可被视为“模板”,允许上传人保留章节、目录结构的情况下删减部份的内容;PDF文档不管是原文档转换或图片扫描而得,本站不作要求视为允许,下载前自行私信或留言给上传者【w****g】。
5、本文档所展示的图片、画像、字体、音乐的版权可能需版权方额外授权,请谨慎使用;网站提供的党政主题相关内容(国旗、国徽、党徽--等)目的在于配合国家政策宣传,仅限个人学习分享使用,禁止用于任何广告和商用目的。
6、文档遇到问题,请及时私信或留言给本站上传会员【w****g】,需本站解决可联系【 微信客服】、【 QQ客服】,若有其他问题请点击或扫码反馈【 服务填表】;文档侵犯商业秘密、侵犯著作权、侵犯人身权等,请点击“【 版权申诉】”(推荐),意见反馈和侵权处理邮箱:1219186828@qq.com;也可以拔打客服电话:4008-655-100;投诉/维权电话:4009-655-100。