关闭

关闭

关闭

封号提示

内容

首页 第20章 偏微分方程的数值解.pdf

第20章 偏微分方程的数值解.pdf

第20章 偏微分方程的数值解.pdf

上传者: Gingerjin 2012-07-26 评分 0 0 0 0 0 0 暂无简介 简介 举报

简介:本文档为《第20章 偏微分方程的数值解pdf》,可适用于工程科技领域,主题内容包含第二十章偏微分方程的数值解自然科学与工程技术中种种运动发展过程与平衡现象各自遵守一定的规律。这些规律的定量表述一般地呈现为关于含有未知函数及其导数的符等。

第二十章偏微分方程的数值解自然科学与工程技术中种种运动发展过程与平衡现象各自遵守一定的规律。这些规律的定量表述一般地呈现为关于含有未知函数及其导数的方程。我们将只含有未知多元函数及其偏导数的方程称之为偏微分方程。方程中出现的未知函数偏导数的最高阶数称为偏微分方程的阶。如果方程中对于未知函数和它的所有偏导数都是线性的这样的方程称为线性偏微分方程否则称它为非线性偏微分方程。初始条件和边界条件称为定解条件未附加定解条件的偏微分方程称为泛定方程。对于一个具体的问题定解条件与泛定方程总是同时提出。定解条件与泛定方程作为一个整体称为定解问题。偏微分方程的定解问题各种物理性质的定常(即不随时间变化)过程都可用椭圆型方程来描述。其最典型、最简单的形式是泊松(Poisson)方程),(yxfyuxuu==Δ()特别地当),(yxf时即为拉普拉斯(Laplace)方程又称为调和方程==Δyuxuu()带有稳定热源或内部无热源的稳定温度场的温度分布不可压缩流体的稳定无旋流动及静电场的电势等均满足这类方程。Poisson方程的第一边值问题为Ω=Γ=Ω=Γ),(|),(),(),(),(yxyxuyxyxfyuxuyxϕ()其中Ω为以Γ为边界的有界区域Γ为分段光滑曲线ΓΩU称为定解区域),(),,(yxyxfϕ分别为ΓΩ,上的已知连续函数。第二类和第三类边界条件可统一表示成),(),(yxunuyxϕα=Γ()其中n为边界Γ的外法线方向。当=α时为第二类边界条件α时为第三类边界条件。在研究热传导过程气体扩散现象及电磁场的传播等随时间变化的非定常物理问题时常常会遇到抛物型方程。其最简单的形式为一维热传导方程)(>=axuatu()方程()可以有两种不同类型的定解问题:初值问题(也称为Cauchy问题)<<=<<>=xxxuxtxuatu)(),(,ϕ()初边值问题===<<<<=TttgtlutgtuxxulxTtxuatu),(),(),(),()(),(,ϕ()其中)(),(),(tgtgxϕ为已知函数且满足连接条件)()(),()(glg==ϕϕ问题()中的边界条件)(),(),(),(tgtlutgtu==称为第一类边界条件。第二类和第三类边界条件为TttgutxuTttgutxulxx====),()(),()(λλ()其中)(,)(ttλλ。当)()(=ttλλ时为第二类边界条件否则称为第三类边界条件。双曲型方程的最简单形式为一阶双曲型方程=xuatu()物理中常见的一维振动与波动问题可用二阶波动方程xuatu=()描述它是双曲型方程的典型形式。方程()的初值问题为<<=<<=<<>==xxtuxxxuxtxuatut)()(),(,φϕ()边界条件一般也有三类最简单的初边值问题为====<<>==Tttgtlutgtulxxtuxxulxtxuatut)(),(),(),()(),(),(,φϕ如果偏微分方程定解问题的解存在唯一且连续依赖于定解数据(即出现在方程和定解条件中的已知函数)则此定解问题是适定的。可以证明上面所举各种定解问题都是适定的。偏微分方程的差分解法差分方法又称为有限差分方法或网格法是求偏微分方程定解问题的数值解中应用最广泛的方法之一。它的基本思想是:先对求解区域作网格剖分将自变量的连续变化区域用有限离散点(网格点)集代替将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替通过用网格点上函数的差商代替导数将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组(称为差分格式)。如果差分格式有解且当网格无限变小时其解收敛于原微分方程定解问题的解则差分格式的解就作为原问题的近似解(数值解)。因此用差分方法求偏微分方程定解问题一般需要解决以下问题:(i)选取网格(ii)对微分方程及定解条件选择差分近似列出差分格式(iii)求解差分格式(iv)讨论差分格式解对于微分方程解的收敛性及误差估计。下面我们只对偏微分方程的差分解法作一简要的介绍。椭圆型方程第一边值问题的差分解法以Poisson方程()为基本模型讨论第一边值问题的差分方法。考虑Poisson方程的第一边值问题()Ω=Γ=Ω=Γ),(|),(),(),(),(yxyxuyxyxfyuxuyxϕ取τ,h分别为x方向和y方向的步长以两族平行线τjyykhxxjk====,),,,,(L=jk将定解区域剖分成矩形网格。节点的全体记为},,,|),{(为整数jijykhxyxRjkjkτ===。定解区域内部的节点称为内点记内点集ΩIR为τhΩ。边界Γ与网格线的交点称为边界点边界点全体记为τhΓ。与节点),(jkyx沿x方向或y方向只差一个步长的点),(jkyx和),(jkyx称为节点),(jkyx的相邻节点。如果一个内点的四个相邻节点均属于ΓΩU称为正则内点正则内点的全体记为)(Ω至少有一个相邻节点不属于ΓΩU的内点称为非正则内点非正则内点的全体记为)(Ω。我们的问题是要求出问题()在全体内点上的数值解。为简便记记),(),,(),(),,(),(,jkjkjkjkyxffyxujkuyxjk===。对正则内点)(),(Ωjk由二阶中心差商公式)(),(),(),(),(hOhjkujkujkuxujk=)(),(),(),(),(ττOjkujkujkuyujk=Poisson方程()在点),(jk处可表示为)(),(),(),(),(),(),(,ττ=hOfjkujkujkuhjkujkujkujk()在式()中略去)(τhO即得与方程()相近似的差分方程jkjkjkjkjkjkjkfuuuhuuu,,,,,,,=τ()式()中方程的个数等于正则内点的个数而未知数jku,则除了包含正则内点处解u的近似值还包含一些非正则内点处u的近似值因而方程个数少于未知数个数。在非正则内点处Poisson方程的差分近似不能按式()给出需要利用边界条件得到。边界条件的处理可以有各种方案下面介绍较简单的两种。(i)直接转移(ii)线性插值由式()所给出的差分格式称为五点菱形格式实际计算时经常取τ=h此时五点菱形格式可化为jkjkjkjkjkjkfuuuuuh,,,,,,)(=()简记为jkjkfuh,,=()其中jkjkjkjkjkjkuuuuuu,,,,,,=。求解差分方程组最常用的方法是同步迭代法同步迭代法是最简单的迭代方式。除边界节点外区域内节点的初始值是任意取定的。例用五点菱形格式求解Laplace方程第一边值问题Ω=Γ=Ω=Γ)lg(|),(),(),(yxyxuyxyuxuyx其中},|),{(=Ωyxyx。取==τh。当τ=h时利用点),(),,(),,(jkjkjk构造的差分格式jkjkjkjkjkjkfuuuuuh,,,,,,)(=()称为五点矩形格式简记为hٱjkjkfu,,=()其中ٱjkjkjkjkjkjkuuuuuu,,,,,,=。抛物型方程的差分解法以一维热传导方程())(>=atuatu为基本模型讨论适用于抛物型方程定解问题的几种差分格式。首先对xt平面进行网格剖分。分别取τ,h为x方向与t方向的步长用两族平行直线),,,(L===kkhxxk),,,(L===jjttjτ将xt平面剖分成矩形网格节点为),,,,,,,)(,(LL==jktxjk。为简便起见记),(),(jkyxjk=),(),(jkyxujku=)(kkxϕϕ=)(jjtgg=)(jjtgg=)(jjtλλ=)(jjtλλ=。微分方程的差分近似在网格内点),(jk处对tu分别采用向前、向后及中心差商公式对xu采用二阶中心差商公式一维热传导方程()可分别表示为)(),(),(),(),(),(hOhjkujkujkuajkujku=ττ)(),(),(),(),(),(hOhjkujkujkuajkujku=ττ)(),(),(),(),(),(hOhjkujkujkuajkujku=ττ由此得到一维热传导方程的不同的差分近似,,,,,=huuuauujkjkjkjkjkτ(),,,,,=huuuauujkjkjkjkjkτ(),,,,,=huuuauujkjkjkjkjkτ()初、边值条件的处理为用差分方程求解定解问题()()等还需对定解条件进行离散化。对初始条件及第一类边界条件可直接得到kkkxuuϕ==),(,),,,,,(nkkLL==或()jjjnijjjgtluugtuu,,),(),(====),,,(=mjL()其中τTmhln==,。对第二、三类边界条件则需用差商近似。下面介绍两种较简单的处理方法。(i)在左边界)(=x处用向前差商近似偏导数xu在右边界)(lx=处用向后差商近似偏导数xu即)(),(),()(),(),(),(),(hOhjnujnuxuhOhjujuxujnj==),,,(mjL=即得边界条件()的差分近似为==jjnjjnjnjjjjjguhuuguhuu,,,,,,λλ),,,(mjL=()(ii)用中心差商近似xu即)(),(),()(),(),(),(),(hOhjnujnuxuhOhjujuxujnj==),,,(mjL=则得边界条件的差分近似为==jjnjjnjnjjjjjguhuuguhuu,,,,,,λλ),,,(mjL=()这样处理边界条件误差的阶数提高了但式()中出现定解区域外的节点),(j和),(jn这就需要将解拓展到定解区域外。可以通过用内节点上的u值插值求出ju,和jnu,也可以假定热传导方程()在边界上也成立将差分方程扩展到边界节点上由此消去ju,和jnu,。几种常用的差分格式下面我们以热传导方程的初边值问题()为例给出几种常用的差分格式。(i)古典显式格式为便于计算令harτ=式()改写成以下形式jkjkjkjkruurruu,,,,)(=将式()与()()结合我们得到求解问题()的一种差分格式========),,,(,),,,(),,,,,,,()(,,,,,,,mjgugunkumjnkruurruujjnjjkkjkjkjkjkLLLLϕ()由于第层)(=j上节点处的u值已知)(,kkuϕ=由式()即可算出u在第一层)(=j上节点处的近似值,ku。重复使用式()可以逐层计算出各层节点的近似值。(ii)古典隐式格式将()整理并与式()()联立得差分格式如下========),,,(,),,,(),,,,,,,)((,,,,,,,,mjgugunkumjnkuuuruujjnjjkkjkjkjkjkjkLLLLϕ()其中harτ=。虽然第层上的u值仍为已知但不能由式()直接计算以上各层节点上的值jku,故差分格式()称为古典隐式格式。(iii)杜福特弗兰克尔(DoFortFrankel)格式DoFortFrankel格式是三层显式格式它是由式()与()()结合得到的。具体形式如下:========),,,(,),,,(),,,,,,,()(,,,,,,,mjgugunkumjnkurruurrujjnjjkkjkjkjkjkLLLLϕ()用这种格式求解时除了第层上的值,ku由初值条件()得到必须先用二层格式求出第层上的值,ku然后再按格式()逐层计算),,,(,mjujkL=。双曲型方程的差分解法对二阶波动方程()xuatu=如果令xuvtuv==,则方程()可化成一阶线性双曲型方程组==xvtvxvatv()记Tvvv),(=则方程组()可表成矩阵形式xvAxvatv==()矩阵A有两个不同的特征值a=λ故存在非奇异矩阵P使得Λ==aaPAP作变换TwwPvw),(==方程组()可化成xwtwΛ=()方程组()由两个独立的一阶双曲型方程联立而成。因此下面主要讨论一阶双曲型方程的差分解法。考虑一阶双曲型方程的初值问题<<=<<>=xxxuxtxuatu)(),(ϕ()与抛物型方程的讨论类似仍将xt平面剖分成矩形网格。取x方向步长为ht方向步长为τ网格线为),,,,(L===kkhxxk),,,(L===jjttjτ。为简便起见记),(),(jkyxjk=),(),(jkyxujku=)(kkxϕϕ=。以不同的差商近似偏导数可以得到方程()的不同的差分近似,,,,=huuauujkjkjkjkτ(),,,,=huuauujkjkjkjkτ(),,,,=huuauujkjkjkjkτ()结合离散化的初始条件可以得到几种简单的差分格式。一维状态空间的偏微分方程的MATLAB解法工具箱命令介绍MATLAB提供了一个指令pdepe用以解以下的PDE方程式),,,()),,,((),,,(xuutxsxuutxfxxxtuxuutxcmm=()其中时间介于fttt之间而位置x则介于,ba有限区域之间。m值表示问题的对称性其可为或分别表示平板(slab)圆柱(cylindrical)或球体(spherical)的情形。因而如果>m则a必等于b也就是说其具有圆柱或球体的对称关系。同时式中),,,(xuutxf一项为流通量(flux)而),,,(xuutxs为来源(source)项。),,,(xuutxc为偏微分方程的对角线系数矩阵。若某一对角线元素为则表示该偏微分方程为椭圆型偏微分方程若为正值(不为)则为拋物型偏微分方程。请注意c的对角线元素一定不全为。偏微分方程初始值可表示为)(),(xvtxu=()而边界条件为),,,(),(),,(=xuutxftxqutxp()其中x为两端点位置即a或b用以解含上述初始值及边界值条件的偏微分方程的MATLAB命令pdepe的用法如下:),,,,,,(optionstspanxmeshbcfunicfunpdepempdepesol=其中m为问题之对称参数xmesh为空间变量x的网格点(mesh)位置向量即,,,NxxxxmeshL=其中ax=(起点)bxN=(终点)。tspan为时间变量t的向量即,,,MttttspanL=其中t为起始时间Mt为终点时间。pdefun为使用者提供的pde函数文件。其函数格式如下:),,,(,,dudxutxpdefunsfc=亦即使用者仅需提供偏微分方程中的系数向量。cf和s均为行(column)向量而向量c即为矩阵c的对角线元素。icfun提供解u的起始值其格式为)(xicfunu=值得注意的是u为行向量。bcfun使用者提供的边界条件函数格式如下:()turxrulxlbcfunqlprqlpl,,,,,,,=其中ul和ur分别表示左边界)(bxl=和右边界)(axr=u的近似解。输出变量中pl和ql分别表示左边界p和q的行向量而pr和qr则为右边界p和q的行向量。sol为解答输出。sol为多维的输出向量):(:,isol为iu的输出即):,(:,isolui=。元素),,(),(ikjsolkjui=表示在)(jtspant=和)(kxmeshx=时iu之答案。options为求解器的相关解法参数。详细说明参见odeset的使用方法。注:MATLABPDE求解器pdepe的算法主要是将原来的椭圆型和拋物线型偏微分方程转化为一组常微分方程。此转换的过程是基于使用者所指定的mesh点以二阶空间离散化(spatialdiscretization)技术为之(KeelandBerzins,)然后以odes的指令求解。采用odes的ode解法主要是因为在离散化的过程中椭圆型偏微分方程被转化为一组代数方程而拋物线型的偏微分方程则被转化为一组联立的微分方程。因而原偏微分方程被离散化后变成一组同时伴有微分方程与代数方程的微分代数方程组故以odes便可顺利求解。x的取点(mesh)位置对解的精确度影响很大若pdepe求解器给出“…hasdifficultyfindingconsistentinitialconsidition”的讯息时使用者可进一步将mesh点取密一点即增加mesh点数。另外若状态u在某些特定点上有较快速的变动时亦需将此处的点取密集些以增加精确度。值得注意的是pdepe并不会自动做xmesh的自动取点使用者必须观察解的特性自行作取点的操作。一般而言所取的点数至少需大于以上。tspan的选取主要是基于使用者对那些特定时间的状态有兴趣而选定。而间距(stepsize)的控制由程序自动完成。若要获得特定位置及时间下的解可配合以pdeval命令。使用格式如下:),,,(,xoutuixmeshmpdevalduoutdxuout=其中m代表问题的对称性。m=表示平板m=表示圆柱体m=表示球体。其意义同pdepe中的自变量m。xmesh为使用者在pdepe中所指定的输出点位置向量。,,,Nxmeshxxx=L。ui即):,,(ijsol。也就是说其为pdepe输出中第i个输出ui在各点位置xmesh处时间固定为)(jtspantj=下的解。xout为所欲内插输出点位置向量。此为使用者重新指定的位置向量。uout为基于所指定位置xout固定时间ft下的相对应输出。duoutdx为相对应的dxdu输出值。refKeel,RDandMBerzins,“AMethodfortheSpatialDiscritizationofParabolicEquationsinOneSpaceVariable”SIAMJSciandSatComput,Vol,pp,以下将以数个例子详细说明pdepe的用法。求解一维偏微分方程例试解以下之偏微分方程式xutu=π其中x且满足以下之条件限制式(i)起始值条件IC:)sin(),(xxuπ=(ii)边界条件BC:),(=tuBC:),(=xtuetπ注:本问题的解析解为)sin(),(xetxutπ=解下面将叙述求解的步骤与过程。当完成以下各步骤后可进一步将其汇总为一主程序exm然后求解。步骤将欲求解的偏微分方程改写成如式的标准式。=xuxxxtuπ此即(),,,π=xuutxc()xuxuutxf=,,,(),,,=xuutxs和=m。步骤编写偏微分方程的系数向量函数。functionc,f,s=expdefun(x,t,u,dudx)c=pi^f=dudxs=步骤编写起始值条件。functionu=exic(x)u=sin(pi*x)步骤编写边界条件。在编写之前先将边界条件改写成标准形式如式()找出相对应的)(p和)(q函数然后写出MATLAB的边界条件函数例如原边界条件可写成BC:,),(),(==xtxtuBC:,),(==xtxuetπ即(,),,plutql==和,tpreqrπ==因而边界条件函数可编写成functionpl,ql,pr,qr=exbc(xl,ul,xr,ur,t)pl=ulql=pr=pi*exp(t)qr=步骤取点。例如x=linspace(,,)x取点t=linspace(,,)时间取点输出步骤利用pdepe求解。m=依步骤之结果sol=pdepe(m,expdefun,exic,exbc,x,t)步骤显示结果。u=sol(:,:,)surf(x,t,u)title('pde数值解')xlabel('位置')ylabel('时间')zlabel('u')若要显示特定点上的解可进一步指定x或t的位置以便绘图。例如欲了解时间为(终点)时各位置下的解可输入以下指令(利用pdeval指令):figure()绘成图M=length(t)取终点时间的下标xout=linspace(,,)输出点位置uout,dudx=pdeval(m,x,u(M,:),xout)plot(xout,uout)绘图title('时间为时,各位置下的解')xlabel('x')ylabel('u')综合以上各步骤可写成一个程序求解例。其参考程序如下:functionex************************************求解一维热传导偏微分方程的一个综合函数程序************************************m=x=linspace(,,)xmesht=linspace(,,)tspan************以pde求解************sol=pdepe(m,expdefun,exic,exbc,x,t)u=sol(:,:,)取出答案************绘图输出************figure()surf(x,t,u)title('pde数值解')xlabel('位置x')ylabel('时间t')zlabel('数值解u')*************与解析解做比较*************figure()surf(x,t,exp(t)'*sin(pi*x))title('解析解')xlabel('位置x')ylabel('时间t')zlabel('数值解u')*****************t=tf=时各位置之解*****************figure()M=length(t)取终点时间的下表xout=linspace(,,)输出点位置uout,dudx=pdeval(m,x,u(M,:),xout)plot(xout,uout)绘图title('时间为时,各位置下的解')xlabel('x')ylabel('u')******************pde函数******************functionc,f,s=expdefun(x,t,u,dudx)c=pi^f=dudxs=******************初始条件函数******************functionu=exic(x)u=sin(pi*x)******************边界条件函数******************functionpl,ql,pr,qr=exbc(xl,ul,xr,ur,t)pl=ulql=pr=pi*exp(t)qr=例试解以下联立的偏微分方程系统)(uuFxutu=)(uuFxutu=其中))(exp())(exp()(uuuuuuF=且x和t。此联立偏微分方程系统满足以下初边值条件。(i)初值条件),(=xu),(=xu(ii)边值条件),(=txu),(=tu),(=tu),(=txu解步骤:改写偏微分方程为标准式=•)()(*uuFuuFxuxuxuut因此=c=xuxuf=)()(uuFuuFs和=m。另外左边界条件(x=处)。写成=•xuxuu即=upl=ql同理右边界条件(x=处)为=•xuxuu即=upr=qr步骤:编写偏微分方程的系数向量函数。functionc,f,s=expdefun(x,t,u,dudx)c='f='*dudxy=u()u()F=exp(*y)exp(*y)s=FF'步骤:编写初始条件函数functionu=exic(x)u='步骤:编写边界条件函数functionpl,ql,pr,qr=exbc(xl,ul,xr,ur,t)pl=ul()'ql='pr=ur()'qr='步骤:取点。由于此问题的端点均受边界条件的限制且时间t很小时状态的变动很大(由多次求解后的经验得知)故在两端点处的点可稍微密集些。同时对于t小处亦可取密一些。例如x=t=以上几个主要步骤编写完成后事实上就可直接完成主程序来求解。此问题的参考程序如下:functionex***************************************求解一维偏微分方程组的一个综合函数程序***************************************m=x=t=*************************************利用pdepe求解*************************************sol=pdepe(m,expdefun,exic,exbc,x,t)u=sol(:,:,)第一个状态之数值解输出u=sol(:,:,)第二个状态之数值解输出*************************************绘图输出*************************************figure()surf(x,t,u)title('u之数值解')xlabel('x')ylabel('t')figure()surf(x,t,u)title('u之数值解')xlabel('x')ylabel('t')***************************************pde函数***************************************functionc,f,s=expdefun(x,t,u,dudx)c='f='*dudxy=u()u()F=exp(*y)exp(*y)s=FF'****************************************初始条件函数****************************************functionu=exic(x)u='****************************************边界条件函数****************************************functionpl,ql,pr,qr=exbc(xl,ul,xr,ur,t)pl=ul()'ql='pr=ur()'qr='化工应用实例例触煤反应装置内温度及转换率的分布以外部热交换式的管形固定层触煤反应装置进行苯加氢反应产生环己烷。此反应系统之质量平衡及热平衡方程式如下:)(=ΔprBApeGCHrrTrrTGCkLTρ=GyMrrfrrfuDLfarBAeρ其中T为温度()f为反应率L为轴向距离r为径向距离。此系统的边界条件为=L,)(rTT=,)(rff==r,==rfrTwrr=,)(wweTThrTk=wrr=,=rf此外式中之相关数据及操作条件如下:(i)反应速率式H苯氢环己烷图反应示意图)(CCBBHHBHBHAPKPKPKPPKkKr=其中P表示分压(atm)而速率参数为RRTKln=RRTKHln=RRTKBln=RRTKCln=上式中下标BH及C分别代表苯氢及环己烷。R为理想气体常数(calmolK)。(ii)操作条件及物性数据总压atmPt=反应管管径cmrw=壁温=wT质量速度hrmkgG=苯对氢之莫耳流量比=m反应管入口的苯之莫耳分率=y反应气体之平均分子量=avM触煤层密度mkgPB=流体平均比热molkgkcalCP=反应热molkgkcalHr=Δ整体传热系数Chrmkcalh=进料温度)(=T反应管管长cmL=流速hrmu=有效热传导系数ChrmkcalKe=壁境膜传热系

用户评论(0)

0/200

精彩专题

上传我的资料

每篇奖励 +2积分

资料评价:

/30
0下载券 下载 加入VIP, 送下载券

意见
反馈

立即扫码关注

爱问共享资料微信公众号

返回
顶部