第8章偏微分方程数值解教育研究

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

2、中，经常遇到抛物型方程。这类方程中最简单、最典型的是热传导方程。中最简单、最典型的是热传导方程。 lxttxfxuatu0, 0),(22其中其中a a是常数。是常数。它表示长度为它表示长度为l l的细杆内，物体温度分布的规律的细杆内，物体温度分布的规律 2章节课堂3双曲型方程双曲型方程lxttxfxuatu0, 0),(22222波动方程波动方程( (波的传播、物体的振动波的传播、物体的振动) )它表示长度为它表示长度为l l的弦振动的规律。的弦振动的规律。 3章节课堂二、定解问题二、定解问题决定方程唯一解所必须给定的决定方程唯一解所必须给定的和和叫做定解条件叫做定解条件. .定解条件由实际

3、问题提出定解条件由实际问题提出. .解条件解条件 抛物型方程边界条件的提法应为物体在端点的温度抛物型方程边界条件的提法应为物体在端点的温度分布为已知，即边界条件分布为已知，即边界条件0)(),()(), 0(tttluttu4章节课堂 双曲型方程初始条件表示弦在两端振动规律为已知：双曲型方程初始条件表示弦在两端振动规律为已知：lxxgxtuxfxu0)()0 ,()()0 ,(5章节课堂npoisson方程反应稳定状态的情况，与时间无关，所以不需要提初始条件。边界条件的提法为：n其中 (x, y)为已知边界，s是区域d的边界。 ),(),(yxyxus6章节课堂n本章主要针对几个典型的微分方程

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

5、用离散化方法（数值微分、在结点上采用离散化方法（数值微分、数值积分、泰勒展开等）将微分方程的初边数值积分、泰勒展开等）将微分方程的初边值问题化成关于离散变量的相应问题，这个值问题化成关于离散变量的相应问题，这个相应问题的解就是方程在点相应问题的解就是方程在点xi上的数值解上的数值解f(x)，或在点（或在点（xi , ti）上的数值解）上的数值解u( xi , ti)。 一般来说，不同的离散化导致不同的方一般来说，不同的离散化导致不同的方法。法。10章节课堂例：取一边长为例：取一边长为1 1的正方形均匀薄板，上下侧面绝热，的正方形均匀薄板，上下侧面绝热，四周保持恒温，求板内各点的稳定温定分布。四

6、周保持恒温，求板内各点的稳定温定分布。u=0u=0yusinu=00 xy) 10, 10:(sin0011002222yxyuuuuyuxuuxyyxlaplace 方程第一边值问题方程第一边值问题11章节课堂222), 1(),(2), 1(hkiukiukiuxuik222) 1,(),(2) 1,(hkiukiukiuyuik),(,kiyyxukhihyxkiki记记u在这些点满足方程在这些点满足方程 02222ikikyuxu12章节课堂得到得到u (i, k)的近似的近似ui,k，所满足的线性代数方程组：，所满足的线性代数方程组： 0411,1, 1, 12ikkikikikiu

7、uuuuh其中其中 )3 , 2 , 1,(0400ikuuuiik3707. 0211707. 04sin; 4kkkkuk用迭代法来解方程组用迭代法来解方程组 1;1; 1; 141kikikikiikuuuuu13章节课堂)(1;)(1,)(; 1)(; 1)1(41nkinkinkinkinikuuuuu简单迭代法简单迭代法)(1,)(, 1)1(1,)1(, 1)1(41nkinkinkinkinikuuuuu高斯高斯赛德尔迭代法赛德尔迭代法14章节课堂表8.1表8.200000k=400000000.35400.707k=300.1510.3540.4530.70700.2500.

8、751k=200.250.4270.751000.35400.707k=100.1510.3540.4530.70700000k=000000i=0i=1i=2i=3i=4u(0)15章节课堂000000.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=416章节课堂000000.7070.4530.2580.151010.573 0.3860.182 00.7070.3810.2430.134000000表表8.4i=0i=1i=2i=3i=4

9、k=0k=1k=2k=3k=417章节课堂用差分法解偏微分方程需要考虑三个问题：用差分法解偏微分方程需要考虑三个问题： 1 1选用网格，将微分方程离散化为差分方程。选用网格，将微分方程离散化为差分方程。2 2当网格步长当网格步长h 0时差分方程的准确解是时差分方程的准确解是否收敛于微分方程的解？否收敛于微分方程的解？ 3 3如何解相应的代数方程组？如何解相应的代数方程组？ 18章节课堂 椭圆型方程最简单的典型问题就是椭圆型方程最简单的典型问题就是02222yuxuu),(2222yxfyuxuu19章节课堂考虑泊松方程第一边值问题：考虑泊松方程第一边值问题： ),(),( ),(),(2222

10、yxyxuyxyxfyuxuu20章节课堂( (一一) ) 矩形网格矩形网格 设设 为为xy平面上一有界区域平面上一有界区域，为其边界，为其边界，是分段光滑曲线。是分段光滑曲线。0 xy正则内点正则内点非正则内点非正则内点边界点边界点21章节课堂（二）（二）五点差分格式五点差分格式 现在假设现在假设（i，k）为为正则内点正则内点。沿着。沿着x，y轴方向分别用轴方向分别用二阶中心差商二阶中心差商代替代替uxx，uyy，则得则得ikkiikkikiikkiikhfhuuuhuuuu221,1,21, 1, 122若以若以uh，fh表示表示网函数网函数，记记),(),(,),(kiikkihikki

11、hyxffyxfuyxu22章节课堂则差分方程可简写成：则差分方程可简写成： hhhfu 利用利用taylor展式展式 ), 1(),(),(! 6! 5! 4! 3! 2),(66615551444133312221, 1kikikixuhxuhxuhxuhxuhxuhuukiikikikikikikki23章节课堂),(),(), 1(! 6! 5! 4! 3! 2),(6661555144413331222121, 1kikikixuhxuhxuhxuhxuhxuhuukiikikikikikikki) 1,(),(),(! 6! 5! 4! 3! 2),(666255524442333

12、2222221,kikikiyuhyuhyuhyuhyuhyuhuukiikikikikikikki24章节课堂),(),() 1,(! 6! 5! 4! 3! 2),(6662555244423332222221,kikikiyuhyuhyuhyuhyuhyuhuukiikikikikikikki这四个式子两两相加便有：这四个式子两两相加便有： 25章节课堂)(36012261664144212221, 1, 1hoxuhxuhxuhuuuikikikkiikki)(360122626642442222221,1,hoyuhyuhyuhuuuikikikkiikki于是可得差分方程的截断误差

13、于是可得差分方程的截断误差 )()(121),(),()(24,44224421hohoyuhxuhyxuyxuurkikihkiik26章节课堂（三）（三）边值条件的处理边值条件的处理 以第一边值条件以第一边值条件 ),(yxu为例为例*h：非正则内点集合：非正则内点集合 h ：边界点集合边界点集合 （1）直接转移法直接转移法对对（xi, yk） ，用边界上距离这点最近的点的值，用边界上距离这点最近的点的值作为作为（xi, yk）的值，即的值，即*h),(kiikyxu27章节课堂（2）线性插值法线性插值法641352h128章节课堂则则u在这些点上的值有近似关系：在这些点上的值有近似关系：

14、 14112huuuu14121huhuu29章节课堂（3）列不等距差分方程列不等距差分方程15132241112211fuuuhuuhuu)(211hf1为为f在在1点的值点的值。 30章节课堂抛物型方程是指如下形式的方程：抛物型方程是指如下形式的方程： ),(222222zyxfzuyuxuatu 很多实际的物理问题都可以用这类方程描述：很多实际的物理问题都可以用这类方程描述：热传导方程热传导方程),(22txfxuatu 31章节课堂现以热传导方程为例，介绍抛物型方程的现以热传导方程为例，介绍抛物型方程的。设热传导方程：设热传导方程：定解条件定解条件),(22txfxuatu (1)0)

