椭圆型方程的差分格式
椭圆型方程的差分方法
椭圆型方程的差分方法差分方法是一种数值计算方法，使用近似的差商来表示微分方程。
椭圆型方程是一类常见的偏微分方程，具有重要的数学和物理应用。
在本文中，我们将介绍椭圆型方程的差分方法，并讨论其优点和缺点。
一、椭圆型方程的差分近似L[u]=-∂(p∂u/∂x)/∂x-∂(q∂u/∂y)/∂y+r(x,y)u=f(x,y)其中，L[u]是一个偏微分算子，u(x,y)是未知函数，p(x,y)，q(x,y)，r(x,y)，f(x,y)是已知函数。
椭圆型方程的解通常在一个区域Ω上求解。
差分方法的主要思想是用网格来离散化区域Ω，将连续的偏微分方程转化为离散的代数方程。
对于椭圆型方程，我们可以选择矩形网格，其中Ω可以被划分为N*M个小矩形，并且网格的步长为Δx和Δy。
假设我们要在网格点(xi, yj)处求解未知函数的值uij，其中i和j分别表示网格的行索引和列索引。
我们可以使用中心差分法来近似x和y方向的偏导数，从而得到离散形式的椭圆型方程:L[u] ≈ -(p(xi+1/2, yj)(ui+1,j - ui-1,j)/Δx^2 + p(xi,yj+1/2)(ui,j+1 - ui,j-1)/Δy^2) + q(xi,yj)uij = f(xi,yj)其中，p(xi+1/2, yj)和p(xi, yj+1/2)分别表示在(xi+1/2, yj)和(xi, yj+1/2)处的系数。
可以通过有限差分方式计算出这些系数。
将上述公式在每个网格点(xi, yj)处形成一个方程，从而得到一个线性方程组。
通过求解这个线性方程组，我们可以得到网格点上的未知函数值。
二、椭圆型方程差分方法的优点和缺点差分方法是一种简单有效的数值计算方法，具有以下优点：1.可以处理任意形状的区域Ω：差分方法可以适应不规则网格和复杂区域，因此适用于各种几何形状的椭圆型方程求解。
2.数值稳定性：差分方法可以确保数值解的稳定性，避免数值上的不稳定问题。
3.线性时间复杂度：差分方法的计算复杂度通常是线性的，即解方程的时间随着网格点数的增加而线性增加。
matlab用p-r迭代格式求解椭圆型方程的五点差分格式的数值解
matlab用p-r迭代格式求解椭圆型方程的五点差分格式的数值解椭圆型方程是指具有椭圆形状的偏微分方程，如拉普拉斯方程和泊松方程。
求解椭圆型方程的五点差分格式采用p-r迭代格式，其中p是点迭代格式，r是剩余项。
五点差分格式是一种常用的离散化方法，它将偏微分方程离散化为代数方程组。
在五点差分格式中，我们将椭圆型方程的二阶导数离散化为中心差分的形式，使用对角线占主导地位的三对角矩阵来表示离散化后的方程组。
为了使用p-r迭代格式求解椭圆型方程的五点差分格式，我们首先将原方程进行离散化。
考虑一个二维平面上的椭圆型方程：∂²u/∂x² + ∂²u/∂y² = f(x, y)在区域Ω上，Ω为矩形区域，边界Γ为Ω的边界，假设Γ上已知u的值。
我们将Ω划分为网格点，使得Ω中的每个内部节点(xi, yj)与周围的四个节点(xi-1, yj)、(xi+1, yj)、(xi, yj-1)、(xi, yj+1)建立联系。
离散化后，我们得到如下的五点差分方程：(1/hx²)(ui-1,j - 2ui,j + ui+1,j) + (1/hy²)(ui,j-1 - 2ui,j + ui,j+1) = fi,j其中，ui,j为节点(xi, yj)上的解，fi,j为f(xi, yj)的离散值。
hx和hy分别是x和y的网格步长。
使用p-r迭代格式，我们可以得到迭代公式：ui,j⁽ⁿ⁺¹⁾ = (1/(2a))(ui-1,j⁽ⁿ⁾ + ui+1,j⁽ⁿ⁾) + (1/(2b))(ui,j-1⁽ⁿ⁾+ ui,j+1⁽ⁿ⁾) - (1/(2a+2b))fi,j其中，a = (hx/2)²，b = (hy/2)²，n表示迭代次数。
在实际计算中，我们可以采用逐次超松弛（SOR）方法作为p-r迭代的形式。
SOR方法会引入松弛因子ω，更新的迭代公式为：ui,j⁽ⁿ⁺¹⁾ = (1-ω)ui,j⁽ⁿ⁾ + (ω/(2a))(ui-1,j⁽ⁿ⁾ + ui+1,j⁽ⁿ⁾) + (ω/(2b))(ui,j-1⁽ⁿ⁾ + ui,j+1⁽ⁿ⁾) - (ω/(2a+2b))fi,j此外，为了获得数值解，我们需要设定合适的初始猜测值和迭代停止准则。
椭圆型方程
§1
差分逼近的基本概念
考虑二阶微分方程边值问题
d 2u Lu 2 qu f , a x b, dx u (a) , u (b) , (1.1) (1.2)
其中 q，f 为 [ a , b ] 上的连续函数， q 0, , 为给定常数. 将其分成等分，分点为
称
uh 收敛到边值问题的解 u .
对于差分方程
Lhvi fi , i 1, 2,3,L , N 1,
定义1.3
v0 vN 0 , 如果存在与网格 I h 及右端 fh 无关的常数
数 M 和 h0 , 使 || vh || M || f h ||R ,
0 h h0
称差分方程关于右端稳定.
第二章
椭圆形方程的有限差分法
有限差分法和有限元方法是解偏微分方程的两种主要数值
方法.
有限差分法：从定解问题的微分或积分形式出发，用数值 微商或数值积分导出相应的线性代数方程组. 有限元方法：从定解问题的変分形式出发，用RitzGalerkin 方法导出相应的线性代数方程组，但基函数要按
特定方式选取.
取 x(1) x0 a, x(2) x1 , 得
2
(2.9) (2.10)
W (a) W ( x1 ) 2 qudx
d2 du hi 1 hi dx 2 ( p dx ) 12 i
d 3u 2 p O ( h ) dx 3 i
于是得逼近方程 (2.1)~(2.2) 的差分方程：
ui 1 ui ui ui 1 2 p 1 Lhui pi 1 i h h h h i i 1 i 1 i 2 2 i i 1, 2,, N 1 ui 1 ui qiui fi , hi hi 1 u0 , uN
第五章 椭圆型方程的差分方法
数学模型解的存在性解的唯一性解的稳定性解的一些性质解的表达式}适定性数值解存储量计算时间区域方程定解条件回顾：问题的离散第5章椭圆型方程的差分方法§1-3Poisson方程一.区域矩形圆环离散(差分法):网格剖分矩形区域i i i i二.差分格式224412(i,j)(i,j-1)(i,j+1)(i+1,j)(i-1,j)2. 九点差分格式x x2222222h2222121222221212三. 边界条件的处理(矩形区域)xy0(,)i x y I 10(,)x y +J+1(,)i x y )j y(注:四. 差分格式的性质:111 .i i i i i i +-解的存在唯一性与边界条件无关0000002(11002. (,)x y i i i ii ih ji i j h u u x y +→→+−−−→差分方程解的收敛性2||||||2h h hji h h jjji i h i D D D u D D a u u u a D x ∂⋃∂≤+∆插入引理：设是定义在上的函数，那么有max max max 其中为矩形区域的方向的边长.x.3差分格式的稳定性五. 极坐标下的差分格式22+x y注：r∂r(,)i j r θθπR六. 一般区域DDyO x第一类边界条件：T QP δyh第三类边界条件：PQQPnnPQ Q PnnPQR1θu u ux y n u ux y ∂∂∂∂∂∂∂∂∂∂ 2θQRT PS注：没有统一的近似，只要合理就好。
§4变系数方程abxy 矩形区域iP1Q 3Q 2Q 3N4Q 2N4N1N,i jD+-i i i i i i 11。
椭圆型方程差分法
(4)
令
2 1 1 2 1 A 1 1 2 ( N 1)( N 1)
系数矩阵A是不可约对角占优阵 A 0
解存在唯一,或直接求A的特征值。
8
习题：计算矩阵
A=
2 1 1 2 1 1 1 2 ( N 1)( N 1)
x
2． Poisson方程五点差分格式
u f u in
其中
(0, a ) (0, b )
建立目标点： a b y h k x 一方向步长： I 1 ； 一方向步长： J 1
21
得
( xi , y j )
1 i I,1 j J
xi ih, y j jk
返回
2
1） 数值计算是否必要？
T '( x) T '(0) f (u )du
x 0
x
T ( x) T '(0) f ( s )ds du
0
0
u
T '(0) x
T '(0) x
x 0 x
x
x
0
u
0
f ( s)dsdu
s
f ( s )duds
从而得到迭代法：
Mxk 1 Nxk b
xk 1 M 1Nxk M 1b Sxk Tb
(*1)
18
阻尼迭代法 (Damped Iterative Method)
k1 Sxk M1b x k1 (1)xk [M1N(1)I]xk M1b (*2) xk1 x
椭圆型方程的差分解法
椭圆型方程的差分解法1.引言考虑问题①二维Poisson 方程2222(,)u u f x y x y ⎛⎫∂∂-+= ⎪∂∂⎝⎭， (,)x y ∈Ω 其中Ω为2R 中的一个有界区域，其边界Γ为分段光滑曲线。
在Γ上u 满足下列边界条件之一：⑴(,)u x y αΓ=（第一边值条件）， ⑵(,)ux y n βΓ∂=∂（第二边值条件）， ⑶(,)uku x y n γΓ∂+=∂（第三边值条件）， (,),(,),(,),(,),(,)f x y x y x y x y k x y αβγ都是连续函数，0k ≥.2.差分格式将区间[,]a b 作m 等分，记为11()/,,0i h b a m x a ih i m =-=+≤≤；将区间[,]c d 作n 等分，记为22()/,,0i h d c n y c jh j n =-=+≤≤.称1h 为x 方向的步长，2h 为y 方向的步长。
2.1 Poisson 方程五点差分格式参考单如图所示：以(,)i j x y 为中心沿y 方向Taylor 展开：2344111111(,)(,)(,)(,)(,)(,)(),2624i j i j x i j xx i j xxx i j xxxx i j h h h u x y u x y h u x y u x y u x y u x y o h +=+++++①2344111111(,)(,)(,)(,)(,)(,)(),2624i j i j x i j xx i j xxx i j xxxx i j h h h u x y u x y h u x y u x y u x y u x y o h -=-+-++②由① + ② 得：42411111(,)(,)2(,)(,)(,)(),12i j i j i j xx i j xxxx i j h u x y u x y u x y h u x y u x y o h -++=+++整理得到：21121121(,)2(,)(,)(,)(,)(),12i j i j i j xx i j xxxx i j u x y u x y u x y h u x y u x y o h h +--+-=--+③同理，以(,)i j x y 为中心沿y 方向Taylor 展开：21122222(,)2(,)(,)(,)(,)(),12i j i j i j yy i j yyyy i j u x y u x y u x y h u x y u x y o h h +--+-=--+④代入原方程得○5：11112212(,)2(,)(,)(,)2(,)(,)(,),i j i j i j i j i j i j i j ij u x y u x y u x y u x y u x y u x y f x y R h h +-+--+-+⎛⎫-+=+ ⎪⎝⎭得到截断误差：44222212441(,)(,)()(),12ij i j h i j iju u R u x y u x y h h o h o h x y ⎡⎤∂∂=∆-∆=-++=⎢⎥∂∂⎣⎦ 其中u 是原方程光滑解，舍去截断误差得到逼近Poisson 方程的五点差分方程：11112212(,)2(,)(,)(,)2(,)(,)(,),i j i j i j i j i j i j i j u x y u x y u x y u x y u x y u x y f x y h h +-+--+-+⎛⎫-+= ⎪⎝⎭○6 考虑到边值条件(,)(,)u x y x y αΓ=，构成差分格式：11112212(,)2(,)(,)(,)2(,)(,)(,),(,)(,),i j i j i j i j i j i j i j u x y u x y u x y u x y u x y u x y f x y h h u x y x y α+-+-Γ⎧-+-+⎛⎫-+=⎪ ⎪⎨⎝⎭⎪=⎩○72.2 Poisson 方程九点差分格式由上式 ③ + ④ 得：11112212442221244222222122222(,)2(,)(,)(,)2(,)(,)(,)1(,)()12(,)(,)1(,)12i j i j i j i j i j i j h i j i j iji j i j i j u x y u x y u x y u x y u x y u x y u x y h h u u u x y h h o h x y u x y u x y u x y h h x y x y +-+--+-+=+⎡⎤∂∂=∆+++⎢⎥∂∂⎣⎦⎛⎫∂∂⎛⎫∂∂=∆+++- ⎪ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭V 422212222242222212122222(,)()12(,)(,)(,)1(,)()1212i j i j i j i j i j u x y h h o h x y f x y f x y u x y h h f x y h h o h x y x y ∂++∂∂⎛⎫∂∂∂+=--+-+ ⎪ ⎪∂∂∂∂⎝⎭○8 又()41122222211111112212311111(,)(,)2(,)(,)()1[(,)2(,)(,)2(,)2(,)(,)(,)2(,)(,)]()i j xx i j xx i j xx i j i j i j i j i j i j i j i j i j i j u x y u x y u x y u x y o h x y h u x y u x y u x y u x y u x y u x y h h u x y u x y u x y o h +-+++-++-+----∂-+=+∂∂=-+--++-++ 则得到：222222121121112112222221211212122222221112111211()(,)(210)(,)()(,)(210)(,)20()(,)(210)(,)(210)(,)()(,)()(,)i j i j i j i j i j i j i j i j i j h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y ---+--++-+++-++--++-+++-+--+-+2212222241222,12(,)(,)1(,)()12i j i j i j h hf x y f x y f x y h h o h x y ⎛⎫∂∂=--++ ⎪ ⎪∂∂⎝⎭○9 舍去截断误差得到逼近Poisson 方程的九点差分方程○10：()()2212,11,,11,1,11,11,11,122122212(,)[42]121(,)(,),12i j i j i j i j i j i j i j i j i j i j ij xx i j yy i j h h u x y u u u u u u u u u h h f h f x y h f x y -++-+---++-++-∆--+++++++''''=++考虑到边值条件(,)(,)u x y x y αΓ=，构成差分格式○11：()()2212,11,,11,1,11,11,11,122122212(,)[42]121(,)(,),12(,)(,),i j i j i j i j i j i j i j i j i j i j ijxx i j yy i j h h u x y u u u u u u u u u h h f h f x y h f x y u x y x y α-++-+---++-+Γ⎧+-∆--+++++++⎪⎪⎪''''=++⎨⎪⎪=⎪⎩3.格式求解3.1 Poisson 方程五点差分格式记122,1,j j j m j m j u u u u u --⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦M ，0.j n ≤≤矩阵格式改写为：11,11j j j j Du Cu Du f j m -+++=≤≤-，其中2221212222112122221121222112(1)111211112111121112m h h h h h h h C h h h h h h h -⎡⎤⎛⎫+-⎢⎥ ⎪⎝⎭⎢⎥⎢⎥⎛⎫⎢⎥-+- ⎪⎢⎥⎝⎭⎢⎥=⎢⎥⎢⎥⎛⎫⎢⎥-+- ⎪⎢⎥⎝⎭⎢⎥⎛⎫⎢⎥-+ ⎪⎢⎥⎝⎭⎣⎦O OO ，22222222(1)1111m h h D h h -⎡⎤-⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥-⎢⎥⎣⎦O，10212212111(,)(,)(,)(,)1(,)(,)j j j j m j m j m j m f x y x y h f x y f f x y f x y x y h ---⎡⎤+Φ⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥+Φ⎢⎥⎣⎦M ， 可进一步写为：110222211(1)*(1).n n n n n n m u f Du C D u f D C D u f D C D u f Du D C -------⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦M M O O O3.2 Poisson 方程九点差分格式记122,1,j j j m j m j u u u u u --⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦M ，0.j n ≤≤矩阵格式改写为：11,11j j j j Du Cu Du f j m -+++=≤≤-，其中2222121222222212121222222212121222221212(1)20()(210)(210)20()(210)(210)20()(210)(210)20()m h h h h h h h h h h C h h h h h h h h h h -⎡⎤+-⎢⎥-+-⎢⎥⎢⎥=⎢⎥-+-⎢⎥⎢⎥-+⎣⎦O O O， 2222211222222212211222222212211222221221(1)(210)()()(210)()()(210)()()(210)m h h h h h h h h h h D h h h h h h h h h h -⎡⎤--+⎢⎥-+--+⎢⎥⎢⎥=⎢⎥-+--+⎢⎥⎢⎥-+-⎣⎦O，22121022221211(,)(210)(,)(,)(,)(,)(210)(,)j j j j m j m j m j m f x y h h x y f x y f f x y f x y h h x y ---⎡⎤--Φ⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥+-Φ⎣⎦M， 可进一步写为：110222211(1)*(1).n n n n n n m u f Du C D u f D C D u f D C D u f Du D C -------⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦M M O O O4.数值例子4.1 Poisson 方程五点差分格式计算如下问题：22220,01,01,(0,)sin cos ,(2,)(sin cos ),01,(,0),(,1)(sin1cos1),0 1.x x u u x y x y u y y y u y e y y y u x e u x e x ⎛⎫∂∂-+=<<<< ⎪∂∂⎝⎭=+=+≤≤==+<<其精确解为：(,)(sin cos ).x u x y e y y =+,11,1,,1,222222122112112()(,),i j i j i j i j i j i j u u u u u f x y h h h h h h -+-++=++++ 考虑到本例中h1=h2，则有2,11,1,,1,(,),4i j i j i j i j i j i j u u u u h f x y u -+-+++++=利用Gauss-Seidel 迭代方法对k=0,1,2,……,计算112,11,1,,11(,),41,2,....,1;1,2,...., 1.k k k k i j i j i j i j i j k ij u u u u h f x y u i m j n ++--+++++++==-=-表1 部分结点处的精确解和取不同步长时所得的数值解表2 取不同步长时部分结点处数值解的误差绝对值图1 取h=1/4时所得的数值解曲线图2 取h=1/4时所得的误差曲线图3 取h=1/16时所得的数值解曲线图4 取h=1/16时所得的误差曲线图5 取h=1/64时所得的数值解曲线图6 精确解曲线图7 取h=1/64时所得的误差曲线4.2 Poisson 方程九点差分格式计算如下问题：22220,01,01,(0,)sin cos ,(2,)(sin cos ),01,(,0),(,1)(sin1cos1),0 1.x x u u x y x y u y y y u y e y y y u x e u x e x ⎛⎫∂∂-+=<<<< ⎪∂∂⎝⎭=+=+≤≤==+<<其精确解为(,)(sin cos ).x u x y e y y =+222222221212121112122222222121112111211211222211120()(,)12(,)()(,)(102)(,)()(,)()(,)()(,)(102)(,)(102)(,)(10i j i j i j i j i j i j i j i j i j h h u x y h h f x y h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y h ----++++--++=+++-+++++++-+-+2212)(,)i j h u x y +-考虑到本例中h1=h2，则有,11,1,,11,11,11,11,1,4(),20i j i j i j i j i j i j i j i j i j u u u u u u u u u -+-+--++-++-+++++++=利用Gauss-Seidel 迭代方法对k=0,1,2,……,计算1111,11,1,,11,11,11,11,11,4(),201,2,....,1;1,2,...., 1.k k k k k k k k i j i j i j i j i j i j i j i j k i j u u u u u u u u u i m j n ++++-+-+--++-++-++++++++==-=-表1 部分结点处的精确解和取不同步长时所得的数值解表2 取不同步长时部分结点处数值解的误差绝对值表3 取不同步长时部分结点处数值解的最大误差图1 取h=1/4时所得的数值解曲线图2 取h=1/16时所得的数值解曲线图3 取h=1/64时所得的数值解曲线图4 取h=1/4时所得的误差曲线图5 取h=1/16时所得的误差曲线图6 取h=1/64时所得的误差曲线5.结论观察Poisson方程五点格式，方程以较快速度迭代收缩。
椭圆型方程的差分方法
通过实验验证理论分析的正确性。
参数调整
根据误差分析结果调整差分方法的参数。
稳定性分析的实例和结果
结果1
通过误差分析和数值实验，验证了差分方 法的数值稳定性和精度。
A 实例1
一维椭圆型方程的差分方法稳定性 分析。
B
C
D
结果2
通过误差分析和数值实验，验证了差分方 法的数值稳定性和精度，并比较了一维和 二维情况下的误差传播特性。
差分方法在椭圆型方程求解中的优势和局限性
优势
差分方法是求解偏微分方程的一种有效 数值方法，特别适用于大规模计算和并 行计算。它能够模拟偏微分方程的解， 并且具有较高的计算效率和精度。
VS
局限性
差分方法在处理边界条件和复杂几何形状 时可能遇到困难，有时需要引入额外的近 似和假设。此外，差分方法对于某些特殊 类型的偏微分方程可能不适用，或者需要 特殊的处理技巧。
04
差分方法的稳定性分析
稳定性分析的基本概念
数值稳定性
差分方法求解偏微分方程时，数值解对初值 和参数的敏感性。
误差传播
差分方法求解过程中误差的累积和扩散现象。
数值解的精度
差分方法得到的数值解与真实解之间的误差 大小。
稳定性分析的方法和步骤
建立数学模型
将偏微分方程转化为差分方程。
误差分析
计算差分方程的截断误差和全局误差。
差分方法的数学基础
离散化
将连续的函数或过程转换为离散的形式，以便于用数 值方法进行计分方程转化为差分 方程。
稳定性
差分方法的稳定性是指当时间步长趋于无穷小时，差 分方法的解收敛于微分方程的解。
差分方法的实现步骤
建立差分方程
根据微分方程和初边值条件，建立离散化的差 分方程。
椭圆型方程的有限差分法
第4章 椭圆型方程的有限差分法§ 2 一维差分格式1、用积分插值法导出逼近微分方程的差分格式。
, d / du 、 du 上 , Lu=- (p )+r +qu=f,a<x<b, dx dx dxu(a)= a ,u(b)=如果系数p,q,r 以及右端f 光滑，则可用中矩形公式计算得解：考虑在［a,b ］内任一小区间［x ⑴，x (2)］，将上式在此区间上积分得 X ⑵d dU曲 d ；(p(x)d ；)dxW(x ⑴)W(x ⑵)x(2)dU(1) r ——dxxdx x ⑵dU 」 xr dxx (2) x (1) qudx X (2)X (1) qudxX (2)x (1)fdXX ⑵X (1) fdX其中，W(x) p吧特别地，取［X ⑴，X ⑵］为对偶单元［x i 1/2 , x i 1/2 ］，则W (x i 1/2)W(x i 1/2)x1/2dUr ——dxxi 1/2 qudxxi 1/ 2xi 1/ 2 fdx 。
Xi 1/ 2-J..将(1.2)改写成—dxW(x)再沿［x i 1/2, x i 1/2］积分，得U i U iX iW(x) dx , Xi 1p(x)矩形公式，得W 1/2a iU i U ih i a iNdx x1 p(x)］xi 1/ 2又qUdxXi 1/ 2h ih i2：d i U i , d ixi 1/ 2 x q(x)dxX1/2x1/2dU , r ——dx.U i 1 U i 1bi "T-2 h i h i 1x 1/2 xr(x)dxxi1/22ih i h i 1x1/ 2 f (x)dxxi 1/2(1.1) (1.2)利用中(1.3)(1.4)(1.5)(1.6)将(1.3) ~ (1.5)代入1.1),即得微分方程的差分格式a iU i U i 1h*h i h i 1)d i U b^Ui21-(h i h i 1) i则逼近阶为O(h 2)。
椭圆形方程的差分法
【实验结论】（结果）
y =
0
0.2000
0.4000
0.6000
0.8000
1.0000
1.2000
1.4000
1.6000
1.8000
2.0000
2.2000
2.4000
2.6000
2.8000
3.0000
3.2000
3.4000
3.6000
3.8000
4.0000
x0=D(1);
xf=D;
dx=(xf-x0)/M;
x=x0+[0:M]*dx;
dy=(yf-y0)/N;
y=y0+[0:N]'*dy;
%边界条件
for m=1:N+1
u(m,[1,M+1])=[feval(bx0,y(m)),feval(bxf,y(m))];
end
for n=1:M+1
u([1,N+1],n)=[feval(by0,x(n)),feval(byf,x(n))];
end
%边界的平均值作为初始值
bvaver=sum([sum(u(2:N,[1,M+1])),sum(u([1,N+1],2:M))]);
u(2:N,2:M)=bvaver/(2*(M+N-2));
for j=2:M
for i=2:N
u(i,j)=ry*(u(i,j+1)+u(i,j-1))+rx*(u(i+1,j)+u(i-1,j))+rxy*(G(i,j)*u(i,j)-F(i,j));
end
五点差分格式
计算科学系 杨韧
椭圆型方程的五点差分格式
第三章 椭圆型方程的差分格式
一般方程
a(x,
y)
2u x 2
c(x,
y)
2u y 2
d (x,
y)
u x
e(x, y) u f (x, y)u g(x, y) y
Laplace方程
2u x2
2u y 2
u g(x, y) (x, y), h 1 2
n
U7
U8
U9 顶点
解 网格节点如图所示
U4
U5 内点
U6 边界点
U1
U2
U3
椭圆型方程的五点差分格式
矩阵方程为
4 2 0 2 0 0 0 0 0 U1
1 4 1 0 2 0
0
0
0
U
§3.3 混合（Robins）边值条件
2u
x
2
2u y 2
0
(x, y)u (x,
y)
u n
(x,
y)
(x, y) (x, y)
其中(x, y),(x, y) 0。
椭圆型方程的五点差分格式
例1 用五点差分格式求解 Laplace方程
2u y 2
0
(x, y) (x, y) | 0 x 1,0 y 1
u
g(x, y) (x, y)
n
u 表示函数 u 沿着边界的外法线方向导数， n
在正方形的四个顶点上法向量没有定义，取平均值代替。
椭圆型方程的五点差分格式
椭圆 数值解
椭圆型方程的差分方法[五点格式] 考虑拉普拉斯方程的第一边值问题式中为定义在D的边界S上的已知函数.采用正方形网格，记，在节点(i, j)上分别用差商代替，对应的差分方程为（1）或即任一节点(i, j)上的值等于周围相邻节点上解的值的算术平均，这种形式的差分方程称为五点格式，在边界节点上取(2)式中是与节点(i, j)最接近的S上的点.于是得到了以所有内节点上的值为未知量的若干个线性代数方程，由于每一个节点都可列出一个方程，所以未知量的个数与方程的个数都等于节点的总数，于是，可用通常的方法（如高斯消去法）解此线性代数方程组，但当步长不很大时，用高斯消去法将会遇到很大困难，可用下面介绍的其他方法求解.若h 0时，差分方程的解收敛于微分方程的解，则称差分方程为收敛的. 在计算过程中，由于进行四则运算引起舍入误差，每一步计算的舍入误差都会影响以后的计算结果，如果这种影响所产生的计算偏差可以控制，而不至于随着计算次数的增加而无限增大，则称差分方程是稳定的.[迭代法解差分方程] 在五点格式的差分方程中，任意取一组初值{}，只要求它们在边界节点(i , j )上取以已知值，然后用逐次逼近法（也称迭代法）解五点格式：逐次求出{}.当(i+1, j )，(i －1, j )，(i , j －1)，(i , j+1)中有一点是边界节点时，每次迭代时，都要在这一点上取最接近的边界点的值.当n →∞时，收敛于差分方程的解，因此n 充分大时，{}可作差分方程的近似解，迭代次数越多，近似解越接近差分方程的解.[用调节余数法求节点上解的近似值] 以差商代替Δu 时，用节点(i+1, j )，(i－1, j )，(i , j+1)，(i , j －1)上u 的近似值来表示u 在节点(i , j )的值将产生的误差，称此误差为余数，即设在(i , j )上给以改变量，从上式可见将减少4，而其余含有的差分方程中的余数将增加，多次调整的值就可将余数调整到许可的有效数字的范围内，这样可获得各节点上u (x , y )的近似值.这种方法比较简单，特别在对称区域中计算更简捷.例 求Δu =0在内节点A ,B ,C ,D 上解的近似值.设在边界节点1,2,3,4上分别取值为1,2,3,4（图14.8） 解 记u (A)=，点A ,B ,C ,D 的余数分别为图14.8－4+ u B+ u c +5=－4 u B + u D+7=R B－4 u c+ u D+3=R Cu B+ u c－4u D+5=R D以边界节点的边值的算术平均值作为初次近似值，即=u B(0)=u C(0)=u D(0)=2.5则相应的余数为：=0, R B=2, R C= －2, R D=0最大余数为±2.先用δu C=－0.5把R C缩减为零，u C相应地变为2，这时,R D也同时缩减(－0.5)，新余数是=－0.5,R B=2,, R D=－0.5.类似地再变更δu B=0.5，从而u B变为3，则得新余数为.这样便可消去各节点的余数，于是u在各节点的近似值为：=2.5, u B=3, u C=2, u D=2.5现将各次近似值及余数列表如下：[解重调和方程的差分方法] 在矩形D(x0≤x≤x0+a,y0≤y≤y0+a)中考虑重调和方程取步长，引直线族(i, j = 0, 1, 2n)作成一个正方形网格.用差商代替偏导数上式表明了以(x, y)为中心时，u(x, y)的函数值与周围各点函数值的关系，但对于邻近边界节点的点(x, y)，如图14.9中的A，就不能直接使用上式，此时将划分网格的直线族延伸，在延伸线上定出与边界距离为h的点，称这些点为外邻边界节点，如图14.9以A为中心时，点E,C为边界节点，点J,K为E,C的外邻边界节点，用下法补充定义外邻边界节点J处函数的近似值u J，便可应用上面的公式.1 边界条件为图14.9时，定义.2︒边界条件为时，定义.[其他与Δu有关的网格]1︒三角网格（图14.10(a)）取P0(x, y)为中心，它的周围6个邻近节点分别为：则式中，R表示余项.2︒六角网格（图14.10(b)）取P0(x, y)为中心，它的三个邻近节点分别为则.图14.103 极坐标系中的网格（图14.10(c)）取P 0(r,)为中心，它的四个邻近节点分别为而拉普拉斯方程的相应的差分方程为。
第五章 椭圆型方程的差分方法
2 . 差分方程解 的收敛性
hx 0 y ) ( , uij h u x i j 0 y
u f , ( x, y ) D , 定理：若 ( x , y ) D , u , 的解在 D D 上有四阶连续的偏导数， 则五点差分格式
j j j j 1 j j 1 u u u u u u 2 2 j j i 1 i i 1 i i i h ui fi 2 2 hx hy
定理(解的存在唯一性)
j j h u i f i , ( xi , y j ) D h , 差分方程边值问题 j j u ( xi , y j ) D h , i , i 的解如果 存在，则唯一 .
证明：只需证明对应的齐 次方程只有零解
j u h i 0, ( xi , y j ) Dh , j u i 0, ( xi , y j ) Dh , 由极值原 理知0 u i j 0, 即 u i j 0, ( xi , y j ) Dh Dh .
uI+2, j uI , j
i,0ui,0 i,0.
注:
1. 当 =0时, 为第二类边界条件,但四边不能都是第二类边界条件; 2. 这里是中心差商,也可用向前/后差商,但中心差商精度高; 3. 其中u1, j , ui,1, uI 2, j , ui,J 2为虚拟节点的值，需进一步处理， 一般通过方程的差分格式去消掉这些项。
Dh j i Dh j i
( ii )如果uij 满足 h uij 0， ( xi , y j ) Dh，则有 min uij min uij .
Dh Dh
证明: (i) 反证法：设在Dh内存在一点P( xi0 , yi0 )及一 个常数M 使得U ( xi0 , yi0 ) M 并有M ui , ( xi , y j ) Dh
椭圆型方程的差分解法
椭圆型方程的差分解法1.引言考虑问题①二维Poisson 方程2222(,)u u f x y x y ⎛⎫∂∂-+= ⎪∂∂⎝⎭， (,)x y ∈Ω 其中Ω为2R 中的一个有界区域，其边界Γ为分段光滑曲线。
在Γ上u 满足下列边界条件之一：⑴(,)u x y αΓ=（第一边值条件）， ⑵(,)ux y n βΓ∂=∂（第二边值条件）， ⑶(,)uku x y n γΓ∂+=∂（第三边值条件）， (,),(,),(,),(,),(,)f x y x y x y x y k x y αβγ都是连续函数，0k ≥.2.差分格式将区间[,]a b 作m 等分，记为11()/,,0i h b a m x a ih i m =-=+≤≤；将区间[,]c d 作n 等分，记为22()/,,0i h d c n y c jh j n =-=+≤≤.称1h 为x 方向的步长，2h 为y 方向的步长。
2.1 Poisson 方程五点差分格式参考单如图所示：以(,)i j x y 为中心沿y 方向Taylor 展开：41)(),j u y o h +①41)(),j u y o h +②41(),u h21(),o h ③22(),o h ④(,),i j ij f x y R -=+(,),i j f x y -=○6 j+1考虑到边值条件(,)(,)u x y x y αΓ=，构成差分格式：11112212(,)2(,)(,)(,)2(,)(,)(,),(,)(,),i j i j i j i j i j i j i j u x y u x y u x y u x y u x y u x y f x y h h u x y x y α+-+-Γ⎧-+-+⎛⎫-+=⎪ ⎪⎨⎝⎭⎪=⎩○72.2 Poisson 方程九点差分格式由上式 ③ + ④ 得：11112212442221244222222122222(,)2(,)(,)(,)2(,)(,)(,)1(,)()12(,)(,)1(,)12i j i j i j i j i j i j h i j i j iji j i j i j u x y u x y u x y u x y u x y u x y u x y h h u u u x y h h o h x y u x y u x y u x y h h x y x y +-+--+-+=+⎡⎤∂∂=∆+++⎢⎥∂∂⎣⎦⎛⎫∂∂⎛⎫∂∂=∆+++- ⎪ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭422212222242222212122222(,)()12(,)(,)(,)1(,)()1212i j i j i j i j i j u x y h h o h x y f x y f x y u x y h h f x y h h o h x y x y ∂++∂∂⎛⎫∂∂∂+=--+-+ ⎪ ⎪∂∂∂∂⎝⎭○8 又()41122222211111112212311111(,)(,)2(,)(,)()1[(,)2(,)(,)2(,)2(,)(,)(,)2(,)(,)]()i j xx i j xx i j xx i j i j i j i j i j i j i j i j i j i j u x y u x y u x y u x y o h x y h u x y u x y u x y u x y u x y u x y h h u x y u x y u x y o h +-+++-++-+----∂-+=+∂∂=-+--++-++ 则得到：222222121121112112222221211212122222221112111211()(,)(210)(,)()(,)(210)(,)20()(,)(210)(,)(210)(,)()(,)()(,)i j i j i j i j i j i j i j i j i j h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y ---+--++-+++-++--++-+++-+--+-+2212222241222,12(,)(,)1(,)()12i j i j i j h hf x y f x y f x y h h o h x y ⎛⎫∂∂=--++ ⎪ ⎪∂∂⎝⎭○9 舍去截断误差得到逼近Poisson 方程的九点差分方程○10：()()2212,11,,11,1,11,11,11,122122212(,)[42]121(,)(,),12i j i j i j i j i j i j i j i j i j i j ij xx i j yy i j h h u x y u u u u u u u u u h h f h f x y h f x y -++-+---++-++-∆--+++++++''''=++考虑到边值条件(,)(,)u x y x y αΓ=，构成差分格式○11：()()2212,11,,11,1,11,11,11,122122212(,)[42]121(,)(,),12(,)(,),i j i j i j i j i j i j i j i j i j i j ijxx i j yy i j h h u x y u u u u u u u u u h h f h f x y h f x y u x y x y α-++-+---++-+Γ⎧+-∆--+++++++⎪⎪⎪''''=++⎨⎪⎪=⎪⎩3.格式求解3.1 Poisson 方程五点差分格式记122,1,j j j m j m j u u u u u --⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦，0.j n ≤≤ 矩阵格式改写为：11,11j j j j Du Cu Du f j m -+++=≤≤-，其中2221212222112122221121222112(1)111211112111121112m h h h h h h h C h h h h h h h -⎡⎤⎛⎫+-⎢⎥ ⎪⎝⎭⎢⎥⎢⎥⎛⎫⎢⎥-+- ⎪⎢⎥⎝⎭⎢⎥=⎢⎥⎢⎥⎛⎫⎢⎥-+- ⎪⎢⎥⎝⎭⎢⎥⎛⎫⎢⎥-+ ⎪⎢⎥⎝⎭⎣⎦，22222222(1)1111m h h D h h -⎡⎤-⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥-⎢⎥⎣⎦，10212212111(,)(,)(,)(,)1(,)(,)j j j j m j m j m j m f x y x y h f x y f f x y f x y x y h ---⎡⎤+Φ⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥+Φ⎢⎥⎣⎦， 可进一步写为：110222211(1)*(1).n n n n n n m u f Du C D u f D C D u f DC D u f Du D C -------⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦3.2 Poisson 方程九点差分格式记122,1,j j j m j m j u u u u u --⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦，0.j n ≤≤ 矩阵格式改写为：11,11j j j j Du Cu Du f j m -+++=≤≤-，其中2222121222222212121222222212121222221212(1)20()(210)(210)20()(210)(210)20()(210)(210)20()m h h h h h h h h h h C h h h h h h h h h h -⎡⎤+-⎢⎥-+-⎢⎥⎢⎥=⎢⎥-+-⎢⎥⎢⎥-+⎣⎦， 2222211222222212211222222212211222221221(1)(210)()()(210)()()(210)()()(210)m h h h h h h h h h h D h h h h h h h h h h -⎡⎤--+⎢⎥-+--+⎢⎥⎢⎥=⎢⎥-+--+⎢⎥⎢⎥-+-⎣⎦，22121022221211(,)(210)(,)(,)(,)(,)(210)(,)j j j j m j m j m j m f x y h h x y f x y f f x y f x y h h x y ---⎡⎤--Φ⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥+-Φ⎣⎦， 可进一步写为：110222211(1)*(1).n n n n n n m u f Du C Du f D C D u f DC D u f Du D C -------⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦4.数值例子4.1 Poisson 方程五点差分格式计算如下问题：22220,01,01,(0,)sin cos ,(2,)(sin cos ),01,(,0),(,1)(sin1cos1),0 1.x x u u x y x y u y y y u y e y y y u x e u x e x ⎛⎫∂∂-+=<<<< ⎪∂∂⎝⎭=+=+≤≤==+<<其精确解为：(,)(sin cos ).x u x y e y y =+,11,1,,1,222222122112112()(,),i j i j i j i j i j i j u u u u u f x y h h h h h h -+-++=++++ 考虑到本例中h1=h2，则有2,11,1,,1,(,),4i j i j i j i j i j i j u u u u h f x y u -+-+++++=利用Gauss-Seidel 迭代方法对k=0,1,2,……,计算112,11,1,,11(,),41,2,....,1;1,2,...., 1.k k k k i j i j i j i j i j k ij u u u u h f x y u i m j n ++--+++++++==-=-表1 部分结点处的精确解和取不同步长时所得的数值解表2 取不同步长时部分结点处数值解的误差绝对值图1 取h=1/4时所得的数值解曲线图2 取h=1/4时所得的误差曲线图3 取h=1/16时所得的数值解曲线图4 取h=1/16时所得的误差曲线图5 取h=1/64时所得的数值解曲线图6 精确解曲线图7 取h=1/64时所得的误差曲线4.2 Poisson 方程九点差分格式计算如下问题：22220,01,01,(0,)sin cos ,(2,)(sin cos ),01,(,0),(,1)(sin1cos1),0 1.x x u u x y x y u y y y u y e y y y u x e u x e x ⎛⎫∂∂-+=<<<< ⎪∂∂⎝⎭=+=+≤≤==+<<其精确解为(,)(sin cos ).x u x y e y y =+222222221212121112122222222121112111211211222211120()(,)12(,)()(,)(102)(,)()(,)()(,)()(,)(102)(,)(102)(,)(10i j i j i j i j i j i j i j i j i j h h u x y h h f x y h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y h ----++++--++=+++-+++++++-+-+2212)(,)i j h u x y +-考虑到本例中h1=h2，则有,11,1,,11,11,11,11,1,4(),20i j i j i j i j i j i j i j i j i j u u u u u u u u u -+-+--++-++-+++++++=利用Gauss-Seidel 迭代方法对k=0,1,2,……,计算1111,11,1,,11,11,11,11,11,4(),201,2,....,1;1,2,...., 1.k k k k k k k k i j i j i j i j i j i j i j i j k i j u u u u u u u u u i m j n ++++-+-+--++-++-++++++++==-=-表1 部分结点处的精确解和取不同步长时所得的数值解表2 取不同步长时部分结点处数值解的误差绝对值表3 取不同步长时部分结点处数值解的最大误差图1 取h=1/4时所得的数值解曲线图2 取h=1/16时所得的数值解曲线图3 取h=1/64时所得的数值解曲线图4 取h=1/4时所得的误差曲线图5 取h=1/16时所得的误差曲线图6 取h=1/64时所得的误差曲线5.结论观察Poisson方程五点格式，方程以较快速度迭代收缩。
偏微分课程课件9椭圆型方程的有限差分方法I
uij
=
ij
,
(xi，y j ) Dh
其中uij u(xi , y j ) = (xi , y j )=ij , (x, y) D.
例：五差分格式求解
2u 2u
x
2
y2
0
(x, y) D
u(
x,
y)
log
(1
x2
)
y2
( x, y) D
D {( x, y) | 0 x, y 1}
hd e dx u( x, y)
hn d n n! dxn )u( x, y)
ex 1 x x2 x3 xn
1 2 3!
n!
u1 =e u0 , u2 =e u0 , u3 =e u0 , u4 =e u0 , u5 =e u0等
u1 =e u0 , u2 =e u0 , u3 =e u0 , u4 =e u0 , u5 =e u0
从小到大顺利排列
i 1, , J; j 1, , I;
按自然顺序排列网点(i,j)
j 1, i 1, , I; j 2, i 1, , I;
j J , i 1, , I;
定义向量
uh u1,1, , uI ,1, u1,2 , , uI ,2 ,
1 于是差分方程为: h2 Huh g
j 1时
4 1 0
I个
j 2 ...
1 0 0
I个 ...
(1,1)对应H的第一行 11
分析系数矩阵H i 1, , I; j 1, , J;
1 h2 [ui, j1 ui1, j 4ui, j ui1, j ui, j1 ] fi, j
对于第二个结点(2,1),
1 h2 [u2,0 u1,1 4u2,1 u3,1 u2,2 ] f2,1
椭圆方程差分格式
由Guass公式有
D
D
u u udxdy D n ds ( l l l ) n ds D l1 2 3 4
从上面的公式把边界积 分转化为在四条边上的 积分， 于是积分分成 段进行，也即转化为定 4 积分。
对于在边l1上，因为在矩形边1上的外法向就是 的负 l y 方向，弧长的微分 dx，于是对此定积分用矩 ds 形 公式近似计算，并且用 差商代替微商，得到
u u x (b( x, y) x )dxdy [b( xi , y j 12 ) x ( xi , y j 12 ) Dij u b( xi , y 1 ) ( xi , y 1 )]h j j 2 x 2
c( x, y)udxdy c u hk, f ( x, y)dxdy f
１：直接差分方法
1 1 (aij xuij ) 2 y (aij y uij ) cijuij f ij 2 x h h
２：有限体积法（积分差分方法）
u x (a( x, y) x )dxdy Dij
y
u u [a( x 1 , y ) ( x 1 , y ) a( x 1 , y ) ( x 1 , y )]dy i i x i 2 x i 2 2 2 y 1
j 2
j
1 2
对上面定积分利用梯形 公式有
u u x (a( x, y) x )dxdy [a( xi 12 , y j ) x ( xi 12 , y j ) Dij u a( x 1 , y j ) ( x 1 , y j )]k i x i 2 2
1.1:五点差分格式
对于poisson方程，考虑在内部节点 xi，y j）取值， （ 于是有 2u 2u [ 2 ]ij [ 2 ]ij [ f ( x, y )]ij x y
偏微(13)椭圆型方程的差分方法
S1 e u0 e u0 e u0 e u0
1 4 2 4u0 h u0 h 2 D 4 u0 h6 12
2
1 4 2 S2 4u0 2h u0 h 4 D 4 u0 h6 6
其中D x , y | 0 x , y 1 .
1 取特殊的网格 h k 。 此时网格点分布见图5.3， 3
1.1 五点差分格式
在内点P , P2 , P3和P4 1 上用差分格式（ ）， 1.6
在其余点，即边界点 取边界条件
1 x 2 y 2 a x, y log
1 u 1 u r u r r r 2 2 f r , r r
2
(1.8)
r
y x y , tan . 域0 r ,0 2 . x
2 2
方程( 8) 1. 的系数当r 0时具有奇异性，因此为了 选出我们感兴趣的解，需补充u在r 0处有界的条 件，可设u满足
1 u xi h, y j 2u xi , y j u xi h, y j （1.2） h2 Nhomakorabea
2u 2u u 2 2 f ( x , y（ ） ) 1.1 x y
1 u x i , y j k 2u x i , y j u x i , y j k 2 k （1.3） 2 2 4 4 u k 2 4 u xi ,1 4 u xi ,2 y y ij 24 y
h uij ui 1, j 2uij ui 1, j h2 ui , j 1 2uij ui , j 1 k2 f ij ,
椭圆型方程的差分格式
在四个顶点上，有
4U0,0 −2U1,0 −2U0,1 = 4hg0,0
4U0,M −2U1,M −2U0,M−1 = 4hg0,M 4UM,0 −2UM,1 −2UM−1,0 =4hgM,0 4U M ,M − 2U M −1,M − 2U M ,M −1 = 4hgM ,M
⎢⎡− ⎢
⎢⎣⎡1+
1 2
h(
p0
+
q0
)⎥⎦⎤
1 2
⎤ ⎥ ⎥
⎢
1
⎢
2
− (2 + hq1 )
1 2
⎥ ⎥
E0 = ⎢
⎥
⎢ ⎢ ⎢
1 2
− (2 + hqM−2 )
⎤
⎢⎢−1 4 −1
⎥ ⎥
B=⎢
⎥
⎢ ⎢
−1 4 −1⎥⎥
⎢⎣
−1 4 ⎥⎦
3.2 Neumann边值问题的差分模拟
现在我们考虑Laplace方程Neumann边值问题，即
⎪⎪⎩⎪⎪⎨⎧∂∂∂∂unx2u2|∂+Ω=∂∂y2gu2(x=, y0) (x, y)∈Ω;Ω={(x, y)| 0< x <1,0< y <1}
其中矩阵A为M 2 阶对称方阵。
⎡E0 K
⎤
⎢ ⎢
K
E1
K
⎥ ⎥
⎢ A=⎢
K E2 K
⎥ ⎥
⎢
⎥
⎢ ⎢
K EM −2
K⎥ ⎥
⎢⎣
K EM −1 ⎥⎦
⎡−(2+ pm) 1
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
⎧x = 1;0 ≤ y ≤ 1
⎨ ⎩
y
=
1;0
≤
x
≤
1
p, q, f0 , f1, g 是给定的函数。
在求解区域Ω内由逼近Laplace方程的五点差分公式
( ) 1
h2
δx2
+δ
2 y
Ul,m
=
0，l,
m
=
1,2,…,
M
−1
给出函数u在结点(lh,mh)的近似值Ulm所满足的差分方程。
对于在x=0上的结点(0,mh)，应用边值条件
其中，ql = q(lh), f1,l = f1(lh) 。
(3.23)
在原点(0,0)上，两边值条件相遇，则
U−1,0 +U1,0 +U0,−1 +U0,1 − 4U0,0 = 0
U1,0 −U−1,0 −2hp0U0,0 = 2hf0,0
U 0,1 − U 0,−1 − 2hq0U 0,0 = 2hf1,0
( ) 2U1,m +U0,m+1 +U0,m−1 − 4 + 2hpm U0,m = 2hf0,m m = 1,2, , M −1
相似地对于y=0上的结点(lh,0)，我们有 ( ) 2Ul,1 +Ul+1,0 +Ul−1,0 − 4 + 2hql Ul,0 = 2hf1,l l =1,2, , M −1
∂u ∂x
−
p(y)u
=
f0(y)
的差分模拟和五点差分公式，即
U−1,m +U1,m +U0,m+1 +U0,m−1 −4U0,m =0
( ) ( ) U1,m −U −1,m 2h − pmU0,m = f0,m pm = p(mh); f0,m = f0 (mh) （3.22）
消去U −1,m ，得
Ul−1,m +Ul+1,m +Ul,m+1 +Ul,m−1 −4Ul,m =0 l, m = 1,2,…, M −1
( ) Ul,1 +Ul+1,0 2 +Ul−1,0 2 − 2 + hql Ul,0 = hf1,l l =1,2, , M −1
( ) U1,m +U0,m+1 2 +U0,m−1 2 − 2 + hpm U0,m = hf0,m m = 1,2, , M −1
⎤
⎢⎢−1 4 −1
⎥ ⎥
B=⎢
⎥
⎢ ⎢
−1 4 −1⎥⎥
⎢⎣
−1 4 ⎥⎦
3.2 Neumann边值问题的差分模拟
现在我们考虑Laplace方程Neumann边值问题，即
⎪⎪⎩⎪⎪⎨⎧∂∂∂∂unx2u2|∂+Ω=∂∂y2gu2(x=, y0) (x, y)∈Ω;Ω={(x, y)| 0< x <1,0< y <1}
=0
它具有截断误差:
1 12
h 2 ⎜⎜⎝⎛
∂ 4u ∂x 4
+
∂ 4u ∂y 4
⎟⎟⎠⎞ l ,m
+
我们引进记号◇，有
◇Ul,m
=
1 h2
(Ul+1,m
+Ul−1,m
+Ul,m+1
+
Ul,m−1
−
4Ul,m
)
因此差分方程(3.6)即◇Ul,m =0。
（3.6） (3.7)
如图3.1所示
在区域Ω的每一内部结点(l,m)上
(3.10)
∂u 表示函数u沿着边界的外法线方向导数。在正方
形的∂n ∂四Ω 个顶点上法向没有定义，事实上，g(x,y)在那里
将不连续，以后将取平均值作为不连续点上的值的定
义。
Neumann边值问题(3.10)的解存在，仅当
∫ g(x, y)dl＝0 ∂Ω
且除了一个任意常数外，解唯一。因为容易看到，如 果u(x,y)是式(3.10)的解，于是， u(x,y)+C(C是一个任 意常数)也是其解。为了唯一性，需要规定u(x,y)在区域 中某一点上的值。
• 3.7 一般二阶线性椭圆型方程差分逼近及其性质 研究
设 Ω 是平面中的具有边界的一个有界区域，本章 考虑如下椭圆型方程的差分解法：
a (x, y ) ∂ 2u
∂x 2
+
2b(x, y ) ∂ 2u
∂x∂y
+ c(x, y ) ∂ 2u
∂y 2
=
d
⎜⎜⎝⎛
x,
y,
u,
∂u ∂x
,
∂u ∂y
⎟⎟⎠⎞
⎢ g8 ⎥ ⎢⎣2g9 ⎥⎦
(3.19)
或简写成 AU=2hg。 显然A是一奇异矩阵。
3.3 混合边值问题
在xy平面的某区域Ω中，未知函数u满足Laplace方 程，将边界 ∂Ω分成若干弧段，要求u在每一弧段上满足 不同类型的边界条件。讨论此类定解问题的差分模拟。
例如，求解如下定解问题：
⎧ ∂2u
,U T M−1,M−1
单位正方形中的内部结点上的 (M −1)2 个线性方程
(3.8)写成矩阵形式为
AU=K
（3.9）
其中，A是 (M −1)2 阶方阵
⎡ B −I
⎤
⎢⎢− I B − I
⎥ ⎥
A=⎢
⎥
⎢ ⎢
−I
B
−
I
⎥ ⎥
⎢⎣
− I B ⎥⎦
I 是(M-1)阶单位方阵；B是(M-1)阶方阵。
⎡ 4 −1
⎡4 −2
⎤
⎢ ⎢
−
1
4
−1
⎥ ⎥
B=⎢
⎥
⎢ ⎢
− 1 4 − 1⎥⎥
⎢⎣
− 2 4 ⎥⎦
方程组(3.18)中的向量U和g由以下给出：
[ ] U=U0,0,U1,0, ,UM,0;U0,1, ,UM,1; ;U0,M, ,UM,M T [ g = 2 g0,0 , g1,0 , , g M −1,0 ,2 g M ,0 ; g0,1,0, ,0, g M ,1;
（3.15） （3.16） （3.17）
由此，正方形 0 ≤ x, y ≤ 1区域的 (M +1)2 个结点上差分 方程解Ulm 满足线性方程组
AU=2hg
(3.18)
这里A是(M +1)2阶方阵
⎡ B − 2I
⎤
⎢ ⎢
−
I
B
−I
⎥ ⎥
A=⎢
⎥
⎢ ⎢
−I
B
−
I
⎥ ⎥
⎢⎣
− 2 I B ⎥⎦
I是(M+1)阶单位方阵；B 是如下(M+1)的阶方阵：
⎨ ⎪
∂u
= g (x , y )
⎪⎩ ∂ n ∂ Ω
解 令h=1/2，应用图3.2中结点次序，则方程(3.18)为
⎡ 4 − 2 0 − 2 0 0 0 0 0 ⎤⎡U1 ⎤
⎢⎢−1 4 −1 0 − 2 0
0
0
0
⎥ ⎥
⎢⎢U
2
⎥ ⎥
⎡2g1 ⎤
⎢ ⎢
g2
⎥ ⎥
⎢0 −2 4 ⎢⎢−1 0 0
0 0 −2 0 0 4 − 2 0 −1 0
代入式(3.13)，则
4U0,m −2U1,m −U0,m−1 −U0,m+1 =2hg0,m m=1, ,M−1
（3.14）
同理，在x=1,y=0,y=1时分别有
4UM,m −2UM−1,m −UM,m−1 −UM,m+1 = 2hgM,m m=1, ,M −1 4Ul,0 − 2Ul,1 −Ul−1,0 −Ul+1,0 = 2hgl,0 l = 1, , M −1
消去 和 ,则 U−1,0
U 0,−1
( ) ( ) 2U1,0 + 2U0,1 − 4 + 2hp0 + 2hq0 U0,0 = 2h f0,0 + f1,0
(3.24)
且对l=0和m=0上成立的方程(3.22)，(3.23)用1/2乘 之，对l=m=0上的方程(3.24)用1/4乘之。
这样在整个计算区域及相应边界网格点上建立了差 分方程：
Neumann边值问题的差分模拟
先在区域Ω中给定一个正方形网格区域，步长为
h，Mh=1,于是必须确定解的结点为（M +1)2个，结点上
的差分方程的解为
U l,m (0 ≤ l, m ≤ M )
（3.12）
在内点上，Laplace方程由差分方程(3.6)代替：
( ) 1
h2
δ
2 x
+
δ
2 y
U l,m
0 0
⎥ ⎥ ⎥
⎢U ⎢⎢U
3 4
⎥ ⎥ ⎥
⎢⎢⎢2gg43
⎥ ⎥ ⎥
⎢0 ⎢
−1 0
−1
4
−1 0
−1
0
⎥ ⎥
⎢U ⎢
5
⎥ ⎥
=
2h⎢ ⎢
0
⎥ ⎥
⎢ 0 0 −1 0 − 2 4 0 0 −1⎥⎢U 6 ⎥ ⎢ g6 ⎥
⎢ ⎢
0
0
0 −2 0
0
4
−2
0
⎥ ⎥
⎢⎢U
7
⎥ ⎥
⎢⎢2
g
7
⎥ ⎥
⎢ 0 0 0 0 − 2 0 −1 4 −1⎥⎢U8 ⎥ ⎢⎣ 0 0 0 0 0 − 2 0 − 2 4 ⎥⎦⎢⎣U 9 ⎥⎦
] ; g0,M −1,0, ,0, gM ,M −1;2g0,M , g1,M , , gM −1,M ,2gM ,M T