第8章-偏微分方程数值解.ppt

第第8章章   偏微分方程数值解偏微分方程数值解一、典型的偏微分方程介绍一、典型的偏微分方程介绍1. 椭圆型方程椭圆型方程: 在研究有热源稳定状态下的热传导，在研究有热源稳定状态下的热传导，有固定外力作用下薄膜的平衡问题时，都会遇到有固定外力作用下薄膜的平衡问题时，都会遇到Poisson方程方程 Laplace 方程方程Poisson方程方程 2021/3/1112. 抛物型方程抛物型方程 热传导方程热传导方程热传导方程热传导方程: :在研究热传导过程、气体扩散现在研究热传导过程、气体扩散现象、电磁场的传播等问题中以及在统计物理、概率象、电磁场的传播等问题中以及在统计物理、概率论和重子力学中，经常遇到抛物型方程这类方程论和重子力学中，经常遇到抛物型方程这类方程中最简单、最典型的是热传导方程中最简单、最典型的是热传导方程 其中其中a a是常数它表示长度为它表示长度为L L的细杆内，物体温度分布的规律的细杆内，物体温度分布的规律 2021/3/1123．．双曲型方程双曲型方程波动方程波动方程( (波的传播、物体的振动波的传播、物体的振动) )它表示长度为它表示长度为L L的弦振动的规律。

的弦振动的规律 2021/3/113二、定解问题二、定解问题决定方程唯一解所必须给定的决定方程唯一解所必须给定的初始条件初始条件初始条件初始条件和和边界条件边界条件边界条件边界条件叫做定解条件叫做定解条件. .定解条件由实际问题提出定解条件由实际问题提出. .解条件解条件 • 抛物型方程边界条件的提法应为物体在端点的温度抛物型方程边界条件的提法应为物体在端点的温度分布为已知，即边界条件分布为已知，即边界条件2021/3/114•双曲型方程初始条件表示弦在两端振动规律为已知：双曲型方程初始条件表示弦在两端振动规律为已知：2021/3/115nPoisson方程反应稳定状态的情况，与时间无关，所以不需要提初始条件边界条件的提法为：n其中 (x, y)为已知边界，s是区域D的边界 2021/3/116n本章主要针对几个典型的微分方程介绍常用的差分方法和有限元方法n这些方法基本思想是：q把一个连续问题离散化q通过各种手法化成有限形式的线性方程组q然后求其解 2021/3/117u计算机只能作有限次的加、减、乘、除运算，计算机只能作有限次的加、减、乘、除运算，计算机只能作有限次的加、减、乘、除运算，计算机只能作有限次的加、减、乘、除运算，它既不能求导数，更不能解偏微分方程它既不能求导数，更不能解偏微分方程它既不能求导数，更不能解偏微分方程它既不能求导数，更不能解偏微分方程u如果想在计算机上求得微分方程数值解，它的如果想在计算机上求得微分方程数值解，它的如果想在计算机上求得微分方程数值解，它的如果想在计算机上求得微分方程数值解，它的主要做法是主要做法是主要做法是主要做法是u把偏微分方程中所有的偏导数分别用差商代替u从而得到一个代数方程组——差分方程组u然后对差分方程求解，并以所求的解作为偏微分方程数值解8.1 差分法简介差分法简介2021/3/118 因此需要对区域进行剖分，用网格点来代替连续区因此需要对区域进行剖分，用网格点来代替连续区域，所以差分法亦称域，所以差分法亦称“网格法网格法网格法网格法”。

0xy2021/3/119 把整体分割成若干个单元来处理问题的把整体分割成若干个单元来处理问题的方法在数学上称为方法在数学上称为““离散化方法离散化方法” ” 在结点上采用离散化方法（数值微分、在结点上采用离散化方法（数值微分、数值积分、泰勒展开等）将微分方程的初边数值积分、泰勒展开等）将微分方程的初边值问题化成关于离散变量的相应问题，这个值问题化成关于离散变量的相应问题，这个相应问题的解就是方程在点相应问题的解就是方程在点xi上的数值解上的数值解f(x)，或在点（，或在点（xi , ti）上的数值解）上的数值解U( xi , ti) 一般来说，不同的离散化导致不同的方一般来说，不同的离散化导致不同的方法2021/3/1110例：取一边长为例：取一边长为1 1的正方形均匀薄板，上下侧面绝热，的正方形均匀薄板，上下侧面绝热，四周保持恒温，求板内各点的稳定温定分布四周保持恒温，求板内各点的稳定温定分布u=0u=0u=00xyLaplace 方程第一边值问题方程第一边值问题2021/3/1111记记u在这些点满足方程在这些点满足方程 2021/3/1112得到得到u (i, k)的近似的近似ui,k，所满足的线性代数方程组：，所满足的线性代数方程组： 其中其中 用迭代法来解方程组用迭代法来解方程组 2021/3/1113简单迭代法简单迭代法高斯高斯—赛德尔迭代法赛德尔迭代法2021/3/1114表8.1表8.200000k=400000000.35400.707k=300.1510.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/1115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=42021/3/1116000000.7070.4530.2580.151010.573 0.3860.182 00.7070.3810.2430.134000000表表8.4i=0i=1i=2i=3i=4k=0k=1k=2k=3k=42021/3/1117用差分法解偏微分方程需要考虑三个问题：用差分法解偏微分方程需要考虑三个问题： 1 1．选用网格，将微分方程离散化为差分方程。