15、, 1 (), 0()(0tutuxut(2)求求（1）满足满足（2）的解。的解。32章节课堂8.3.1 矩形网格矩形网格用两组平行直线族用两组平行直线族xj = jh, , tk = k (j = 0, 1, ，k = 0， 1, )构成的矩形构成的矩形网覆盖了网覆盖了xt平面，网格点平面，网格点（xj, tk）称为结点，简记为称为结点，简记为（j, k），h、 为常数，分别为常数，分别称为称为，或称或称h为沿为沿x方向的步长，称方向的步长，称 为沿为沿t方向的步长，方向的步长，n为正整数。在为正整数。在t = 0上的结点称为边界结点，其余所有属于上的结点称为边界结点，其余所有属于内的结点称

16、为内部结点。内的结点称为内部结点。 ttxc 0,txoh (xj, tk)33章节课堂8.3.2古典差分格式古典差分格式于平面区域于平面区域 上考虑传导方程：上考虑传导方程： ttxtxg 0, 10),(22( , )uuuaf x ttx抖d=-=a 为正常数为正常数 （3） tttutuxxxu00), 1 (), 0(10)() 0 ,( (4)34章节课堂于结点于结点（j, k）处偏导数与差商之间有如下近似的关系：处偏导数与差商之间有如下近似的关系：)(),(),(1 otutxutxukjkjkj )(),(),(2)(222211hoxuhtxutxutxukjkjkjkj 利

17、用上述表达式得到利用上述表达式得到 lu 在在 (j, k) 处的关系式：处的关系式： 2111),(),(2)(),(),(htxutxutxuatxutxukjkjkjkjkj 2()kjuoht=d+(5)35章节课堂),(kjkjtxff 视为视为 u (xj, tk) 的近似值的近似值。 kju令令，j = 1, 2, n 1；k = 0, 1, 2, kjkjkjkjkjkjfhuuuauu 21112 则有：则有：（6）差分方程差分方程（6）称为解热传导方程称为解热传导方程（3）的古典显格式，的古典显格式，它所用到的结点如下图：它所用到的结点如下图： * * * * （j, k）

18、36章节课堂将将（6）写成便于计算的格式：写成便于计算的格式： kjkjkjkjkjkjfuuuruu 1112（7）2/har 称为称为，利用，利用（7）及初边值条件及初边值条件（4）在网格上的值在网格上的值 0)(00knkjjjuuxu （8）即可算出即可算出k = 1, 2,，各层上的值各层上的值 。kju截断误差阶为截断误差阶为 0( + h2)。 37章节课堂为了提高截断误差的阶，可以利用中心差商：为了提高截断误差的阶，可以利用中心差商：)(2),(),(211 otutxutxukjkjkj j = 1, 2, n 1；k = 0, 1, 2, kjkjkjkjkjkjfhuuu

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

20、k层的值层的值,39章节课堂如果利用向后差商如果利用向后差商 )(),(),(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)，与古典显格式相同。，与古典显格式相同。 40章节课堂8.3.3 六点对称格式六点对称格式取该点的中心差商，从而取该点的中心差商，从而对于方程对于方程（3）式式，

