第8章偏微分方程数值解

《第8章偏微分方程数值解》由会员分享，可在线阅读，更多相关《第8章偏微分方程数值解（64页珍藏版）》请在装配图网上搜索。

1、2021/3/111一、典型的偏微分方程介绍一、典型的偏微分方程介绍1.椭圆型方程椭圆型方程:在研究有热源稳定状态下的热传导，在研究有热源稳定状态下的热传导，有固定外力作用下薄膜的平衡问题时，都会遇到有固定外力作用下薄膜的平衡问题时，都会遇到Poisson方程方程 Dyxyxfyuxu),(),(222202222yuxuLaplace 方程方程Poisson方程方程 2021/3/1122.抛物型方程抛物型方程 :在研究热传导过程、气体扩散现在研究热传导过程、气体扩散现象、电磁场的传播等问题中以及在统计物理、概率象、电磁场的传播等问题中以及在统计物理、概率论和重子力学中，经常遇到抛物型方程。

2、这类方程论和重子力学中，经常遇到抛物型方程。这类方程中最简单、最典型的是热传导方程。中最简单、最典型的是热传导方程。Lxttxfxuatu0,0),(22其中其中a a是常数。是常数。它表示长度为它表示长度为L L的细杆内，物体温度分布的规律的细杆内，物体温度分布的规律 2021/3/1133双曲型方程双曲型方程Lxttxfxuatu0,0),(22222波动方程波动方程(波的传播、物体的振动波的传播、物体的振动)它表示长度为它表示长度为L L的弦振动的规律。的弦振动的规律。2021/3/114二、定解问题二、定解问题决定方程唯一解所必须给定的决定方程唯一解所必须给定的和和叫做定解条件叫做定解

3、条件.定解条件由实际问题提出定解条件由实际问题提出.解条件解条件 抛物型方程边界条件的提法应为物体在端点的温度抛物型方程边界条件的提法应为物体在端点的温度分布为已知，即边界条件分布为已知，即边界条件0)(),()(),0(tttLuttu2021/3/115 双曲型方程初始条件表示弦在两端振动规律为已知：双曲型方程初始条件表示弦在两端振动规律为已知：Lxxgxtuxfxu0)()0,()()0,(2021/3/116nPoisson方程反应稳定状态的情况，与时间无关，所以不需要提初始条件。边界条件的提法为：n其中(x,y)为已知边界，s是区域D的边界。),(),(yxyxus2021/3/11

4、7n本章主要针对几个典型的微分方程介绍常用的差分方法和有限元方法。n这些方法基本思想是：q把一个连续问题离散化q通过各种手法化成有限形式的线性方程组q然后求其解。2021/3/118u把偏微分方程中所有的偏导数分别用差商代替u从而得到一个代数方程组差分方程组u然后对差分方程求解，并以所求的解作为偏微分方程数值解8.1 差分法简介差分法简介2021/3/119 因此需要对区域进行剖分，用网格点来代替连续区因此需要对区域进行剖分，用网格点来代替连续区域，所以差分法亦称域，所以差分法亦称“”。0 xy2021/3/1110 把整体分割成若干个单元来处理问题的把整体分割成若干个单元来处理问题的方法在数

5、学上称为方法在数学上称为“离散化方法离散化方法”。在结点上采用离散化方法（数值微分、在结点上采用离散化方法（数值微分、数值积分、泰勒展开等）将微分方程的初边数值积分、泰勒展开等）将微分方程的初边值问题化成关于离散变量的相应问题，这个值问题化成关于离散变量的相应问题，这个相应问题的解就是方程在点相应问题的解就是方程在点xi上的数值解上的数值解f(x)，或在点（或在点（xi,ti）上的数值解）上的数值解U(xi,ti)。一般来说，不同的离散化导致不同的方一般来说，不同的离散化导致不同的方法。法。2021/3/1111例：取一边长为例：取一边长为1 1的正方形均匀薄板，上下侧面绝热，的正方形均匀薄板

6、，上下侧面绝热，四周保持恒温，求板内各点的稳定温定分布。四周保持恒温，求板内各点的稳定温定分布。u=0u=0yusinu=00 xy)10,10:(sin0011002222yxyuuuuyuxuuxyyxLaplace 方程第一边值问题方程第一边值问题2021/3/1112222),1(),(2),1(hkiukiukiuxuik222)1,(),(2)1,(hkiukiukiuyuik),(,kiyyxukhihyxkiki记记u在这些点满足方程在这些点满足方程 02222ikikyuxu2021/3/1113得到得到u(i,k)的近似的近似ui,k，所满足的线性代数方程组：，所满足的线性

7、代数方程组：0411,1,1,12ikkikikikiuuuuuh其中其中 )3,2,1,(0400ikuuuiik3707.0211707.04sin;4kkkkuk用迭代法来解方程组用迭代法来解方程组 1;1;1;141kikikikiikuuuuu2021/3/1114)(1;)(1,)(;1)(;1)1(41nkinkinkinkinikuuuuu简单迭代法简单迭代法)(1,)(,1)1(1,)1(,1)1(41nkinkinkinkinikuuuuu高斯高斯赛德尔迭代法赛德尔迭代法2021/3/1115表8.1表8.200000k=400000000.35400.707k=300.1

8、510.3540.4530.70700.2500.751k=200.250.4270.751000.35400.707k=100.1510.3540.4530.70700000k=000000i=0i=1i=2i=3i=4u(0)2021/3/1116000000.7070.4530.2580.151010.583 0.4270.182 00.7070.453 0.2580.151000000表表8.3i=0i=1i=2i=3i=4k=0k=1k=2k=3k=42021/3/1117000000.7070.4530.2580.151010.573 0.3860.182 00.7070.3810

9、.2430.134000000表表8.4i=0i=1i=2i=3i=4k=0k=1k=2k=3k=42021/3/1118用差分法解偏微分方程需要考虑三个问题：用差分法解偏微分方程需要考虑三个问题：1 1选用网格，将微分方程离散化为差分方程。选用网格，将微分方程离散化为差分方程。2 2当网格步长当网格步长h 0时差分方程的准确解是时差分方程的准确解是否收敛于微分方程的解？否收敛于微分方程的解？3 3如何解相应的代数方程组？如何解相应的代数方程组？2021/3/1119 椭圆型方程最简单的典型问题就是椭圆型方程最简单的典型问题就是02222yuxuu),(2222yxfyuxuu2021/3/1

10、120考虑泊松方程第一边值问题：考虑泊松方程第一边值问题：),(),(),(),(2222yxyxuyxyxfyuxuu2021/3/1121(一一)矩形网格矩形网格 设设 为为xy平面上一有界区域平面上一有界区域，为其边界，为其边界，是分段光滑曲线。是分段光滑曲线。0 xy正则内点正则内点非正则内点非正则内点边界点边界点2021/3/1122（二）（二）五点差分格式五点差分格式 现在假设现在假设（i，k）为为正则内点正则内点。沿着。沿着x，y轴方向分别用轴方向分别用二阶中心差商二阶中心差商代替代替uxx，uyy，则得则得ikkiikkikiikkiikhfhuuuhuuuu221,1,21,

11、1,122若以若以uh，fh表示表示网函数网函数，记记),(),(,),(kiikkihikkihyxffyxfuyxu2021/3/1123则差分方程可简写成：则差分方程可简写成：hhhfu 利用利用Taylor展式展式),1(),(),(!6!5!4!3!2),(66615551444133312221,1kikikixuhxuhxuhxuhxuhxuhuukiikikikikikikki2021/3/1124),(),(),1(!6!5!4!3!2),(6661555144413331222121,1kikikixuhxuhxuhxuhxuhxuhuukiikikikikikikki)1