．选用网格，将微分方程离散化为差分方程2 2．当网格步长．当网格步长h 0时差分方程的准确解是否时差分方程的准确解是否收敛于微分方程的解？收敛于微分方程的解？ 3 3．如何解相应的代数方程组？．如何解相应的代数方程组？ 2021/3/11188.2 8.2 8.2 8.2 椭圆型方程的差分解法椭圆型方程的差分解法椭圆型方程的差分解法椭圆型方程的差分解法 椭圆型方程最简单的典型问题就是椭圆型方程最简单的典型问题就是拉普拉斯方程拉普拉斯方程拉普拉斯方程拉普拉斯方程泊松方程泊松方程泊松方程泊松方程 2021/3/1119考虑泊松方程第一边值问题：考虑泊松方程第一边值问题：2021/3/1120( (一一) ) 矩形网格矩形网格          设设 为为xy平面上一有界区域平面上一有界区域，，为其边界，为其边界，是分段光滑曲线是分段光滑曲线0xy正则内点正则内点非正则内点非正则内点边界点边界点2021/3/1121（二）（二）五点差分格式五点差分格式        现在假设现在假设（（i，，k））为为正则内点正则内点沿着x，，y轴方向分别用轴方向分别用二阶中心差商二阶中心差商代替代替uxx，，uyy，，则得则得若以若以uh，，fh表示表示网函数网函数，，记记2021/3/1122则差分方程可简写成：则差分方程可简写成： 利用利用Taylor展式展式 2021/3/11232021/3/1124这四个式子两两相加便有：这四个式子两两相加便有： 2021/3/1125于是可得差分方程的截断误差于是可得差分方程的截断误差 2021/3/1126（三）（三）边值条件的处理边值条件的处理 以第一边值条件以第一边值条件 为例为例：非正则内点集合：非正则内点集合 h ：：边界点集合边界点集合 （（1））直接转移法直接转移法对对（（xi, yk））  ，用边界上距离这点最近的点的值，用边界上距离这点最近的点的值作为作为（（xi, yk））的值，即的值，即2021/3/1127（（2））线性插值法线性插值法641352h12021/3/1128则则u在这些点上的值有近似关系：在这些点上的值有近似关系： 2021/3/1129（（3））列不等距差分方程列不等距差分方程f1为为f在在1点的值点的值。

 2021/3/11308.3   8.3   抛物型方程的差分解法抛物型方程的差分解法抛物型方程的差分解法抛物型方程的差分解法   抛物型方程是指如下形式的方程：抛物型方程是指如下形式的方程： 很多实际的物理问题都可以用这类方程描述：很多实际的物理问题都可以用这类方程描述：热传导方程热传导方程2021/3/1131现以热传导方程为例，介绍抛物型方程的现以热传导方程为例，介绍抛物型方程的有限差分格式有限差分格式有限差分格式有限差分格式设热传导方程：设热传导方程：定解条件定解条件(1)(2)求求（（1））满足满足（（2））的解2021/3/1132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上的结点称为边界结点，其余所有属于上的结点称为边界结点，其余所有属于内的结点称为内部结点内的结点称为内部结点txoh (xj, tk)2021/3/11338.3.2．古典差分格式古典差分格式于平面区域于平面区域 上考虑传导方程：上考虑传导方程： a 为正常数为正常数 （（3）） (4)2021/3/1134于结点于结点（（j, k））处偏导数与差商之间有如下近似的关系：处偏导数与差商之间有如下近似的关系：利用上述表达式得到利用上述表达式得到 LU 在在 (j, k) 处的关系式：处的关系式： (5)2021/3/1135视为视为 u (xj, tk) 的近似值的近似值 令令，j = 1, 2,…, N – 1；；k = 0, 1, 2, …则有：则有：（（6））差分方程差分方程（（6））称为解热传导方程称为解热传导方程（（3））的古典显格式，的古典显格式，它所用到的结点如下图：它所用到的结点如下图： * *        *        * （（j, k））2021/3/1136将将（（6））写成便于计算的格式：写成便于计算的格式：（（7））称为称为网比网比网比网比，利用，利用（（7））及初边值条件及初边值条件（（4））在网格上的值在网格上的值（（8））即可算出即可算出k = 1, 2,…，，各层上的值各层上的值 。

