有限差分法求解偏微分方程MATLAB

任课教师评语：签名：年月日南京理工大学课程考核论文课程名称：高等数值分析论文题目：有限差分法求解偏微分方程姓名：罗晨学号：115104000545成绩：有限差分法求解偏微分方程一、主要内容1.有限差分法求解偏微分方程，偏微分方程如一般形式的一维抛物线型方程：22(,)()uufxttx其中为常数具体求解的偏微分方程如下：22001(,0)sin()(0,)(1,)00uuxtxuxxututt2.推导五种差分格式、截断误差并分析其稳定性；3.编写MATLAB程序实现五种差分格式对偏微分方程的求解及误差分析；4.结论及完成本次实验报告的感想。二、推导几种差分格式的过程：有限差分法（finite-differencemethods）是一种数值方法通过有限个微分方程近似求导从而寻求微分方程的近似解。有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。推导差分方程的过程中需要用到的泰勒展开公式如下：()2100000000()()()()()()()......()(())1!2!!nnnfxfxfxfxfxxxxxxxoxxn(2-1)求解区域的网格划分步长参数如下：11kkkkttxxh(2-2)2.1古典显格式2.1.1古典显格式的推导由泰勒展开公式将(,)uxt对时间展开得2,(,)(,)()()(())iikikkkuuxtuxtttottt(2-3)当1ktt时有21,112,(,)(,)()()(())(,)()()ikikikkkkkikikuuxtuxtttotttuuxtot(2-4)得到对时间的一阶偏导数1,(,)(,)()=()ikikikuxtuxtuot(2-5)由泰勒展开公式将(,)uxt对位置展开得223,,21(,)(,)()()()()(())2!kikikiikiiuuuxtuxtxxxxoxxxx(2-6)当11iixxxx和时，代入式(2-6)得2231,1,1122231,1,1121(,)(,)()()()()(())2!1(,)(,)()()()()(())2!ikikikiiikiiiiikikikiiikiiiiuuuxtuxtxxxxoxxxxuuuxtuxtxxxxoxxxx(2-7)因为1kkxxh，代入上式得2231,,22231,,21(,)(,)()()()2!1(,)(,)()()()2!ikikikikikikikikuuuxtuxthhohxxuuuxtuxthhohxx(2-8)得到对位置的二阶偏导数2211,22(,)2(,)(,)()()ikikikikuxtuxtuxtuohxh(2-9)将式(2-5)、(2-9)代入一般形式的抛物线型偏微分方程得21112(,)(,)(,)2(,)(,)(,)()ikikikikikikuxtuxtuxtuxtuxtfxtohh(2-10)为了方便我们可以将式(2-10)写成11122kkkkkkiiiiiiuuuuufh(2-11)11122kkkkkkiiiiiiuuuuufh(2-12)最后得到古典显格式的差分格式为111(12)kkkkkiiiiiurauruuf(2-13)2rh其中，古典显格式的差分格式的截断误差是2()oh。2.1.2古典显格式稳定性分析古典显格式(2-13)写成矩阵形式为112kkkhhhuraIraCuf(2-14)12212,(,,......,,)kkkkkhNNruuuuuh其中。(1)(1)01010100101010NNC上面的C矩阵的特征值是：2cos()1,2,......,1CjhjN12HraIraC212=122cos()121cos()14sin1,2,......,12HjCrarararajhrajhjhrajN(2-15)使()1H，即2114sin12jhra102ra结论：当102ra时，所以古典显格式是稳定的。2.2古典隐格式2.2.1古典隐格式的推导将1ktt代入式(2-3)得21,11(,)(,)()()(())jkjkjkkkkkuuxtuxtttottt(2-16)21,(,)(,)()()jkjkjkuuxtuxtot(2-17)得到对时间的一阶偏导数1,(,)(,)()=()jkjkjkuxtuxtuot(2-18)将式(2-9)、(2-18)原方程得到11122(,)(,)(,)2(,)(,)(,)()jkjkjkjkjkjkuxtuxtuxtuxtuxtfxtohh(2-19)为了方便把(2-19)写成11122kkkkkjjjjjkjuuuuufh(2-20)11122kkkkkkjjjjjjuuuuufh(2-21)最后得到古典隐格式的差分格式为111(12)kkkkkjjjjjrauruuuf(2-22)2rh其中，古典隐格式的差分格式的截断误差是2()oh。2.2.2古典隐格式稳定性分析将古典隐格式(2-22)写成矩阵形式如下1212()kkkhhhraIraCuufrh(2-23)误差传播方程112kkhhraIraCvv(2-24)12,AraIraCBI所以误差方程的系数矩阵为1112HAraIraC11,2,......,1122cosHjjNrarajh使()1H，显然21122cos()112(1cos())114sin2Hjrarajhrajhjhra1Hj恒成立。结论：对于0r，即任意网格比下，古典隐格式是绝对稳定的。2.3Richardson格式2.3.1Richardson格式的推导将11kktttt和，代入式(2-3)得21,1121,11(,)(,)()()(())(,)(,)()()(())ikikikkkkkikikikkkkkuuxtuxtttotttuuxtuxtttottt(2-25)即21,21,(,)(,)()()(,)(,)()()ikikikikikikuuxtuxtotuuxtuxtot(2-26)由此得到可得211,(,)(,)()()2ikikikuxtuxtuot(2-27)将式(2-9)、(2-27)代入原方程得到下式2211112(,)(,)(,)2(,)(,)(,)()2ikikikikikikuxtuxtuxtuxtuxtfxtohh(2-28)为了方便可以把式(2-28)写成1111222kkkkkkiiiiiiuuuuufh(2-29)即111122kkkkkkiiiiiiuuuuufh(2-30)最后得到Richardson显格式的差分格式为1111222kkkkkkiiiiiiuruuuuf(2-31)2rh其中，古典显格式的差分格式的截断误差是22()oh。2.3.2Richardson稳定性分析将Richardson显格式(2-31)写成如下矩阵形式11222kkkkhhhhurCIuuf(2-32)误差传播方程矩阵形式1122kkkhhhkkhhvrCIvvvv(2-33)再将上面的方程组写成矩阵形式112(2)0kkkhkraCIIvwwIv(2-34)系数矩阵的特征值是4cos()4110jrajhra228sin102jhra(2-35)解得特征值为22241,28sin64sin4222jhjhrara(2-36)222412,=4sin+1+16sin122jhjhMaxrara(恒成立)(2-37)结论：上式对任意的网比都恒成立，即Richardson格式是绝对不稳定的。4.Crank-Nicholson格式3.4.1Crank-Nicholson格式的推导将12ktt代入式(2-9)得1112221112222,2+1,1+1+1(,)(,)()()(())(,)(,)()()(())jjkjkkkkkkjjkjkkkkkkuuxtuxtttotttuuxtuxtttottt(2-40)即12122,2+1,1(,)(,)()()2(,)(,)()()2jjkjkkjjkjkkuuxtuxtotuuxtuxtot(2-41)得到如下方程12122,12,12(,)(,)()=()2(,)(,)()=()jjkkjkjjkkjkuxtuxtuotuxtuxtuot(2-42)所以12ktt处的一阶偏导数可以用下式表示：112212,1,122(,)(,)2(,)(,)11()()()22(,)(,)()jjkjjkkkjkjkjkjkuxtuxtuxtuxtuuottuxtuxto(2-43)将ktt,11jjxxxx和代入式(2-6)可以得到式(2-9)；同理1ktt,11jjxxxx和代入式(2-6)可以得到2111112,122(,)2(,)(,)()()jkjkjkjkuxtuxtuxtuohxh(2-44)所以12ktt，jx处的二阶偏导数用式(2-6)、(2-44)表示：22,,12211111112221()()2(,)2(,)(,)(,)2(,)(,)1()2jkjkjkjkjkjkjkjkuuxxuxtuxtuxtuxtuxtuxtohhh(2-45)所以12ktt，jx处的函数值可用下式表示：11(,)(,)2jkjkfxtfxt(2-46)原方程变为：2212,1,,,12211()()()()()222kkjkjkjkjkjjuuuuffohttxx(2-47)将差分格式代入上式得：1111111122221(,)(,)(,)2(,)(,)(,)2(,)(,)21(,)(,)()2jkjkjkjkjkj