12、,(),(),(!6!5!4!3!2),(6662555244423332222221,kikikiyuhyuhyuhyuhyuhyuhuukiikikikikikikki2021/3/1125),(),()1,(!6!5!4!3!2),(6662555244423332222221,kikikiyuhyuhyuhyuhyuhyuhuukiikikikikikikki这四个式子两两相加便有：这四个式子两两相加便有：2021/3/1126)(36012261664144212221,1,1hOxuhxuhxuhuuuikikikkiikki)(360122626642442222221,1,hO

13、yuhyuhyuhuuuikikikkiikki于是可得差分方程的截断误差于是可得差分方程的截断误差 )()(121),(),()(24,44224421hOhOyuhxuhyxuyxuuRkikihkiik2021/3/1127（三）（三）边值条件的处理边值条件的处理 以第一边值条件以第一边值条件 ),(yxu为例为例*h：非正则内点集合：非正则内点集合 h：边界点集合边界点集合（1）直接转移法直接转移法对对（xi,yk），用边界上距离这点最近的点的值，用边界上距离这点最近的点的值作为作为（xi,yk）的值，即的值，即*h),(kiikyxu2021/3/1128（2）线性插值法线性插值法6

14、41352h12021/3/1129则则u在这些点上的值有近似关系：在这些点上的值有近似关系：14112huuuu14121huhuu2021/3/1130（3）列不等距差分方程列不等距差分方程15132241112211fuuuhuuhuu)(211hf1为为f在在1点的值点的值。2021/3/1131抛物型方程是指如下形式的方程：抛物型方程是指如下形式的方程：),(222222zyxfzuyuxuatu 很多实际的物理问题都可以用这类方程描述：很多实际的物理问题都可以用这类方程描述：热传导方程热传导方程),(22txfxuatu 2021/3/1132现以热传导方程为例，介绍抛物型方程的现

15、以热传导方程为例，介绍抛物型方程的。设热传导方程：设热传导方程：定解条件定解条件),(22txfxuatu (1)0),1(),0()(0tutuxut(2)求求（1）满足满足（2）的解。的解。2021/3/11338.3.1 矩形网格矩形网格用两组平行直线族用两组平行直线族xj=jh,tk=k (j=0,1,，k=0，1,)构成的矩形构成的矩形网覆盖了网覆盖了xt平面，网格点平面，网格点（xj,tk）称为结点，简记为称为结点，简记为（j,k），h、为常数，分别为常数，分别称为称为，或称或称h为沿为沿x方向的步长，称方向的步长，称 为沿为沿t方向的步长，方向的步长，N为正整数。在为正整数。在t

16、=0上的结点称为边界结点，其余所有属于上的结点称为边界结点，其余所有属于内的结点称为内部结点。内的结点称为内部结点。Ttxc 0,txoh(xj,tk)2021/3/11348.3.2古典差分格式古典差分格式于平面区域于平面区域 上考虑传导方程：上考虑传导方程：TtxtxG 0,10),(22(,)uuUaf x ttx抖D=-=a 为正常数为正常数（3）Tttutuxxxu00),1(),0(10)()0,(4)2021/3/1135于结点于结点（j,k）处偏导数与差商之间有如下近似的关系：处偏导数与差商之间有如下近似的关系：)(),(),(1 otutxutxukjkjkj )(),(),

