差分法求解偏微分方程MAAB

南京理工大学课程考核论文课程名称：高等数值分析论文题目：有限差分法求解偏微分方程姓名：罗晨学号：成绩：有限差分法求解偏微分方程一、主要内容1.有限差分法求解偏微分方程，偏微分方程如一般形式的一维抛物线型方程：具体求解的偏微分方程如下：2.推导五种差分格式、截断误差并分析其稳定性；3.编写MATLAB程序实现五种差分格式对偏微分方程的求解及误差分析；4.结论及完成本次实验报告的感想。

二、推导几种差分格式的过程：有限差分法（finite-differencemethods ）是一种数值方法通过有限个微分方程近似求导从而寻求微分方程的近似解。

有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。

推导差分方程的过程中需要用到的泰勒展开公式如下：()2100000000()()()()()()()......()(())1!2!!n n n f x f x f x f x f x x x x x x x o x x n +'''=+-+-++-+-(2-1)求解区域的网格划分步长参数如下：11k k k kt t x x h τ++-=⎧⎨-=⎩(2-2) 2.1古典显格式2.1.1古典显格式的推导由泰勒展开公式将(,)u x t 对时间展开得2,(,)(,)()()(())i i k i k k k uu x t u x t t t o t t t∂=+-+-∂(2-3) 当1k t t +=时有21,112,(,)(,)()()(())(,)()()i k i k i k k k k k i k i k uu x t u x t t t o t t tuu x t o tττ+++∂=+-+-∂∂=+⋅+∂(2-4)得到对时间的一阶偏导数1,(,)(,)()=()i k i k i k u x t u x t uo t ττ+-∂+∂(2-5) 由泰勒展开公式将(,)u x t 对位置展开得223,,21(,)(,)()()()()(())2!k i k i k i i k i i u uu x t u x t x x x x o x x x x∂∂=+-+-+-∂∂(2-6)当11i i x x x x +-==和时，代入式(2-6)得2231,1,1122231,1,1121(,)(,)()()()()(())2!1(,)(,)()()()()(())2!i k i k i k i i i k i i i i i k i k i k i i i k i i i iu u u x t u x t x x x x o x x x xu u u x t u x t x x x x o x x x x ++++----⎧∂∂=+-+-+-⎪⎪∂∂⎨∂∂⎪=+-+-+-⎪∂∂⎩(2-7) 因为1k k x x h +-=，代入上式得2231,,22231,,21(,)(,)()()()2!1(,)(,)()()()2!i k i k i k i k i k i k i k i ku u u x t u x t h h o h x xu u u x t u x t h h o h x x +-⎧∂∂=+⋅+⋅+⎪⎪∂∂⎨∂∂⎪=-⋅+⋅+⎪∂∂⎩(2-8) 得到对位置的二阶偏导数2211,22(,)2(,)(,)()()i k i k i k i k u x t u x t u x t uo h x h+--+∂=+∂(2-9) 将式(2-5)、(2-9)代入一般形式的抛物线型偏微分方程得(2-10)为了方便我们可以将式(2-10)写成11122k kk k k k i i i i i i u u u u u f h ατ++-⎡⎤--+-=⎢⎥⎣⎦(2-11) ()11122k k k k k k i i i i i i u u uu u f hτατ++----+=(2-12)最后得到古典显格式的差分格式为()111(12)k k k k k i i i i i u ra u r u u f ατ++-=-+++(2-13)2r hτ=其中，古典显格式的差分格式的截断误差是2()o h τ+。

2.1.2古典显格式稳定性分析古典显格式(2-13)写成矩阵形式为 ()112k k k h h h u ra I raC u f τ+=-++⎡⎤⎣⎦(2-14)上面的C 矩阵的特征值是：2cos()1,2,......,1C j h j N λπ==-()()()212=122cos()121cos()14sin 1,2,......,12H j C ra ra ra ra j h ra j h j hra j N λλπππ=-+-+=--=-=-(2-15)使()1H ρ≤，即结论：当102ra ≤≤时，所以古典显格式是稳定的。

2.2古典隐格式2.2.1古典隐格式的推导 将1k t t -=代入式(2-3)得21,11(,)(,)()()(())j k j k j k k k k k uu x t u x t t t o t t t---∂=+-+-∂(2-16) 21,(,)(,)()()j k j k j k uu x t u x t o tττ-∂=-⋅+∂(2-17)得到对时间的一阶偏导数1,(,)(,)()=()j k j k j k u x t u x t u o t ττ--∂+∂(2-18) 将式(2-9)、(2-18)原方程得到(2-19)为了方便把(2-19)写成11122k k k k kj jj j j k j u u u u u f h ατ-+-⎡⎤--+-=⎢⎥⎢⎥⎣⎦(2-20) ()11122k k kk k kj jj j j j u u uu u f h τατ-+----+=(2-21)最后得到古典隐格式的差分格式为()111(12)k k k k k j j j jj ra u r u u u f ατ-+-+-+=+(2-22) 2r h τ=其中，古典隐格式的差分格式的截断误差是2()o h τ+。

2.2.2古典隐格式稳定性分析将古典隐格式(2-22)写成矩阵形式如下()1212()k k khh h ra I raC u u f r hττ++-=+=⎡⎤⎣⎦(2-23)误差传播方程()112k k h h ra I raC v v ++-=⎡⎤⎣⎦(2-24) 所以误差方程的系数矩阵为 使()1H ρ≤，显然1H j λ≤恒成立。

结论：对于0r ∀>，即任意网格比下，古典隐格式是绝对稳定的。

2.3Richardson 格式2.3.1Richardson 格式的推导 将11k k t t t t +-==和，代入式(2-3)得21,1121,11(,)(,)()()(())(,)(,)()()(())i k i k i k k k k k i k i k i k k k k ku u x t u x t t t o t t t u u x t u x t t t o t t t +++---∂⎧=+-+-⎪⎪∂⎨∂⎪=+-+-⎪∂⎩(2-25) 即21,21,(,)(,)()()(,)(,)()()i k i k i k i k i k i ku u x t u x t o t u u x t u x t o t ττττ+-∂⎧=+⋅+⎪⎪∂⎨∂⎪=-⋅+⎪∂⎩(2-26) 由此得到可得211,(,)(,)()()2i k i k i k u x t u x t uo t ττ++-∂=+∂(2-27) 将式(2-9)、(2-27)代入原方程得到下式 (2-28)为了方便可以把式(2-28)写成1111222k k k k k k i i i i i i u u u u u f h ατ+-+-⎡⎤--+-=⎢⎥⎣⎦(2-29) 即()111122k k kk k k i i i i i i u u uu u f h τατ+-+----+=(2-30)最后得到Richardson 显格式的差分格式为()1111222k k k k k k i i i i i i u r u u u u f ατ+-+-=-+++(2-31)2r hτ=其中，古典显格式的差分格式的截断误差是22()o h τ+。

2.3.2Richardson 稳定性分析将Richardson 显格式(2-31)写成如下矩阵形式 ()11222k k k k h h h h u r C I u u f ατ+-=-++(2-32)误差传播方程矩阵形式()1122k k k h h hkkh hv r C I v v v v α+-⎧=-+⎪⎨=⎪⎩(2-33) 再将上面的方程组写成矩阵形式112(2)0k k k hk ra C I I v ww I v ++-⎛⎫⎡⎤== ⎪⎢⎥⎣⎦⎝⎭(2-34) 系数矩阵的特征值是228sin 102j hra πλλ+-=(2-35) 解得特征值为1,2λ={}212,=4sin 12j h Max ra πλλ>(恒成立)(2-37) 结论：上式对任意的网比都恒成立，即Richardson 格式是绝对不稳定的。

4.Crank-Nicholson 格式3.4.1Crank-Nicholson 格式的推导 将12k t t +=代入式(2-9)得1112221112222,2+1,1+1+1(,)(,)()()(())(,)(,)()()(())j j k j k k k k k k j j k j k k k k k k u u x t u x t t t o t t t u u x t u x t t t o t t t +++++++∂⎧=+-+-⎪⎪∂⎨∂⎪=+-+-⎪∂⎩(2-40) 即12122,2+1,1(,)(,)()()2(,)(,)()()2j j k j k k j j k j k k u u x t u x t o t u u x t u x t o t ττττ+++∂⎧=+⋅+⎪⎪∂⎨∂⎪=-⋅+⎪∂⎩(2-41) 得到如下方程()()12122,12,12(,)(,)()=()2(,)(,)()=()j j k k j k j j k k j k u x t u x t u o t u x t u x t u o tττττ++++⎧⎡⎤-∂⎪⎢⎥+⎪⎢⎥∂⎢⎥⎪⎣⎦⎨⎡⎤⎪-∂⎢⎥⎪-+⎢⎥∂⎪⎢⎥⎣⎦⎩(2-42) 所以12k t t +=处的一阶偏导数可以用下式表示：()()112212,1,122(,)(,)2(,)(,)11()()()22(,)(,)()j j k j j k k k j k j k j k j k u x t u x t u x t u x t u u o t t u x t u x t o τττττ+++++⎡⎤--∂∂⎡⎤⎢⎥+=-+⎢⎥⎢⎥∂∂⎣⎦⎢⎥⎣⎦-=+(2-43)将k t t =,11j j x x x x +-==和代入式(2-6)可以得到式(2-9)； 同理1k t t +=,11j j x x x x +-==和代入式(2-6)可以得到2111112,122(,)2(,)(,)()()j k j k j k j k u x t u x t u x t u o h x h+++-++-+∂=+∂(2-44) 所以12k t t +=，j x 处的二阶偏导数用式(2-6)、(2-44)表示：(2-45)所以12k t t +=，j x 处的函数值可用下式表示： 11(,)(,)2j k j k f x t f x t +⎡⎤+⎣⎦(2-46) 原方程变为：()2212,1,,,12211()()()()()222k k j k j k j k j k j j u u u u f f o h t t x x α+++⎡⎤∂∂∂∂⎡⎤+-+=++⎢⎥⎢⎥∂∂∂∂⎣⎦⎣⎦(2-47) 将差分格式代入上式得：1111111122221(,)(,)(,)2(,)(,)(,)2(,)(,)21(,)(,)()2j k j k j k j k j k j k j k j k j k j k u x t u x t u x t u x t u x t u x t u x t u x t h h f x t f x t o h αττ++++-++-+--+-+⎡⎤-+⎢⎥⎣⎦⎡⎤=+++⎣⎦(2-48)为了方便写成：()11111111112222k k k k k k k k k k j j j j j j j j j j r u u u u u u u u f f ατ++++++-+-⎡⎤⎡⎤---++-+=+⎣⎦⎢⎥⎣⎦(2-49) 最后得到Crank-Nicholson 格式的差分格式为()()()()()1111111111=1+222k k k k k k k k j j j j j j j j r r r u u u r u u u f f αααατ+++++-+-⎡⎤+-+-+++⎢⎥⎣⎦(2-50) 2r hτ=其中，Crank-Nicholson 格式的差分格式的截断误差是22()o h τ+。