21、在在 点列方程点列方程， ， 21, kjtu 22xu )(212/1 ouuxukjkjkj )(22/12/12/12/12/122hohxuxuxukjkjkj )(122/112/12/12/11hohuuhuuhkjkjkjkj 22( , )uuuaf x ttx抖d=-=a 为正常数为正常数 （3） 41章节课堂 kjkjkjuuu1112/1121 kjkjkjuuu 12/121 kjkjkjfff 12/121将以上各式代入将以上各式代入（3）式得到差分方程：式得到差分方程： 2/12/112/12/11212 kjkjkjkjkjkjfuuuhauu 2/1111111

22、122222 kjkjkjkjkjkjkjfuuuuuuha42章节课堂 2/1111111121222 kjkjkjkjkjkjkjkjkjfuuuuuuhauu 111112)1(2 kjkjkjururur2/1112)1(2 kjkjkjkjfururur整理整理, ,得得 此即六点对称格式，也称为此即六点对称格式，也称为crank-nicolson格式，所用结点图为：格式，所用结点图为： * * * k + 1* * * k j + 1 j j 1（13）43章节课堂8.3.4 稳定性稳定性（1）当步长无限缩小时，差分方程的解是否逼近于微分方程当步长无限缩小时，差分方程的解是否逼近于

23、微分方程（2）计算过程中产生的误差在以后的计算中是无限增加，计算过程中产生的误差在以后的计算中是无限增加，还是可以控制？（稳定性）还是可以控制？（稳定性）的解？（收敛性）的解？（收敛性） 44章节课堂考察考察 richardson 格式的稳定性格式的稳定性。 用用 表示计算表示计算 所产生的误差，如果右端所产生的误差，如果右端 无误差存在，无误差存在，kje1 kjukjf则则 满足：满足：kje 111122 kjkjkjkjkjeeeere取取21 r11112 kjkjkjkjkjeeeee（14）假设假设k - 1层之前无误差存在。即层之前无误差存在。即 ，而在第，而在第k层产生了层产

24、生了01 kje误差误差。 ，这一层其它点也无误差，而且在计算过程这一层其它点也无误差，而且在计算过程 kje0中不再产生新的误差，利用（中不再产生新的误差，利用（14）式算出误差）式算出误差 的传播如下表：的传播如下表： 45章节课堂 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

25、-109 1311 -109 641 -260 71 46章节课堂 r 1/2 时古典显格式的误差传播时古典显格式的误差传播 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 47章节课堂差分格式关于初值稳定的实际含义是：如果其解差分格式关于初

26、值稳定的实际含义是：如果其解在某一层存在误差，则由它引起的以后各层上的误差不在某一层存在误差，则由它引起的以后各层上的误差不超过原始误差的超过原始误差的m倍（倍（m为与为与 无关的常数）。无关的常数）。因此，在稳定的条件下，只要初始误差足够小，因此，在稳定的条件下，只要初始误差足够小，以后各层的误差也能足够小。以后各层的误差也能足够小。以上构造的几种差分格式中，以上构造的几种差分格式中，古典显格式：古典显格式：r 1/2时稳定时稳定古典隐格式：绝对稳定古典隐格式：绝对稳定richardson格式：绝对不稳定格式：绝对不稳定六点对称格式：绝对稳定。六点对称格式：绝对稳定。：48章节课堂8.4 双

27、曲型方程的差分解法双曲型方程的差分解法 一阶线性双曲型方程最简单的形式为一阶线性双曲型方程最简单的形式为0 xuatulu（8.4.1)当给定初始条件当给定初始条件xxxu- )()0 ,(8.4.2)以后，容易验证，双曲型方程以后，容易验证，双曲型方程(8.4.1)的解为：的解为：xtatxtxu , 0 )(),(8.4.3)49章节课堂也就是说，在平面也就是说，在平面 x t上，沿着上，沿着katx（k 是常数是常数） 这样的直线，这样的直线，u 的值保持不变。这种直线叫做的值保持不变。这种直线叫做。 0 xta00 xta0时，波形时，波形 （x）沿沿x轴方向传播轴方向传播，为右传播波

28、为右传播波，a 0时，恒有时，恒有 ，格式，格式（8.4.7）不稳定不稳定 ；1),(2kg当当a0且且ar 1时时， ，格式格式（8.4.7）稳定稳定。 1),(2kg格式格式（8.4.8）在在a 0且且ar 1时稳定。时稳定。 将迎风格式写为统一形式：将迎风格式写为统一形式： )(2)(2111njnjnjnjnjnjuuaauuaaruu稳定性条件为：稳定性条件为：1ra(8.4.9)56章节课堂b) lax-friedrichs格式格式该格式构造于该格式构造于1954年，用到年，用到 )(2111njnjnjuuu的技巧，截断误差为的技巧，截断误差为 ： )(22hoho节点分布节点分布图：图： * * * * （j, n）02)(2111111huuuuunjnjnjnjnj(8.4.10)57章节课堂传播因子传播因子 khiarkhkgsincos),(khrakhrakhkg22222222sin)1 (1sincos),(时稳定时稳定。当当 时时， ，即格式即格式（8.4.10）在在1ra1),(kg1ra58章节课堂c）lax-wendroff格式格式)2(