17、(2)(222211hoxuhtxutxutxukjkjkjkj 利用上述表达式得到利用上述表达式得到 LU 在在 (j,k)处的关系式：处的关系式：2111),(),(2)(),(),(htxutxutxuatxutxukjkjkjkjkj 2()kjUoht=D+(5)2021/3/1136),(kjkjtxff 视为视为 u(xj,tk)的近似值的近似值。kju令令，j=1,2,N 1；k=0,1,2,kjkjkjkjkjkjfhuuuauu 21112 则有：则有：（6）差分方程差分方程（6）称为解热传导方程称为解热传导方程（3）的古典显格式，的古典显格式，它所用到的结点如下图：它所用

18、到的结点如下图：*（j,k）2021/3/1137将将（6）写成便于计算的格式：写成便于计算的格式：kjkjkjkjkjkjfuuuruu 1112（7）2/har 称为称为，利用，利用（7）及初边值条件及初边值条件（4）在网格上的值在网格上的值 0)(00kNkjjjuuxu （8）即可算出即可算出k=1,2,，各层上的值各层上的值 。kju截断误差阶为截断误差阶为 0(+h2)。2021/3/1138为了提高截断误差的阶，可以利用中心差商：为了提高截断误差的阶，可以利用中心差商：)(2),(),(211 otutxutxukjkjkj j=1,2,N 1；k=0,1,2,kjkjkjkjk

19、jkjfhuuuauu 2111122（9）得到得到 Richardson 格式，其结点图为：格式，其结点图为：*（j,k）*2021/3/1139截断误差阶为截断误差阶为o(2+h2)，较古典显格式高较古典显格式高。将将（9）式改写成适于计算的形式：式改写成适于计算的形式：j=1,2,N 1；k=1,2,r=a /h2 称为网比，（称为网比，（10）式中出现了三层网格上的值，式中出现了三层网格上的值，kjkjkjkjkjkjfuuuuru 2221111 （10）才能逐层计算。才能逐层计算。故需要事先求得第故需要事先求得第k1层的值层的值 1 kjukju和第和第k层的值层的值,2021/3

20、/1140如果利用向后差商如果利用向后差商)(),(),(1 oxutxutxukjkjkj j=1,2,N 1；k=0,1,2,kjkjkjkjkjkjfhuuuauu 21112(11)kjkjkjkjkjkjfuuuuru 1112(12)j=1,2,N 1；k=0,1,2,，其结点图为：，其结点图为：（j,k）*截断误差为截断误差为o(+h2)，与古典显格式相同。，与古典显格式相同。2021/3/11418.3.3 六点对称格式六点对称格式取该点的中心差商，从而取该点的中心差商，从而对于方程对于方程（3）式式，在在 点列方程点列方程，21,kjtu 22xu )(212/1 ouuxu

21、kjkjkj )(22/12/12/12/12/122hohxuxuxukjkjkj )(122/112/12/12/11hohuuhuuhkjkjkjkj 22(,)uuUaf x ttx抖D=-=a 为正常数为正常数（3）2021/3/1142 kjkjkjuuu1112/1121 kjkjkjuuu 12/121 kjkjkjfff 12/121将以上各式代入将以上各式代入（3）式得到差分方程：式得到差分方程：2/12/112/12/11212 kjkjkjkjkjkjfuuuhauu 2/1111111122222 kjkjkjkjkjkjkjfuuuuuuha2021/3/1143

22、2/1111111121222 kjkjkjkjkjkjkjkjkjfuuuuuuhauu 111112)1(2 kjkjkjururur2/1112)1(2 kjkjkjkjfururur整理整理,得得 此即六点对称格式，也称为此即六点对称格式，也称为Crank-Nicolson格式，所用结点图为：格式，所用结点图为：*k+1*k j+1 j j 1（13）2021/3/11448.3.4 稳定性稳定性（1）当步长无限缩小时，差分方程的解是否逼近于微分方程当步长无限缩小时，差分方程的解是否逼近于微分方程（2）计算过程中产生的误差在以后的计算中是无限增加，计算过程中产生的误差在以后的计算中是无

23、限增加，还是可以控制？（稳定性）还是可以控制？（稳定性）的解？（收敛性）的解？（收敛性）2021/3/1145考察考察 Richardson 格式的稳定性格式的稳定性。用用 表示计算表示计算 所产生的误差，如果右端所产生的误差，如果右端 无误差存在，无误差存在，kje1 kjukjf则则 满足：满足：kje 111122 kjkjkjkjkjeeeere取取21 r11112 kjkjkjkjkjeeeee（14）假设假设k-1层之前无误差存在。即层之前无误差存在。即 ，而在第，而在第k层产生了层产生了01 kje误差误差。，这一层其它点也无误差，而且在计算过程这一层其它点也无误差，而且在计算

24、过程 kje0中不再产生新的误差，利用（中不再产生新的误差，利用（14）式算出误差）式算出误差 的传播如下表：的传播如下表：2021/3/1146 r=时时 Richardson 格式的误差传播格式的误差传播 j j0 4 j0 3 j0 2 j0 1 j0 j0+1 j0+2 j0+3 j0+4k -2 -4 7 4 -6 17 -24 17 -6 -8 31 -68 89 -68 31 -8 -10 49 -144 277 -388 277 -144 49 -10 71 -260 641 -109 1311 -109 641 -260 71 2021/3/1147 r 1/2 时古典显格式

25、的误差传播时古典显格式的误差传播 j j0 4 j0 3 j0 2 j0 1 j0 j0+1 j0+2 j0+3 j0+4k 0.5 0 0.5 0.25 0 0.5 0 0.25 0.125 0 0.375 0 0.375 0 0.125 0.0625 0 0.25 0 0.375 0 0.25 0 0.0625 如果选用如果选用 r=时的古典显格式，误差方程为时的古典显格式，误差方程为：kjkjkjeee11121 2021/3/1148差分格式关于初值稳定的实际含义是：如果其解差分格式关于初值稳定的实际含义是：如果其解在某一层存在误差，则由它引起的以后各层上的误差不在某一层存在误差，则由

26、它引起的以后各层上的误差不超过原始误差的超过原始误差的M倍（倍（M为与为与 无关的常数）。无关的常数）。因此，在稳定的条件下，只要初始误差足够小，因此，在稳定的条件下，只要初始误差足够小，以后各层的误差也能足够小。以后各层的误差也能足够小。以上构造的几种差分格式中，以上构造的几种差分格式中，古典显格式：古典显格式：r 1/2时稳定时稳定古典隐格式：绝对稳定古典隐格式：绝对稳定Richardson格式：绝对不稳定格式：绝对不稳定六点对称格式：绝对稳定。六点对称格式：绝对稳定。：2021/3/11498.4 双曲型方程的差分解法双曲型方程的差分解法 一阶线性双曲型方程最简单的形式为一阶线性双曲型方

27、程最简单的形式为0 xuatuLu（8.4.1)当给定初始条件当给定初始条件xxxu-)()0,(8.4.2)以后，容易验证，双曲型方程以后，容易验证，双曲型方程(8.4.1)的解为：的解为：xtatxtxu ,0 )(),(8.4.3)2021/3/1150也就是说，在平面也就是说，在平面 x t上，沿着上，沿着katx（k 是常数是常数）这样的直线，这样的直线，u 的值保持不变。这种直线叫做的值保持不变。这种直线叫做。0 xta00 xta0时，波形时，波形（x）沿沿x轴方向传播轴方向传播，为右传播波为右传播波，a 0时，恒有时，恒有 ，格式，格式（8.4.7）不稳定不稳定 ；1),(2k

28、G当当a0且且ar 1时时，格式格式（8.4.7）稳定稳定。1),(2kG格式格式（8.4.8）在在a 0且且ar 1时稳定。时稳定。将迎风格式写为统一形式：将迎风格式写为统一形式：)(2)(2111njnjnjnjnjnjuuaauuaaruu稳定性条件为：稳定性条件为：1ra(8.4.9)2021/3/1157b)Lax-Friedrichs格式格式该格式构造于该格式构造于1954年，用到年，用到)(2111njnjnjuuu的技巧，截断误差为的技巧，截断误差为 ：)(22hoho节点分布节点分布图：图：*（j,n）02)(2111111huuuuunjnjnjnjnj(8.4.10)20

29、21/3/1158传播因子传播因子 khiarkhkGsincos),(khrakhrakhkG22222222sin)1(1sincos),(时稳定时稳定。当当 时时，即格式即格式（8.4.10）在在1ra1),(kG1ra2021/3/1159c）Lax-Wendroff格式格式)2(2)(21122111njnjnjnjnjnjnjuuurauurauu截断误差为截断误差为 )(22ho节点分布图：节点分布图：*（j,n）传播因子传播因子khiarkhrakGsin2sin21),(2222sin)1(41),(422222khrarakG当当 时有时有 ，即格式在，即格式在 条件下稳定

30、。条件下稳定。1ra1G1ra2021/3/1160d）古典隐式格式古典隐式格式ut用向后差商代替用向后差商代替，ux用中心差商代替用中心差商代替得得02111huuauunjnjnjnj截断误差为：截断误差为：)(2ho传播因子：传播因子：khrakharikG222sin1sin1),(khrakG2222sin11),(对任意的网格比，均有对任意的网格比，均有 ，故古典隐格式绝对稳定。，故古典隐格式绝对稳定。1G2021/3/1161e）Grank-Nicholson格式格式02221111111huuhuuauunjnjnjnjnjnj在在（）处展开，由处展开，由21,njtx2/)(

31、121njnjnjuuu及中心差商以式而得到：及中心差商以式而得到：截断误差为截断误差为 ：)(22ho绝对稳定绝对稳定2021/3/11628.4.3 二阶双曲型方程的差分格式二阶双曲型方程的差分格式直接构造方程直接构造方程 的差分格式的差分格式 22222xuatuutt,uxx均用中心差商代替之，得均用中心差商代替之，得21121122huuuauuunjnjnjnjnjnj11122221)()1(2njnjnjnjnjuuuraurau其中网格比其中网格比 :h截断误差截断误差 )(22ho2021/3/1163b）隐格式隐格式利用关系利用关系21221222),(),(21),(xtxuxtxuxtxunjnjnj可得三层隐式格式：可得三层隐式格式：21111121111122112222huuuhuuuauuunjnjnjnjnjnjnjnjnj)()1(24)1(211112212211221221122njnjnjnjnjnjnjuuauauuauaua截断误差截断误差 :)(21ho绝对稳绝对稳定定2021/3/1164考试就要到了*n祝同学们期末考试愉快！