截断误差阶为截断误差阶为 0(  + h2) 2021/3/1137为了提高截断误差的阶，可以利用中心差商：为了提高截断误差的阶，可以利用中心差商：j = 1, 2,……, N – 1；；k = 0, 1, 2,…… （（9））得到得到 Richardson  格式，其结点图为：格式，其结点图为： * * * * （j, k） *2021/3/1138截断误差阶为截断误差阶为o ( 2 + h2)，，较古典显格式高较古典显格式高将将（（9））式改写成适于计算的形式：式改写成适于计算的形式：j = 1, 2,……, N – 1；； k = 1, 2,…… r = a  /h2  称为网比，（称为网比，（10））式中出现了三层网格上的值，式中出现了三层网格上的值，（（10））才能逐层计算才能逐层计算故需要事先求得第故需要事先求得第k－－1层的值层的值 和第和第k层的值层的值,2021/3/1139如果利用向后差商如果利用向后差商 j = 1, 2,……, N – 1；；k = 0, 1, 2,…… (11)(12)j = 1, 2,……, N – 1；；k = 0, 1, 2,…… 古典隐格式古典隐格式古典隐格式古典隐格式，其结点图为：，其结点图为： （（j, k））  *       *       *           *截断误差为截断误差为o (  + h2)，与古典显格式相同。

与古典显格式相同 2021/3/11408.3.3 六点对称格式六点对称格式取该点的中心差商，从而取该点的中心差商，从而对于方程对于方程（（3））式式，在在                  点列方程点列方程， ，a 为正常数为正常数 （（3）） 2021/3/1141将以上各式代入将以上各式代入（（3））式得到差分方程：式得到差分方程： 2021/3/1142整理整理, ,得得 此即六点对称格式，也称为此即六点对称格式，也称为Crank-Nicolson格式，所用结点图为：格式，所用结点图为： *       *      *    k + 1*       *      *    k            j + 1   j j – 1（（13））2021/3/11438.3.4 稳定性稳定性（（1））当步长无限缩小时，差分方程的解是否逼近于微分方程当步长无限缩小时，差分方程的解是否逼近于微分方程（（2））计算过程中产生的误差在以后的计算中是无限增加，计算过程中产生的误差在以后的计算中是无限增加，还是可以控制？（稳定性）还是可以控制？（稳定性）的解？（收敛性）的解？（收敛性）稳定性问题是研究抛物型差分方程的一个中心课题！稳定性问题是研究抛物型差分方程的一个中心课题！稳定性问题是研究抛物型差分方程的一个中心课题！稳定性问题是研究抛物型差分方程的一个中心课题！ 2021/3/1144考察考察 Richardson 格式的稳定性格式的稳定性。

用用 表示计算表示计算 所产生的误差，如果右端所产生的误差，如果右端 无误差存在，无误差存在，则则 满足：满足：取取（（14））假设假设k - 1层之前无误差存在即层之前无误差存在即 ，而在第，而在第k层产生了层产生了误差误差 ，这一层其它点也无误差，而且在计算过程这一层其它点也无误差，而且在计算过程中不再产生新的误差，利用（中不再产生新的误差，利用（14）式算出误差）式算出误差  的传播如下表：的传播如下表： 2021/3/1145 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/1146 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 = ½ 时的古典显格式，误差方程为时的古典显格式，误差方程为：2021/3/1147差分格式关于初值稳定的实际含义是：如果其解差分格式关于初值稳定的实际含义是：如果其解在某一层存在误差，则由它引起的以后各层上的误差不在某一层存在误差，则由它引起的以后各层上的误差不超过原始误差的超过原始误差的M倍（倍（M为与为与  无关的常数）。

