第九章-偏微分方程差分方法汇总
- 格式:doc
- 大小:386.00 KB
- 文档页数:14
偏微分方程数值方法偏微分方程(Partial Differential Equation, PDE)是数学中的一种重要的方程类型,它描述了一个函数的多个变量的变化关系。
解决偏微分方程的数值方法在科学和工程领域有着广泛的应用。
本文将介绍几种常见的偏微分方程数值方法,并对其进行详细阐述。
1. 差分法(Finite Difference Method):差分法是最早也是最直接的一种数值方法,它基于连续函数在一些点的导数可以用它的前向、后向或中心的差商来近似的思想。
偏微分方程的差分格式包括向前差分法、向后差分法和中心差分法等。
对于二维的偏微分方程,可以采用网格化的方式将空间离散化,然后利用差分法进行近似求解。
2. 有限元法(Finite Element Method):有限元法是一种基于原始形式或变分形式对偏微分方程进行离散化的方法。
在有限元法中,将求解域分割成许多小的、简单的几何单元,然后在每个单元上构建近似解函数和试验函数。
通过构建弱形式并应用基本的变分原理,可以得到离散化的方程组,并通过求解这个方程组来得到数值解。
3. 有限差分法(Finite Difference Method):有限差分法是一种将连续的偏微分方程离散化成差分方程的方法。
它与差分法的主要区别在于有限差分法不需要对求解域进行网格化,而是直接在连续的求解域上进行离散化。
将偏微分方程中的导数通过差商来近似,然后通过求解离散化的差分方程来得到数值解。
4. 有限体积法(Finite Volume Method):有限体积法是一种将偏微分方程离散化为离散体积元的方法。
在有限体积法中,将求解域划分成离散的控制体积,然后通过对控制体积的积分运算,将偏微分方程转化为离散的代数方程组。
然后通过求解得到的代数方程组,可以得到数值解。
以上介绍的只是几种常见的偏微分方程数值方法,实际上还有很多其他的方法,如边界元法(Boundary Element Method)、谱方法(Spectral Method)、逆问题方法(Inverse Problem Method)等。
偏微分方程的求解方法偏微分方程(Partial Differential Equation,简称PDE)是一类重要的数学问题,其应用范围遍及自然科学、工程技术以及金融等领域。
如何求解偏微分方程是一个具有挑战性的问题,通常需要采用多种方法结合起来进行求解。
本文将简要介绍几种常见的偏微分方程求解方法。
1. 分离变量法分离变量法是一种简单而重要的偏微分方程求解方法。
该方法基于以下假设:偏微分方程的一个解可以写成一系列单一变量的函数乘积的形式。
具体地说,对于一个偏微分方程u(x, y) = 0(其中x, y为自变量),假设其解可以表示为u(x, y) = X(x)Y(y),其中X(x)和Y(y)分别是关于x和y的单一变量函数。
将u(x, y)代入原方程,得到X(x)Y(y) = 0。
由于0的任何一侧都是0,因此可得到两个单一变量方程:X(x) = 0和Y(y) = 0。
这两个方程的部分解(即使其中一个变量为常数时的解)可以结合在一起,形成原偏微分方程的一般解。
2. 特征线法特征线法是另一种重要的偏微分方程求解方法。
该方法的基本思想是将原方程转化为常微分方程,进而求解。
具体地说,对于一个二阶线性偏微分方程:a(x, y)u_xx + 2b(x, y)u_xy + c(x, y)u_yy + d(x, y)u_x + e(x, y)u_y + f(x, y)u = g(x, y),通过变量的代换,可以将该方程化为一个与一次微分方程组相关的形式。
进一步地,可以选择沿着特定的方向(例如x或y方向)进行参数化,从而得到关于变量的一阶微分方程。
该微分方程的解通常可以通过传统的常微分方程求解技巧来获得。
3. 数值方法数值方法是目前应用最广泛的偏微分方程求解方法之一。
由于大多数偏微分方程的解析解很难获得,因此数值方法成为了一种有效的、可行的替代方法。
常见的数值方法包括有限差分法、有限元法和边界元法等。
这些方法通过将偏微分方程离散化为一个有限维的计算问题,然后使用数值方法求解这个问题的解。
第9章 偏微分方程的差分方法含有偏导数的微分方程称为偏微分方程。
由于变量的增多和区域的复杂性,求偏微分方程的精确解一般是不可能的,经常采用数值方法求方程的近似解。
偏微分方程的数值方法种类较多,最常用的方法是差分方法。
差分方法具有格式简单,程序易于实现,计算量小等优点,特别适合于规则区域上偏微分方程的近似求解。
本章将以一些典型的偏微分方程为例,介绍差分方法的基本原理和具体实现方法。
9.1椭圆型方程边值问题的差分方法9.1.1 差分方程的建立最典型的椭圆型方程是Poisson (泊松)方程G y x y x f yux u u ∈=∂∂+∂∂-≡∆-),(),,()(2222 (9.1)G 是x ,y 平面上的有界区域,其边界Γ为分段光滑的闭曲线。
当f (x ,y )≡0时,方程(9.1)称为Laplace(拉普拉斯)方程。
椭圆型方程的定解条件主要有如下三种边界条件第一边值条件 ),(y x u α=Γ (9.2) 第二边值条件),(y x nuβ=∂∂Γ (9.3) 第三边值条件 ),()(y x ku nuγ=+∂∂Γ (9.4) 这里,n 表示Γ上单位外法向,α(x,y ),β(x,y ),γ(x,y )和k (x,y )都是已知的函数,k (x,y )≥0。
满足方程(9.1)和上述三种边值条件之一的光滑函数u (x ,y )称为椭圆型方程边值问题的解。
用差分方法求解偏微分方程,就是要求出精确解u (x ,y )在区域G 的一些离散节点(x i ,y i )上的近似值u i ,j ≈(x i ,y i )。
差分方法的基本思想是,对求解区域G 做网格剖分,将偏微分方程在网格节点上离散化,导出精确解在网格节点上近似值所满足的差分方程,最终通过求解差分方程,通常为一个线性方程组,得到精确解在离散节点上的近似值。
设G ={0<x <a , 0<y <b }为矩形区域,在x ,y 平面上用两组平行直线x =ih 1, i =0,1,…,N 1, h 1=a /N 1 y =jh 2, j =0,1,…,N 2, h 2=b /N 2将G 剖分为网格区域,见图9-1。
数值模拟偏微分方程的三种方法介绍(有限差分方法、有限元方法、有限体积方法)I.三者简介有限差分方法(Finite Difference Methods)是数值模拟偏微分方程最早采用的方法,至今仍被广泛使用O该方法包括区域剖分和差商代替导数两个步骤。
首先将求解区域划分为差分网格,用有限个网格节点代替连续的求解区域。
其次,利用Taylor级数展开等方法将偏微分方程中的导数项在网格节点上用函数值的差商代替进行离散,从而建立以网格节点上的值为未知量的代数方程组。
该方法是一种直接将微分问题变为代数问题的近似数值解法,数学概念直观,表达简单,是发展较早且十分成熟的数值方法。
差商代替导数后的格式称为有限差分格式,从格式的精度来考虑,有一阶格式、二阶格式和高阶格式。
从差分的空间离散形式来考虑,有中心格式和迎风格式。
对于瞬态方程,考虑时间方向的离散,有显格式、隐格式、交替显隐格式等。
目前常见的差分格式,主要是以上几种格式的组合,不同的组合构成不同的差分格式。
差分方法主要适用于结构网格,网格的大小一般根据问题模型和Courant稳定条件来决定。
有限元方法(Finite Element Methods)的基础是虚位移原理和分片多项式插值。
该方法的构造过程包括以下三个步骤。
首先,利用虚位移原理得到偏微分方程的弱形式,将计算区域划分为有限个互不重叠的单元(三角形、四边形、四面体、六面体等),在每个单元上选择合适的节点作为求解函数的插值点,将偏微分方程中的变量改写成由各变量或其导数的节点值与所选用的分片插值基函数组成的线性表达式,得到微分方程的离散形式。
利用插值函数的局部支集性质及数值积分可以得到未知量的代数方程组。
有限元方法有较完善的理论基础,具有求解区域灵活(复杂区域)、单元类型灵活(适于结构网格和非结构网格)、程序代码通用(数值模拟软件多数基于有限元方法)等特点。
有限元方法最早应用于结构力学,随着计算机的发展已经渗透到计算物理、流体力学与电磁学等各个数值模拟领域。
偏微分方程的求解方法偏微分方程是研究自然现象中具有变化性、互相联系的物理量之间的关系的数学工具。
例如流体力学、电磁学、量子力学等领域中,大量问题都可以用偏微分方程来描述。
因此,研究偏微分方程求解方法是数学领域中一个重要的研究方向。
偏微分方程的一般形式为$$F(x, u, \frac{\partial u}{\partial x}, \frac{\partial^2 u}{\partial x^2}, ..., \frac{\partial^n u}{\partial x^n})=0$$其中,$x$是自变量,$u(x)$是未知函数,$\frac{\partialu}{\partial x}, \frac{\partial^2 u}{\partial x^2}, ..., \frac{\partial^n u}{\partial x^n}$是$u(x)$的各阶导数,$F$是给定的函数。
偏微分方程的求解方法主要有分离变量法、变量代换法、特征线法、有限差分法、有限元法等。
一、分离变量法分离变量法是偏微分方程最常用的求解方法之一。
分离变量法的基本思路是,假设$u(x)$可以表示为几个只与$x$有关的函数的积的形式,通过代入偏微分方程中,再根据对称性和正交性等特征来推导出每个函数的具体形式。
例如,考虑一维热传导方程$$\frac{\partial u}{\partial t}=\alpha\frac{\partial^2 u}{\partialx^2}$$其中,$u(x, t)$表示在位置$x$和时间$t$上的温度分布,$\alpha$为热传导系数。
假设$u(x, t)$可以表示为$$u(x,t)=X(x)T(t)$$将$u(x,t)$代入热传导方程中,得到$$\frac{1}{\alpha}\frac{T'(t)}{T(t)}=\frac{X''(x)}{X(x)}=-\lambda$$其中,$\lambda$为常数。
微分方程数值解差分法微分方程是自然科学和工程技术中广泛使用的工具,它们描述了许多物理过程的动力学行为。
对于复杂的微分方程,解析解往往很难或者不可能得到。
此时我们需要数值解差分法来解决问题。
一、微分方程数值解的方法1.分裂法分裂法是将一个复杂的微分方程分解为多个简单的方程。
例如,将一个偏微分方程分解成几个常微分方程,从而可以方便地使用数值方法计算解。
2.有限差分法有限差分法是一种常见的微分方程数值计算方法。
它将一维或多维的连续函数离散为一系列离散点,然后使用差分方程近似微分方程,最后用迭代法计算数值解。
3.有限元法有限元法是一种广泛使用的数值计算方法,它可以用于求解各种类型的微分方程。
该方法将求解区域分割成多个小区域,然后对每个小区域进行离散化和近似处理。
二、数值解差分法数值解差分法是微分方程数值解的基本方法之一。
它是一种基于差分方程的离散化方法,可以对微分方程进行近似,并将微分方程转化为一个差分方程。
数值解的差分法可以分为前向差分、后向差分和中心差分三种方法。
1.前向差分法前向差分法使用前一时间步的值,计算当前时间步的值。
它的近似误差随着时间步长的增大而增大。
前向差分的公式如下:y_i+1 = y_i + hf_i(x_i,y_i)其中,h是时间步长,f_i是微分方程的左侧。
2.后向差分法后向差分法使用后一时间步的值,计算当前时间步的值。
它的近似误差随着时间步长的增大而减小。
后向差分的公式如下:y_i+1=y_i + hf_i(x_i+1,y_i+1)3.中心差分法中心差分法使用前一时间步和后一时间步的值,计算当前时间步的值。
它的近似误差随着时间步长的增大而增大。
中心差分的公式如下:y_i+1=y_i + 1/2hf_i(x_i,y_i) + 1/2hf_i(x_i+1,y_i+1)三、差分法的优缺点差分法作为微分方程数值解的一种基本方法,具有以下优缺点:1.优点(1)简单易实现:差分法的实现很简单,只需要计算微分方程的离散值和靠近值即可。
--以有限差分法为例偏微分方程数值求解1. 偏微分方程求解问题的描述教材P653[12.1.1]椭圆型教材P653[12.1.2]教材P664[12.2.1]双曲型教材P665[12.2.4]拉普拉斯泊松对流波动教材P684[12.3.1]抛物型教材P685[12.3.6]扩散对流扩散教材P686[12.3.8]二维扩散教材P678[12.2.23]二维对流⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧≥≤≤==≥≤≤==≤≤=>≥≤≤≤≤⎪⎪⎭⎫ ⎝⎛∂∂+∂∂=∂∂0,0, ),(),,(),(),0,(0,0,),(),,(),(),,0(,0,),()0,,(0,0 , 0 , 0 21212222t L x t x v t L x u t x v t x u t L y t y t y L u t y t y u L y x y x y x u b t L y L x y u x u b t u μμϕΩ求解域初值条件边值条件),,(t y x u 未知函数⎪⎪⎪⎪⎪⎭⎪⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧≥<<-==≥<<==≥≤≤-==≥≤≤==≤≤==≤≤≤≤≤≤⎪⎪⎭⎫ ⎝⎛∂∂+∂∂=∂∂0 , 50 , sin 255sin ),(),5,(0 , 50 , 0),(),0,(0 , 50 , 5sin sin 25),(),,5(0 , 50 , 0),(),,0(5,0,0),()0,,( 10000 , 50 , 50 001.022********t x x x t x v t x u t x t x v t x u t y y y t y t y u t y t y t y u y x y x y x u t y x y u x u t u μμϕΩ求解域初值条件边值条件以具体问题为例演示具体的求解过程),,(t y x u 未知函数0x 1x 2x 3x 4x 0y 1y 2y 3y 4y 0t 1t 2t 3t 4t x j jh x =y k kh y =τn t n =xh x 区间的剖分步长τ区间的剖分步长t y h y 区间的剖分步长y x h h h ==0x 1x 2x 3x 4x 0y 1y 2y 3y 4y 0t 1t 2t 3t 4t jh x j =kh y k =τn t n =xh x 区间的剖分步长τ区间的剖分步长t y h y 区间的剖分步长y x h h h ==1x 2x 3x 4x 0y 1y 2y 3y 4y 0t 1t 2t 3t 4t 0x),,(n k j t y x ⎩⎨⎧===4..0 , 4..04..0j k n 0x 1x 2x 3x 4x 0y 1y 2y 3y 4y 0t 1t 2t 3t 4t ),,(211t y x ⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧≤≤≤≤≤≤=Ω100005050),,(t y x t y x ⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧====Ω4..04..04..0),,(j k n t y x n k j n kj0x 1x 2x 3x 4x 0y 1y 2y 3y 4y 0t 1t 2t 3t 4t ⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧≤≤≤≤≤≤=Ω100005050),,(t y x t y x ⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧====Ω4..04..04..0),,(j k n t y x n k j n kj ?),,(),,(=Ω∈t y x t y x u ?),,(),,(=Ω∈nkjn k j t y x n k j t y x u 求解目标求解目标离散化n kju4040====k k j j 或或或边界点:1x 2x 3x 4x 0y 1y 2y 3y4y 0t 1t 2t 3t 4t 0x4040≠≠≠≠k k j j 且且且内点:1x 2x 3x 4x 0y 1y 2y 3y4y 0t 1t 2t 3t 4t 0x1x 2x 3x 4x 0y 1y 2y 3y 4y 0t 1t 2t 3t 4t 0x 5,0,0)0,,(≤≤=y x y x u 初值条件0),,(04..04..00====t y x u uk j k j kj 0kju000u001u002u003u004u1x 2x 3x 4x 0y 1y 2y 3y 4y 0t 1t 2t 3t 4t 0x000u001u002u003u004u1x 2x 3x 4x 0y 1y 2y 3y 4y 0t 1t 2t 3t 4t 0x 010u011u012u013u014u000u 001u002u003u 004u1x 2x 3x 4x 0y 1y 2y 3y 4y 0t 1t 2t 3t 4t 0x 010u011u 012u 013u014u 020u021u022u023u024u000u 001u 002u 003u 004u 1x 2x 3x 4x 0y 1y 2y 3y 4y 0t 1t 2t 3t 4t 0x 010u 011u 012u 013u 014u 020u 021u 022u 023u 024u 030u031u032u033u034u000u 001u 002u 003u 004u 1x 2x 3x 4x 0y 1y 2y 3y 4y 0t 1t 2t 3t 4t 0x 010u 011u 012u 013u 014u 020u 021u 022u 023u 024u 030u 031u 032u033u 034u040u041u042u043u044u的存储设计计算数据0kju 1x 2x 3x 4x 0y 1y 2y 3y 4y 0t 1t 2t 3t 4t 0x 000u 001u 002u 003u 004u 010u 011u 012u 013u 014u 020u 021u 022u 023u 024u 030u 031u 032u033u 034u040u041u042u043u044u1234512345行号列号MATLAB矩阵U0的存储设计计算数据0kju 1x 2x 3x 4x 0y 1y 2y 3y 4y 0t 1t 2t 3t 4t 0x 000u 001u 002u 003u 004u 010u 011u 012u 013u 014u 020u 021u 022u 023u 024u 030u 031u 032u033u 034u040u041u042u043u044u012341234行号列号C 语言矩阵U0的图像计算结果可视化)0,,( :y x u 1x 2x 3x 4x 0y 1y 2y 3y 4y 0t 1t 2t 3t 4t 0x 000u 001u 002u 003u 004u 010u 011u 012u 013u 014u 020u 021u 022u 023u 024u 030u031u 032u 033u034u 040u041u042u043u044u1234512345行号列号MATLAB 矩阵U00kju),,()0,,(0kjk j u y x y x u 上的点4..0,4..0 ),,( :211===k j t y x u u k j kj求步第边值条件11104t x 103t x y 102t x y 101t x y 100t x y 0,50 , 0),,0(≥≤≤=t y t y u 0),,(104..00===t y x u uk k k边值条件1104tx 103t x y 102tx y 101tx y 10t x y 140u 130u 120u 110u 100u 0),,(104..010===t y x u uk k k 0,50 , 0),,0(≥≤≤=t y t y u边值条件25sin )()sin(25),,(2144..014k k k k k y y t y x u u-===14t 14t 14t x 141t x 140t x y 0,50 , 5sin sin 25),,5(2≥≤≤-=t y y y t y u边值条件25sin )()sin(25),,(2144..014k k k k k y y t y x u u-===14t 14t 14tx 141t x 14t x y 144u134u124u 114u104u 0,50 , 5sin sin 25),,5(2≥≤≤-=t y y y t y u边值条件3),,(103..110===t y x u uj j j 0,50 , 0),0,(≥<<=t x t x u 110t x y 120t x y 130t x y边值条件3),,(103..110===t y x u uj j j 110t x y 120tx y 130t x y 101u102u 103u 0,50 , 0),0,(≥<<=t x t x u边值条件411t x 12t x 13t x )sin(255sin )(),,(2143..114j j j j j x x t y x u u-===0,50 , sin 255sin ),5,(2≥<<-=t x x x t x u边值条件411t x 12t x 13t x )sin(255sin )(),,(2143..114j j j j j x x t y x u u-===0,50 , sin 255sin ),5,(2≥<<-=t x x x t x u 141u 142u143u141u 142u 143u 101u 102u 103u 144u 134u 124u 114u 104u 140u 130u 120u 110u 100u1,,+n j k t x y nj k t x y ,,1+nj k t x y ,,1-nj k t x y ,,1-nj k t x y ,,1+nj k t x y ,,策略”1,,+n jk t x y nj k t x y ,,1+n j k tx y ,,1-njk t x y ,,1-nj k t x y,,1+nj kt x y ,,2,1,121,1,22hu u ubhu u ubuujk kj jk j k kj j k kjkj-+-++-++-=-τ策略”1+n kj u nkj u nj k u 1,-nj k u 1,+nj u ,njk u ,1+策略”1+n kjun kjun j k u1,-n j k u 1,+n jk u,1-n jk u,1+2,1,121,1,22hu u ubhu u ubuujk kj jk j k kj j k kjkj-+-++-++-=-τ策略”1+n kjun kjun j k u1,-n j k u 1,+n jk u,1-n jk u,1+2,1,121,1,22hu u ubhu u ubuujk kj jk j k kj j k kjkj-+-++-++-=-τ策略”1+n kjun kjun j k u1,-n j k u 1,+n jk u,1-n jk u,1+2,1,121,1,22hu u ubhu u ubuujk kj jk j k kj j k kjkj-+-++-++-=-τ)(, 22.3.12,41:22h PDE h b +O =−−−→−≤ττ误差估计的解原偏微分方程求出的近似解按显式差分格式当可证收敛并稳定{}1..1,1-=+M j k n kj u {}Mj k nkj u ..0,=目标{}Mj M k n kju 或或或001==+“隐式差1,,+n j k t x y 11,,++n j k t x y 11,,+-n j k t x y 11,,+-n j k t x y 11,,++n j k t x y nj k t x y ,,“隐式差1,,+n jk t x y 11,,++n j k t x y 11,,+-n j k t x y 11,,+-n j k t x y 11,,++n j k t x yn j k t x y ,,“隐式差11,+-n j k u11,++n j k u 1,1++n jk u 1,1+-n jk unjk u ,1+n kj u“隐式差11,+-n j k u11,++n j k u 1,1++n jk u 1,1+-n jk unjk u,1+n kjuτn kjn kjuu-+1=211,1,11,2h uuub n j k n j k n j k +-++++-21,11,1,12h uuub n j k n j k n jk +-++++-+2hbc τ=标准化nj k n j k n j k n j k n j k n j k u u c u c u c u c u c ,1,111,1,11,1,1)41(=∙-∙-++∙-∙-++++++-+-“隐式差11,+-n j k u11,++n j k u 1,1++n jk u 1,1+-n jk un jk u,1+n kjun jk n jk n j k n jk n j k n jk uuc uc uc uc uc ,1,111,1,11,1,1)41(=∙-∙-++∙-∙-++++++-+-2hbc τ=“隐式差11,+-n j k u11,++n j k u 1,1++n jk u 1,1+-n jk un jk u,1+n kjun jk n jk n j k n jk n j k n jk uuc uc uc uc uc ,1,111,1,11,1,1)41(=∙-∙-++∙-∙-++++++-+-2hbc τ=+1kjn“隐式差+1kjn n n n n n n uuc uc uc uc uc 11111)41(=∙-∙-++∙-∙-+++++111+n u110+n u112+n u 121+n u 101+n unu11列差分方程层的内点值基于例如111:+n u t n n n n n n uuc uc uc uc uc 11121112111110101)41(=∙-∙-++∙-∙-+++++列差分方程基于内点值111+n u “隐式差n n n n n n uuc uc uc uc uc 11111)41(=∙-∙-++∙-∙-+++++112+n u111+n u113+n u122+n u102+n u n u 12列差分方程层的内点值基于例如121:+n u t n n n n n n uuc uc uc uc uc 12122113112111102)41(=∙-∙-++∙-∙-+++++列差分方程基于内点值112+n u “隐式差n n n n n n uuc uc uc uc uc 11111)41(=∙-∙-++∙-∙-+++++112+n u111+n u113+n u 122+n u132+n u121+n u131+n u123+n u133+n un 12un 11un 13un 22un 32un 21un 31un 23un 33u线性方程组“隐式差n n n n n n uuc uc uc uc uc 11111)41(=∙-∙-++∙-∙-+++++红色标志方程组的未知量绿色标志方程组的已知量个差分方程列出个内点值层的基于9)3..1,3..1(911==++j k ut n kjn112+n u111+n u113+n u 122+n u132+n u121+n u131+n u123+n u133+n un12u n 11u n 13u n 22u n 32un 21u n 31u n 23u n 33u线性方程组“隐式差n n n n n n uuc uc uc uc uc 11111)41(=∙-∙-++∙-∙-+++++红色标志方程组的未知量绿色标志方程组的已知量个差分方程列出个内点值层的基于9)3..1,3..1(911==++j k ut n kjn )( 32.3.12:2h PDE +O =−−−−→−τ误差估计的解原偏微分方程求出的近似解按隐式差分格式可证收敛并绝对稳定。
第9章 偏微分方程的差分方法含有偏导数的微分方程称为偏微分方程。
由于变量的增多和区域的复杂性,求偏微分方程的精确解一般是不可能的,经常采用数值方法求方程的近似解。
偏微分方程的数值方法种类较多,最常用的方法是差分方法。
差分方法具有格式简单,程序易于实现,计算量小等优点,特别适合于规则区域上偏微分方程的近似求解。
本章将以一些典型的偏微分方程为例,介绍差分方法的基本原理和具体实现方法。
9.1椭圆型方程边值问题的差分方法9.1.1 差分方程的建立最典型的椭圆型方程是Poisson (泊松)方程G y x y x f yux u u ∈=∂∂+∂∂-≡∆-),(),,()(2222 (9.1)G 是x ,y 平面上的有界区域,其边界Γ为分段光滑的闭曲线。
当f (x ,y )≡0时,方程(9.1)称为Laplace(拉普拉斯)方程。
椭圆型方程的定解条件主要有如下三种边界条件第一边值条件 ),(y x u α=Γ (9.2) 第二边值条件),(y x nuβ=∂∂Γ (9.3) 第三边值条件 ),()(y x ku nuγ=+∂∂Γ (9.4) 这里,n 表示Γ上单位外法向,α(x,y ),β(x,y ),γ(x,y )和k (x,y )都是已知的函数,k (x,y )≥0。
满足方程(9.1)和上述三种边值条件之一的光滑函数u (x ,y )称为椭圆型方程边值问题的解。
用差分方法求解偏微分方程,就是要求出精确解u (x ,y )在区域G 的一些离散节点(x i ,y i )上的近似值u i ,j ≈(x i ,y i )。
差分方法的基本思想是,对求解区域G 做网格剖分,将偏微分方程在网格节点上离散化,导出精确解在网格节点上近似值所满足的差分方程,最终通过求解差分方程,通常为一个线性方程组,得到精确解在离散节点上的近似值。
设G ={0<x <a , 0<y <b }为矩形区域,在x ,y 平面上用两组平行直线x =ih 1, i =0,1,…,N 1, h 1=a /N 1 y =jh 2, j =0,1,…,N 2, h 2=b /N 2将G 剖分为网格区域,见图9-1。
h 1,h 2分别称为x 方向和y 方向的剖分步长,网格交点(x i ,y i )称为剖分节点(区域内节点集合记为G h ={(x i ,y i ); (x i ,y i )∈G }),网格线与边界Γ的交点称为边界点,边界点集合记为Γh 。
现在将微分方程(9.1)在每一个内节点(x i ,y i )上进行离散。
在节点(x i ,y i )处,方程(9.1)为h i i i i i i i i G y x y x f y x yuy x x u ∈=∂∂+∂∂-),(),,()],(),([2222 (9.5) 需进一步离散(9.5)中的二阶偏导数。
为简化记号,简记节点(x i ,y i )=(i ,j ),节点函数值u (x i ,y i )=u (i ,j )。
利用一元函数的Taylor 展开公式,推得二阶偏导数的差商表达式)(0)]1,(),(2)1,([1),()(0)],1(),(2),1([1),(222222212122h j i u j i u j i u h j i y u h j i u j i u j i u h j i x u +-+-++=∂∂+-+-++=∂∂代入(9.5)式中,得到方程(9.1)在节点(i ,j )处的离散形式h j i G j i h h f j i u j i u j i u h j i u j i u j i u h ∈++=-+-+--+-+-),(),(0)]1,(),(2)1,([1)],1(),(2),1([12221,2221其中),(,i i j i y x f f =。
舍去高阶小项)(02221h h +,就导出了u (i ,j )的近似值u i ,j 所满足的差分方程h j i j i j i j i j i j i j i G j i f u u u h u u u h ∈=+--+---+-+),(,]2[1]2[1,1,,1,22,1,,121 (9.6) 在节点(i ,j )处方程(9.6)逼近偏微分方程(9.1)的误差为)(2221h h O +,它关于剖分步长是二阶的。
这个误差称为差分方程逼近偏微分方程的截断误差,它的大小将影响近似解的精度。
在差分方程(9.6)中,每一个节点(i ,j )处的方程仅涉及五个节点未知量u i ,j ,u i +1,j ,u i -1,j ,u i ,j +1,u i ,j -1,因此通常称(9.6)式为五点差分格式,当h 1= h 2=h 时,它简化为h j i j i j i j i j i j i G j i f u u u u u h∈=-+++--+-+),(,]4[1,,1,1,,1,12 差分方程(9.6)中,方程个数等于内节点总数,但未知量除内节点值u i ,j ,(i ,j )∈G h 外,还包括边界点值。
例如,点(1,j )处方程就含有边界点未知量u 0,j 。
因此,还要利用给定的边值条件补充上边界点未知量的方程。
对于第一边值条件式(9.2),可直接取u i ,j =α(x i ,y i ), (i ,j )∈Γh (9.7) 对于第三(k =0时为第二)边值条件式(9.4), 以左边界点(1,j )为例,见图9-2, 利用一阶差商公式)(),1(),0(),0(11h O h j u j u j n u +-=∂∂ 则得到边界点(0,j )处的差分方程j j j jj r u k h u u ,0,0,01,1,0=+- (9.8)联立差分方程(9.6)与(9.7)或(9.8)就形成了求解Poisson 方程边值问题的差分方程组,它实质上是一个关于未知量{u i ,j }的线性代数方程组,可采用第2,3章介绍的方法进行求解。
这个方程组的解就称为偏微分方程的差分近似解,简称差分解。
考虑更一般形式的二阶椭圆型方程G y x y x f Eu yu D x u C y u B y x u A x ∈=+∂∂+∂∂+∂∂∂∂+∂∂∂∂-),(),,(])()([(9.9) 其中A (x ,y )≥A m in >0, B (x ,y ) ≥B m in >0, E(x ,y ) ≥0。
引进半节点,12121h x xi i ±=±,22121h y yi i ±=±利用一阶中心差商公式,在节点(i ,j )处可有)(2),1(),1(),()(]),1(),(),(),1([1)()],21)((),21)([(1),)((211211,211,211211h O h j i u j i u j i x u h O h j i u j i u A h j i u j i u A h h O j i x u A j i x u A h j i x u A x j i j i +--+=∂∂+----+=+-∂∂-+∂∂=∂∂∂∂-+对yu y u B y ∂∂∂∂∂∂),(类似处理,就可推得求解方程(9.9)的差分方程 hj i j i j i j i j i j i j i j i j i j i G j i j i f u a u a u a u a u a ∈=-+++---++---+),(),,(][,,1,1,1,1,,1,1,1,1 (9.10)其中⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎨⎧++++=-=+=-=+=-+--+----+-+---+-+ji j i j i j i j i j i j i j i j i j i j i j i j i j i j i j i j i j i E B B h A A h a D h B h a D h B h a C h A h a C h A h a ,21,21,22,21,2121,,221,221,,221,221,,1,2121,1,1,2121,1)()()2()2()2()2( (9.11) 显然,当系数函数A (x ,y )=B (x ,y )=1, C (x ,y )=D (x ,y )=E (x ,y )=0时,椭圆型方程(9.9)就成为Poisson 方程(9.1),而差分方程(9.10)就成为差分方程(9.6)。
容易看出,差分方程(9.10)的截断误差为)(2221h h O +阶。
9.1.2 一般区域的边界条件处理前面已假设G 为矩形区域,现在考虑G 为一般区域情形,这里主要涉及边界条件的处理。
考虑Poisson 方程第一边值问题⎩⎨⎧Γ∈=∈=∆-),(),,(),(),,(y x y x u Gy x y x f u α (9.12) 其中G 可为平面上一般区域,例如为曲边区域。
仍然用两组平行直线:x =x 0+ih 1,y =y 0+jh 2,i ,j =0,±1,…,对区域G 进行矩形网格剖分,见图9-3。
如果一个内节点(i ,j )的四个相邻节点(i +1,j ),(i -1,j ),(i ,j +1)和(i ,j -1)属于Γ⋃=G G ,则称其为正则内点,见图9-3中打“。
”号者;如果一个节点(i ,j )属于G 且不为正则内点,则称其为非正则内点,见图9-3中打“.”号者。
记正则内点集合为hG ',非正则内点集合为h Γ'。
显然,当G 为矩形区域时,h h h hG G Γ=Γ'=',成立。
在正则内点(i ,j )处,完全同矩形区域情形,可建立五点差分格式h j i j i j i j i j i j i j i G j i f u u u h u u u h '∈=+--+---+-+),(,]2[1]2[1,1,,1,22,1,,121 (9.13) 在方程(9.13)中,当(i ,j )点临近边界时,将出现非正则内点上的未知量,因此必须补充非正则内点处的方程。
若非正则内点恰好是边界点,如图9-4中 D 点,则利用边界条件可取u D =α(D)对于不是边界点的非正则内点,如图9-4中B 点,一般可采用如下两种处理方法。
a.直接转移法.取与点B 距离最近的边界点(如图9-4中E 点)上的u 的值作为u (B )的近似值u B ,即u B =u (E)=α(E)直接转移法的优点是简单易行,但精度较低,只为一阶近似。
b .线性插值法.取B 点的两个相邻点(如图9-4中边界点A 和正则内点C 作为插值节点对u (B )进行线性插值)()()()(21h O C u x x x x A u x x x x B u AC AB AC B C +--+--=则得到点B 处的方程 A B C B x x u h A h h u -=+++=δδδαδ,)(111 线性插值法精度较高,为二阶近似。