还剩25页未读,继续阅读
本资源只提供10页预览,全部文档请下载后查看!喜欢就下载吧,查找使用更方便
文本内容:
第十六章偏微分方程的数值解法科学研究和工程技术中的许多问题可建立偏微分方程的数学模型包含多个自变量的微分方程称为偏微分方程partialdifferentialequation,简称PDE偏微分方程问题,其求解是十分困难的除少数特殊情况外,绝大多数情况均难以求出精确解因此,近似解法就显得更为重要本章仅介绍求解各类典型偏微分方程定解问题的差分方法
16.1几类偏微分方程的定解问题一个偏微分方程的表示通常如下
16.
1.1式中,是常数,称为拟线性quasilinear数通常,存在3种拟线性方程双曲型hyperbolic方程;抛物线型parabolic方程;椭圆型ellliptic方程
16.
1.2双曲型方程最简单形式为一阶双曲型方程
16.
1.2物理中常见的一维振动与波动问题可用二阶波动方程
16.
1.3描述,它是双曲型方程的典型形式方程的初值问题为
16.
1.4边界条件一般有三类,最简单的初边值问题为
16.
1.
516.
1.3抛物型方程其最简单的形式为一维热传导方程
16.
1.8方程可以有两种不同类型的定解问题1初值问题
16.
1.62初边值问题
16.
1.7其中,,为已知函数,且满足连接条件
16.
1.8边界条件为第一类边界条件第二类和第三类边界条件为
16.
1.9其中,当时,为第二类边界条件,否则称为第三类边界条件
16.
1.4椭圆型方程其最典型、最简单的形式是泊松Poisson方程
16.
1.10特别地,当时,即为拉普拉斯Lapla__方程,又称为调和方程
16.
1.11Poisson方程的第一边值问题为
16.
1.12其中为以为边界的有界区域,为分段光滑曲线,称为定解区域,,分别为,上的已知连续函数第二类和第三类边界条件可统一表示为
16.
1.13其中为边界的外法线方向当时为第二类边界条件,时为第三类边界条件
16.2差分方法的基本概念差分方法又称为有限差分方法或网格法,是求偏微分方程定解问题的数值解中应用最广泛的方法之一它的基本思想是先对求解区域作网格剖分,将自变量的连续变化区域用有限离散点网格点集代替;将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替;通过用网格点上函数的差商代替导数,将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组称为差分格式如果差分格式有解,且当网格无限变小时其解收敛于原微分方程定解问题的解,则差分格式的解就作为原问题的近似解数值解因此,用差分方法求偏微分方程定解问题一般需要解决以下问题1选取网格;2对微分方程及定解条件选择差分近似,列出差分格式;3求解差分格式;4讨论差分格式解对于微分方程解的收敛性及误差估计下面,用一个简单的例子来说明用差分方法求解偏微分方程问题的一般过程及差分方法的基本概念设有一阶双曲型方程初值问题
16.
2.1选取网格图
16.
2.1差分示意图首先对定解区域作网格剖分,最简单常用一种网格是用两族分别平行于轴与轴的等距直线
16.
2.2将分成许多小矩形区域这些直线称为网格线,其交点称为网格点,也称为节点,和分别称作方向和方向的步长这种网格称为矩形网格1对微分方程及定解条件选择差分近似,列出差分格式如果用向前差商表示一阶偏导数,即
16.
2.3其中方程
16.
2.4在节点处可表示为
16.
2.5其中由于当足够小时,在式中略去,就得到一个与方程相近似的差分方程
16.
2.6此处,可看作是问题的解在节点处的近似值同初值条件
16.
2.7结合,就得到求问题的数值解的差分格式式
16.
2.8称为差分方程的截断误差如果一个差分方程的截断误差为,则称差分方程对是阶精度,对是阶精度的显然,截断误差的阶数越大,差分方程对微分方程的逼近越好若网格步长趋于0时,差分方程的截断误差也趋于0,则称差分方程与相应的微分方程是相容的这是用差分方法求解偏微分方程问题的必要条件如果当网格步长趋于0时,差分格式的解收敛到相应微分方程定解问题的解,则称这种差分格式是收敛的
16.3双曲型方程的差分解法
16.
3.1一阶双曲型方程的差分格式考虑一阶双曲型方程的初值问题
16.
3.1将平面剖分成矩形网格,取方向步长为方向步长为,网格线为为简便,记以不同的差商近似偏导数,可以得到方程的不同的差分近似
16.
3.
216.
3.
316.
3.4截断误差分别为,与结合离散化的初始条件,可以得到几种简单的差分格式
16.
3.
516.
3.
616.
3.7其中如果已知第层节点上的值,按上面三种格式就可求出第层上的值因此,这三种格式都是显式格式如果对采用向后差商,采用向前差商,则方程可化成
16.
3.8相应的差分格式为
16.
3.9此差分格式是一种隐式格式,必须通过解方程组才能由第层节点上的值求出第层节点上的值例
16.
3.1对初值问题其中,用差分格式求其数值解,取解记,由初始条件按差分格式计算公式为计算结果略如果用差分格式求解,计算公式为计算结果略与准确解比较知,按前一个差分格式所求得的数值解不收敛到初值问题的解,而后一个差分格式的解收敛到原问题的解
16.
3.2波动方程的差分格式对二阶波动方程
16.
3.1如果令,,则方程可化成一阶线性双曲型方程组
16.
3.2记,则方程组可表成矩阵形式
16.
3.3矩阵有两个不同的特征值,故存在非奇异矩阵,使得
16.
3.4作变换,方程组可化为
16.
3.5方程组由二个__的一阶双曲型方程联立而成因此本节主要讨论一阶双曲型方程的差分解法下面给出如下波动方程和边界条件的差分格式
16.
3.
616.
3.7将矩形划分成个小矩形,长宽分别为,,形成一个网格,如图
16.
3.1所示图
16.
3.1网格图可通过差分方程的方法计算近似值在连续行内,网格点的真实值为求和的中心差分为
16.
3.
816.
3.9在每一行的网格间距是均匀的,而且同时,它在每一列也是均匀的,,接下来,将和去掉,用
16.
3.8和
16.
3.9中的近似,并按顺序代入
16.
3.6式,可得到差分方程
16.
3.10可用它来近似
16.
3.6式为方便起见,可将代入上式,可得
16.
3.11设行和的近似值已知,可用上式求网格的行
16.
3.12式中,根据上式左边的4个已知值可得到近似值,如图
16.
3.2所示图
16.
3.2波动方程的空格样板用上式时,必须注意,如果计算的某个阶段带来的误差最终会越来越小,则方法是稳定的为了保证上式的稳定性,必然使还存在其他一些差分方程方法以,称为隐格式法,它们更难实现,但对无限制
16.
3.3差分方法求解波动方程的__TLAB程序求解区间,以
16.
3.7为边界条件的波动方程的差分方法程序**********************************************************functionU=finediffgabcnm%Input-f=ux0asastringf%-g=utx0asastringg%-aandbrightendpointsof[0a]and[0b]%-ctheconstantinthew__eequation%-nandmnumberofgridpointsover[0a]and[0b]%Output-Usolution__trix;____ogoustoTable
10.1%InitializeparametersandUh=a/n-1;k=b/m-1;r=c*k/h;r2=r^2;r22=r^2/2;s1=1-r^2;s2=2-2*r^2;U=zerosnm;%Computfirstandsecondrowsfori=2:n-1Ui1=fevalfh*i-1;Ui2=s1*fevalfh*i-1+k*fevalgh*i-
1...+r22*fevalfh*i+fevalfh*i-2;end%Computere__iningrowsofUforj=3:mfori=2:n-1Uij=s2*Uij-1+r2*Ui-1j-1+Ui+1j-1-Uij-2;endendU=U;*******************************************************************
16.4抛物型方程的差分解法以一维热传导方程
16.
4.1为基本模型讨论适用于抛物型方程定解问题的几种差分格式
16.
4.1差分格式的建立首先对平面进行网格剖分分别取为方向与方向的步长,用两族平行直线,,将平面剖分成矩形网格,节点为为简便,记,,,,,,
16.
4.2微分方程的差分近似在网格内点处,对分别采用向前、向后及中心差商公式
16.
4.2一维热传导方程可分别表示为
16.
4.
316.
4.
416.
4.5由此得到一维热传导方程的不同差分近似
16.
4.
616.
4.
716.
4.8上述差分方程所用到的节点各不相同其截断误差分别为,和因此,它们都与一维热传导方程相容如果将式中的用代替,则可得到又一种差分近似
16.
4.9差分方程用到四个节点由Taylor公式容易得出
16.
4.10故其的截断误差为因而不是对任意的,此差分方程都能逼近热传导方程,仅当时,才成立综上可知,用不同的差商公式可以得到微分方程的不同的差分近似构造差分格式的关键在于使其具有相容性、收敛性和稳定性前面三个方程都具有相容性,而此方程则要在一定条件下才具有相容性
16.
4.3初、边值条件的处理为用差分方法求解定解问题初值问题
16.
4.11初边值问题
16.
4.12还需对定解条件进行离散化对初始条件及第一类边界条件,可直接得到,,式中对第
二、三类边界条件
16.
4.13需用差分近似下面介绍两种较简单的处理方法在左边界处用向前差商近似偏导数,在右边界处用向后差商近似,即,
16.
4.14则得边界条件的差分近似为,
16.
4.15其截断误差为2用中心差商近似即,
16.
4.16则得边界条件的差分近似为,
16.
4.17其截断误差为误差的阶数提高了,但出现定解区域外的节点和,这就需要将解拓展到定解区域外可以通过用内节点上的值插值求出和,也可以假定热传导方程在边界上也成立,将差分方程扩展到边界节点上,由此消去和
16.
4.4几种常用的差分格式以热传导方程的初边值问题
16.
4.18为例给出几种常用的差分格式1古典显式格式令,则
16.
4.19可改写成将其与初始条件及第一类边界条件,,结合,我们得到求解此问题的一种差分格式
16.
4.20由于第0层上节点处的值已知,由此即可算出在第一层上节点处的近似值重复使用此式,可以逐层计算出所有的,因此此差分格式称为典显式格式又因式中只出现相邻两个时间层的节点,故此式是二层显式格式2古典隐式格式将式
16.
4.21整理并与初始条件及第一类边界条件式联立,得差分格式如下
16.
4.22其中虽然第0层上的值仍为已知,但不能由上式直接计算以上各层节点上的值,必须通过解下列线性方程组
16.
4.23式中,才能由计算,故此差分格式称为古典隐式格式此方程组是三对角方程组,且系数矩阵严格对角占优,故解存在唯一3Richardson格式Richardson格式是将式整理后与初始条件及第一类边界条件式联立其计算公式为
6.
4.24这种差分格式中所涉及的节点出现在三层上,故为三层显式格式Richardson格式是一种完全不稳定的差分格式,因此它在实际计算中是不能采用的4杜福特-弗兰克尔DoFort-Frankel格式DoFort-Frankel格式也是三层显式格式,它是由式与初始条件及第一类边界条件式结合得到的具体形式如下
16.
4.25用这种格式求解时,除了第0层上的值由初值条件得到,必须先用二层格式求出第1层上的值,然后再按上式逐层计算5六点隐式格式对二阶中心差商公式如果用在与处的二阶中心差商的平均值近似在处的值,即同时在点处的值也用中心差商近似,即这样又得到热传导方程的一种差分近似
16.
4.26其截断误差为,将上式与初始条件及第一类边界条件式联立并整理,得差分格式
14.
5.27此格式涉及到六个节点,它又是隐式格式,故称为六点隐式格式与古典隐式格式类似,用六点格式由第层的值计算第层的值时,需求解三对角方程组
16.
4.28式中,此方程组的系数矩阵严格对角占优,故仍可用追赶法求解例
16.
4.1用古典显式格式求初边值问题的数值解,取解这里,,由格式可得到将初值代入上式,即可算出将边界条件,及上述结果代入又可求得如此逐层计算,得全部节点上的数值解为
16.
4.5前向差分法求解热传导方程的__TLAB程序设,其中,,而且,,其中,求解区间内的近似解*****************************************************functionU=forwdiffc1c2abcnm%Input-f=ux0asastringf%-c1=u0tandc2=uat%-aandbrightendpointsof[0a]and[0b]%-ctheconstantintheheatequation%-nandmnumberofgridpointsover[0a]and[0b]%Output-Usolution__trix;____ogoustoTable
10.4%InitializeparametersandUh=a/n-1;k=b/m-1;r=c^2*k/h^2;s=1-2*r;U=zerosnm;%BoundaryconditionsU11:m=c1;Un1:m=c2;%GeneratefirstrowU2:n-11=fevalfh:h:n-2*h;%Generatere__iningrowsofUforj=2:mfori=2:n-1Uij=s*Uij-1+r*Ui-1j-1+Ui+1j-1;endendU=U;*****************************************************
16.
4.6Crank-Nicholson求解热传导方程的__TLAB程序设,其中,,而且,,其中,求解区间内的近似解*****************************************************functionU=crnichfc1c2abcnm%Input-f=ux0asastringf%-c1=u0tandc2=uat%-aandbrightendpointsof[0a]and[0b]%-ctheconstantintheheatequation%-nandmnumberofgridpointsover[0a]and[0b]%Output-Usolution__trix;____ogoustoTable
10.5%InitializeparametersandUh=a/n-1;k=b/m-1;r=c^2*k/h^2;s1=2+2/r;s2=2/r-2;U=zerosnm;%BoundaryconditionsU11:m=c1;Un1:m=c2;%GeneratefirstrowU2:n-11=fevalfh:h:n-2*h;%Formthediagonalandoff-diagonalelemntsofAand%theconstantvectorBandsolvetridiagonalsystemAX=BVd11:n=s1*ones1n;Vd1=1;Vdn=1;Va=-ones1n-1;Van-1=0;Vc=-ones1n-1;Vc1=0;Vb1=c1;Vbn=c2;forj=2:mfori=2:n-1Vbi=Ui-1j-1+Ui+1j-1+s2*Uij-1;endX=trisysVaVdVcVb;U1:nj=X;endU=U*****************************************************
16.5椭圆型方程第一边值问题的差分解法本节以Poisson方程为基本模型讨论第一边值问题的差分方法
16.
5.1差分格式的建立考虑Poisson方程的第一边值问题
16.
5.1图
16.
5.1取分别为方向和方向的步长,如图所示,以两族平行线,将定解区域剖分成矩形网格节点的全体记为定解区域内部的节点称为内点,记内点集为边界与网格线的交点称为边界点,边界点全体记为与节点沿方向或方向只差一个步长的点和称为节点的相邻节点如果一个内点的四个相邻节点均属于,称为正则内点,正内点的全体记为,至少有一个相邻节点不属于的内点称为非正则内点,非正则内点的全体记为问题是要求出第一边值问题在全体内点上的数值解为简便,记,,对正则内点,由二阶中心差商公式
16.
5.2Poisson方程在点处可表示为
16.
5.3其中
16.
5.4为其截断误差表示式,略去,即得与方程相近似的差分方程
16.
5.5式中方程的个数等于正则内点的个数,而未知数则除了包含正则内点处解的近似值外,还包含一些非正则内点处的近似值,因而方程个数少于未知数个数在非正则内点处Poisson方程的差分近似不能按上式给出,需要利用边界条件得到边界条件的处理可以有各种方案,下面介绍较简单的两种1直接转移用最接近非正则内点的边界点上的值作为该点上值的近似,这就是边界条件的直接转移例如,点为非正则内点,其最接近的边界点为点,则有
16.
3.6上式可以看作是用零次插值得到非正则内点处的近似值,容易求出,其截断误差为将上式代入,方程个数即与未知数个数相等2线性插值这种方案是通过用同一条网格线上与点相邻的边界点与内点作线性插值得到非正则内点处值的近似由点与的线性插值确定的近似值,得
16.
5.7其中,其截断误差为将其与方程相近似的差分方程联立,得到方程个数与未知数个数相等的方程组,求解此方程组可得Poisson方程第一边值问题的数值解上面所给出的差分格式称为五点菱形格式,
16.
5.8实际计算时经常取,此时五点菱形格式可化为
16.
5.9简记为
16.
5.10其中例
16.
5.1用五点菱形格式求解拉普拉斯Lapla__方程第一边值问题其中取解网格中有四个内点,均为正则内点由五点菱形格式,得方程组代入边界条件其解为,,,当时,对利用点,,构造的差分格式,称为五点矩形格式,简记为
16.
5.11其中,其截断误差为
16.
5.12五点菱形格式与矩形格式的截断误差均为,称它们具有二阶精度如果用更多的点构造差分格式,其截断误差的阶数可以提高,如利用菱形格式及矩形格式所涉及的所有节点构造出的九点格式就是具有四阶精度的差分格式
16.
5.2用Dirichlet法求解Lapla__方程的__TLAB程序求解区间内的近似解而且满足条件,其中,,而且,,其中,而且存在整数和,使得*****************************************************functionU=dirichf1f2f3f4abhtol__x1%Input-f1f2f3f4areboundaryfunctionsinputasstrings%-aandbrightendpointsof[0a]and[0b]%-hstepsize%-tolisthetoleran__%Output-Usolution__trix;____ogoustoTable
10.6%InitializeparametersandUn=fixa/h+1;m=fixb/h+1;__e=a*fevalf10+fevalf
20...+b*fevalf30+fevalf40/2*a+2*b;U=__e*onesnm;%BoundaryconditionsU11:m=fevalf30:h:m-1*h;Un1:m=fevalf40:h:m-1*h;U1:n1=fevalf10:h:n-1*h;U1:nm=fevalf20:h:n-1*h;U11=U12+U21/2;U1m=U1m-1+U2m/2;Un1=Un-11+Un2/2;Unm=Un-1m+Unm-1/2;%SORparameterw=4/2+sqrt4-cospi/n-1+cospi/m-1^2;%Refineapproxi__tionsandsweepoperatorthroughoutthegriderr=1;cnt=0;whileerrtolcnt=__x1err=0;forj=2:m-1fori=2:n-1relx=w*Uij+1+Uij-1+Ui+1j+Ui-1j-4*Uij/4;Uij=Uij+relx;iferr=absrelxerr=absrelx;endendendcnt=cnt+1;endU=flipudU;*****************************************************。