无关的常数）因此，在稳定的条件下，只要初始误差足够小，因此，在稳定的条件下，只要初始误差足够小，以后各层的误差也能足够小以后各层的误差也能足够小以上构造的几种差分格式中，以上构造的几种差分格式中，古典显格式：古典显格式：r ≤1/2时稳定时稳定古典隐格式：绝对稳定古典隐格式：绝对稳定Richardson格式：绝对不稳定格式：绝对不稳定六点对称格式：绝对稳定六点对称格式：绝对稳定稳定性概念稳定性概念稳定性概念稳定性概念：：2021/3/11488.4  双曲型方程的差分解法双曲型方程的差分解法 一阶线性双曲型方程最简单的形式为一阶线性双曲型方程最简单的形式为（（8.4.1)当给定初始条件当给定初始条件(8.4.2)以后，容易验证，双曲型方程以后，容易验证，双曲型方程(8.4.1)的解为：的解为：(8.4.3)2021/3/1149也就是说，在平面也就是说，在平面 x t上，沿着上，沿着（（k 是常数是常数）） 这样的直线，这样的直线，u 的值保持不变这种直线叫做的值保持不变这种直线叫做特征线特征线特征线特征线 0xta>00xta<0这是个单向的传播波，这是个单向的传播波，a>0时，波形时，波形 （（x））沿沿x轴方向传播轴方向传播，，为右传播波为右传播波，，a < 0时，为左传播波，在传播过程中，波形时，为左传播波，在传播过程中，波形均不发生变化。

均不发生变化2021/3/1150在物理上常见的双曲型偏微分方程最简单模型是波动方程在物理上常见的双曲型偏微分方程最简单模型是波动方程 其中其中 ，如果引进变量，如果引进变量则得到与波动方程等价的方程组则得到与波动方程等价的方程组(8.4.4)(8.4.6)(8.4.5)2021/3/11518.4.1 矩形网格矩形网格用两组平行直线族用两组平行直线族xj = jh, , tn = n  (j = 0,  1, …，，n = 0，，1,2 …)构成的矩形构成的矩形网覆盖了网覆盖了xt平面，网格点平面，网格点（（xj, tn））称为结点，简记为称为结点，简记为（（j, n）），，h、、 为常数，分别为常数，分别称为空间步长及时间步长，称为空间步长及时间步长，或称或称h为沿为沿x方向的步长，称方向的步长，称 为沿为沿t方向的步长，，方向的步长，，N为正整数在为正整数在t = 0上的结点称为边界结点，其余所有属于上的结点称为边界结点，其余所有属于内的结点称为内部结点内的结点称为内部结点txoh (xj, tn)2021/3/1152a））迎风格式迎风格式 ut (xj, tn)用向前差商代替用向前差商代替，，ux（（xj, tn））用向前或向后用向前或向后差商代替，得差商代替，得8.4.2 一阶双曲方程的差分法一阶双曲方程的差分法或或 2021/3/1153令令r=   / h，，得得截断误差均为截断误差均为 节点分布节点分布图：图： * *        *        * （（j, n））(8.4.7)(8.4.8)2021/3/1154稳定性稳定性的讨论：的讨论：令令代入代入（（8.4.7））得传播因子得传播因子2021/3/1155当当a>0时，恒有时，恒有 ，格式，格式（（8.4.7））不稳定不稳定 ；；当当a<0且且ar   1时时，， ，，格式格式（（8.4.7））稳定稳定。

 格式格式（（8.4.8））在在a <0时不稳定时不稳定，，在在a >0且且ar   1时稳定 将迎风格式写为统一形式：将迎风格式写为统一形式： 稳定性条件为：稳定性条件为：(8.4.9)2021/3/1156b)  Lax-Friedrichs格式格式该格式构造于该格式构造于1954年，用到年，用到 的技巧，截断误差为的技巧，截断误差为 ：： 节点分布节点分布图：图： * *        *        * （（j, n））(8.4.10)2021/3/1157传播因子传播因子 时稳定时稳定当当               时时，，                     ，，即格式即格式（（8.4.10））在在2021/3/1158c））Lax-Wendroff格式格式截断误差为截断误差为 节点分布图：节点分布图： * *        *        * （（j, n））传播因子传播因子当当 时有时有 ，即格式在，即格式在 条件下稳定。

条件下稳定 2021/3/1159d））古典隐式格式古典隐式格式ut用向后差商代替用向后差商代替，，ux用中心差商代替用中心差商代替得得截断误差为：截断误差为： 传播因子：传播因子：对任意的网格比，均有对任意的网格比，均有 ，故古典隐格式绝对稳定故古典隐格式绝对稳定2021/3/1160e））Grank-Nicholson格式格式在在（（ ））处展开，由处展开，由及中心差商以式而得到：及中心差商以式而得到：截断误差为截断误差为 ：：绝对稳定绝对稳定2021/3/11618.4.3  二阶双曲型方程的差分格式二阶双曲型方程的差分格式直接构造方程直接构造方程 的差分格式的差分格式 utt, uxx均用中心差商代替之，得均用中心差商代替之，得其中网格比其中网格比 : :截断误差截断误差 2021/3/1162b））隐格式隐格式利用关系利用关系可得三层隐式格式：可得三层隐式格式： 截断误差截断误差 : :绝对稳绝对稳定定2021/3/1163考试就要到了*******n祝同学们期末考试愉快！2021/3/1164。