解一维二阶双曲方程的一个新型四阶紧致差分格式
4-迎风格式
§4. 迎风格式 4.1 迎风原则先考虑线化的模型方程 - 对流方程0u u a t x抖+=抖 其中 a 为常数。
迎风原则的示意图4.2 迎风格式向后差分格式和向前差分格式都可以看成是对中心差分格式的修改111112 22n n n nn n nj jj j j j j u u u u u u u aatxx++-+----++=D D D （用于 0a > ）11111222n nn nn n nj jj j j j j u u u u u u u aatxx++-+----++=-D D D （用于 0a < ）不仅如此，经过这样的改写，我们发现它们还可以统一起来，写成11111222n nn nn n nj jj j j j j u u u u u u u aatxx++-+----++=D D D如果再经过进一步的改写11111022n nn n n nn n n nj jj j j j j j j j u u u u u u u u u u aatxx++-+---+---++-=D D D并重新组合格式中的各项，就得到111022n nn nn nj jj j j ja a a a u u u u u u txx+-++----++=D D D记 2a a a ++=，2a a a --=，则上述格式可写成1110n nn nn nj j j j j ju u u u u u a a t xx+-++----++=D D D注意到，无论 0a > 还是 0a < ，都有02a a a ++=> ， 02a a a --=<而且22a a a a a a a +-+-+=+=于是我们可以对上面得到的格式做出新的解释：将 a 分解成了“正”和“负”两个部分a a a +-=+ ， 0a +> ， 0a -<对流项也相应地分裂成了两项。
计算电磁学-第4章-有限差分法
同样对微分方程的解y(x)在点(xn,yn)进行泰勒展开
yn1 yn hf ( xn , yn )
1 ' 2 1 '' 3 y ( xn 1 ) y ( xn ) f n h f n h f n h 2! 3!
比较上面两式，只要它们前面项的系数尽可能多的相等，就 保证了截断精度。
1、差分与差商
用差分代替微分，是有限差分法的基本出发点。 这一点由微分原理保证的，当自变量的差分趋于 零时，差分变成微分
f ( x) f ( x h) f ( x)， h x
df f ( x) f ( x) lim dx x 0 x
'
f ( x) f ( x h) f ( x) f ( x) x h
龙格-库塔法
选取α、β、ω系数，使两式项的系数相等
1 fn , 2 f , 3 f , 4 f ,
' n '' n ''' n
如果该关系式能够一直维持到第m阶仍能成立， 但m＋1阶不再成立，就称为m阶龙格-库塔法
cem@
cem@
cem@
cem@
cem@
cem@
cem@
CST粒子仿真
Pierce Gun
ＭＡＧＩＣ
cem@
dy f ( x, y ) dx y x x 0 y0
y( x) y0 f (t , y(t )dt
x0
x
欧拉近似法在函数图上用阶梯的折线代替曲线
f(x) y(x)
yn+1 yn y(x n+1)1) f(n+
计算地球流体力学：第2讲 有限差分格式的构造
y
fu
0
(2.9)
t
u
x
v
y
0
其中 u，v 和 分别为纬向风，经向风和位势高度，
f 为科氏参数。
10
发展方程的算子形式
计算地球流体力学中所要求解的发展方程一般都可以 化为如下算子形式：
F L F G t
(2.10)
其中， L L (F, x,t)是一个线性或非线性算子，F F(x,t) 是
u t
a
u x
0
(0 x l, 0 t T ), (a)
u(x, 0) (x)
(0 x l),
(b)
u(0,t) 1(t),u(l,t) 2 (t) (0 t T ).
(c)
(2.4)
第一步，先对求解区域 {(x,t) | 0 x l,0 t T} 作如图 2.1 所示的网格剖分。为简单起见，假设时空都为等距剖分，格 距分别为 t 和 x 。
u t
n i
O(t 2
)
(2.12)
37
(2.11)
u(xi
,
tn 1 )
u ( xi
,
tn
)
u t
n i
t
1
2!
2u t 2
n
i
t
2
O
t 3
u(xi
,
tn1
)
u ( xi
,
tn
)
u t
n i
t
1
2!
2u t 2
n
i
t
2
O
t 3
38
u(
xi
,
tn
1
)
u(
偏微分方程数值解法(1)
第十章 偏微分方程数值解法一、 典型的偏微分方程介绍 1．椭圆型方程 科学技术中经常遇到一些重要的、典型的偏微分方程。
在研究有热源稳定状态下的热传导，有固定外力作用下薄膜的平衡问题时，都会遇到Poisson 方程D y x y x f yux u ∈=∂∂+∂∂),(),(2222（10.1）其中D 表示平面区域。
特别在没有热源或没有外力时，就得到Laplace 方程02222=∂∂+∂∂y ux u （10.2）此外，当研究不可压缩理想流体无旋流动的速度势以及静电场的电位等，也会遇到（10.1）或（10.2）类型的方程。
2．抛物型方程 在研究热传导过程、气体扩散现象、电磁场的传播等问题中以及在统计物理、概率论和重子力学中，经常遇到抛物型方程。
这类方程中最简单、最典型的是热传导方程。
L x t xu a t u <<>=∂∂-∂∂0,0,022（10.3）其中a 是常数。
它表示长度为L 的细杆内，物体温度分布的规律。
3．双曲型方程 在研究波的传播、物体的振动时，常遇到双曲型方程。
这类方程中最简单、最典型的是波动方程L x t xu a t u <<>=∂∂-∂∂0,0,022222（10.4）它表示长度为L 的弦振动的规律。
二、定解问题偏微分方程（10.1）~（10.4）是描述物理过程的普遍规律的。
要使它们刻划某一特定的物理过程，必须给出附加条件。
把决定方程唯一解所必须给定的初始条件和边界条件叫做定解条件。
定解条件由实际问题提出。
对方程（10.3）来说，初始条件的提法应为)()0,(x f x u =，其中f (x )为已知函数，它表示物体在初始状态下温度分布是已知的。
边界条件的提法应为物体在端点的温度分布为已知，即⎩⎨⎧≥==0)(),()(),0(t t t L u t t u ψϕ （10.5）其中ϕ（t ）和ψ（t ）为已知函数。
对（10.4）来说，边界条件的提法和（10.5）形式一样，它表示弦在两端振动规律为已知。
求解波动方程的时空四阶算法
数，且 φ（ａ）＝α（０）、φ（ｂ）＝β（０）、ψ（ａ）＝α′（０）、
ψ（ｂ）＝β′（０）。 本文假定方程（１）—（３）的解满足
ｕ（ｘ，ｔ）∈ Ｃ６，４（［ａ，ｂ］）
（４）
并且ｆ（ｘ，ｔ）在解的ε０ 邻 域 内 关 于 第 一、第 二 分 量
有一阶连续导数，其中ε０ 是一个正数。
对于 二 阶 线 性 边 值 问 题，常 用 的 数 值 方 法 有 打
双曲型偏微分方程在实际问题中有着广泛的应 用，但一般情况下，直 接 分 析 和 求 解 很 困 难，因 此 常 常通过对数值格式稳定性进行间接分析得到数值格 式的精度。本文利 用 时 空 紧 算 子，推 导 了 一 类 新 的 差分格式，利用 对 角 占 优 矩 阵 性 质 证 ［５］ 明 差 分 格 式 解的存在 唯 一 性；另 外 利 用 Ｆｏｕｒｉｅｒ级 数 分 析 了 该 格式的条件稳定 性，并 利 用 Ｌａｘ－Ｍｉｌｇｒａｍ 定 理 以 及 相容性条件得到了 数 值 格 式 的 收 敛 性，为 数 值 求 解 双曲型方程提供一种有效解法。
靶法和有限差分法等。很多学者都在此领域进行了
深入广泛的 研 究。 例 如，孙 志 忠 以 ［１］ 波 动 方 程 为 模
型，对第一边值问题 分 别 建 立 了 时 空 方 向 均 有 二 阶
精度的显格式和 隐 格 式。 随 后，他 又 建 立 了 时 间 二
阶空间四阶的紧致差分格式。刘相国等 利 ［２］ 用有 限
关 键 词 ：波 动 方 程 ；紧 差 分 格 式 ；可 解 性 ；收 敛 性 ；稳 定 性 中图分类号：Ｏ２４４ 文献标志码：Ａ 文章编号：１６７３－３８５１ （２０１９）０３－０２４９－０６
Ｔｉｍｅ－ｓｐａｃｅ ｆｏｕｒｔｈ ｏｒｄｅｒ ａｌｇｏｒｉｔｈｍ ｆｏｒ ｓｏｌｖｉｎｇ ｗａｖｅ ｅｑｕａｔｉｏｎｓ
一维Euler方程的特征有限体积格式
文章编号 :100020887 (2009) 0320291210 Ζ 应用数学和力学编委会 , I SSN 100020887Ξ一维 Euler 方程的特征有限体积格式郭 彦 , 刘儒勋(中国科学技术大学 数学系 ,合肥 230026)(周哲玮推荐)摘要 : 提出了一种用于求解一维标量方程和无粘 Eu ler 方程组的高阶有限体积格式 1 其中时间 离散采用 S imp son 数值积分公式从而实现时间上的高阶 1 利用特征线理论得到网格节点在各个时 间层沿着特征线的位置 ,而积分公式中的节点值通过三阶和五阶的中心加权本质无震荡重构得到 1 最后 ,给出了几个数值算例验证此方法的高精度和收敛性以及捕获激波的能力 1关 键 词 : 双曲方程 ; 有限体积方法 ; 特征理论 ; WE NO 重构 ; Rung e 2Ku tta 方法 中图分类号 : O241. 8 ;O352文献标识码 : A引 言当今求解流体动力学问题需要一些高精度的计算方法 1 最近几年已经提出很多高阶格式求解守恒型双曲方程组 ,这些格式一般都是基于高阶重构数值流通量来构造高阶格式 ,而通过 这种方法构造的高阶格式在光滑区域能够达到很高精度 ,但是在出现间断区域就会出现明显的数值震荡 1 高阶无震荡格式中具有代表性的格式是本质无震荡格式 ( E NO ) 1 22,加权本质无 震荡格式 (WE NO ) 3 和中心加权本质无震荡格式 ( C WE NO ) 4 25 1 还有其它高阶格式如最近所提出的 RK D G 有限元格式 6 ,激波捕获差分格式紧7 ,紧致格式8 ,分段有理方法 ( PR M ) 9和 CIP/ MM FV M 格式10 1 其中 ,CIP/ MM FV M 格式是通过对一维 E uler 守恒律方程运用特征线理 论 ,并沿着特征曲线求解 R iem ann 不变量的半 Lagrange 数值解 1 该格式主要利用节点值 ,单元 平均值和一阶导数值构建 3 次多项式 ,从而利用特征线理论所求出的节点位置代入所求的三 阶多项式就可以得到各个变量在节点的值 1 从格式分析及其数值结果中可以看到该格式具有 三阶精度 1本文对 E uler 守恒律方程组提出一种新的高阶有限体积格式 1 采用对时间积分进行 Sim p 2 son 积分公式离散得到时间上的高精度 1 我们运用 C WE NO 重构 4 25 得到的多项式和由特征方 程求得的节点位置得到有限体积格式中的各时间层沿特征线的变量的节点值 1 这种处理结合 特征线方法的高分辨率性质 、C W E NO 格式的高精度和本质无震荡性质 1 最后给出一些经典的 数值算例验证该格式的数值精度和格式的性质 1本文如下安排各节内容 : 第 1 节中 ,对一维的标量方程和欧拉方程组分别构造所提的格 式 ;第 2 节将该格式应用于一些经典算例 ,从算例中可以看出格式的高精度收敛性及其捕获激Ξ 收稿日期 : 2008207204 ; 修订日期 : 2009202212基金项目 : 国家自然科学基金资助项目 (10771134)作者简介 : 郭彦 (1982 —) ,男 ,安徽人 ,博士 (联系人 . E 2mail : g ysx @mail . ustc . edu. cn) .291292郭 彦 刘 儒 勋波的能力 1 最后对本文做简短的总结 11 数 值 方 法1 . 1 标量守恒方程本节考虑如下标量守恒方程 :u t + f ( u ) x = 0 ,其中 , u ( x , 0) = u 0 ( x ) 1为近似求解方程 (1) ,分别离散空间和时间并取步长 h (1)= Δx 和Δt , 记空间网格节点为 x i= x i - 1/ 2 , x i +1/ 2 ] 为空间剖分单元 ,, N , 时间层为 t n +1 = t n +Δt , 同时记 I i = i × h , i = 0 , x i +1/ 21 h∫x¯u in ≈ u ( x , t n) d x , i = 1 ,, N 为变量在时间 t n处的单元平均值 , 记 u n≈ u ( x , t n) , i =i i i - 1/ 2, N 为节点值 1 对方程 (1) 在单元片 I i ×[ t n , t n +1 ] 上积分可以得到如下公式 :0 ,n +1t1 h∫tn¯u n +1 n= ¯u i - , t ) ) ]d t 1 (2)[ f u (( xi +1/ 2 , t) ) - f ( u ( x ii - 1/ 2 由各变量在时间 t n 层处的单元平均值通过该方程求得在时间 t n +1 处的值 1 为得到方程 (2) 的 全离散形式 ,下面就要对时间积分进行离散 11 . 1 . 1 时间积分本文利用 Sim p s on 积分公式对方程 (2) 进行时间离散实现时间上的高精度3d t d x 6^ ^ ¯u n +1nn n= ¯u i- N l [ f ( u ( x i +1/ 2 , t + βl d t ) ) - f ( u ( x i - 1/ 2 , t + βl d t ) ) , (3)i l = 1^其中 , N = ( 1/ 6 , 2/ 3 , 1/ 6) T 为权 ,β = ( 0 , 1/ 2 , 1) T 为积分节点 1 u 为所求变量的节点值 11. 1 . 2 网格点位置类似 CIP/ MM FV M 10,我们运用半 Lagrange 方法求解如下特征方程 ,从而得到节点沿特征线的位移 :d x ′ = f ( u ) 1 (4)d t求解该方程有很多方法 1 本文分别采用四阶 Runge 2K utta 方法和中点方法求解该特征方程 1四阶 Runge 2K utta 方法求解单元边界点沿特征线在时间 t n +1 处的位置 `x j +1/ 2 可以表示为4x n +1 = x n+ d t6 bkg ,( k )(5)k = 1其中 , b = ( 1/ 6 , 2/ 6 , 2/ 6 , 1/ 6) , g ( k) 是 Runge 2K utta 通量逼近 1g (1)g(2) g (3)g (4)f ′( u ( t n , x n ) ) ,f ′( u ( t n+ d t / 2 , x n+ d t / (2 g(1)) ) ) ,f ′( u ( t n + d t/ 2 , x n + d t/ (2g (2) ) ) ) , f ′( u ( t n + d t , x n + d t g (3) ) ) 1= -= - = - = - 类似得到单元边界点在时间 t n +1/ 2 处沿特征线方向位于时间 t n 处位置 `x j +1/ 24xn +1/ 2= x n+ d t6B k g,( k )(6)k = 1其中 , B = (5/ 24 , 4/ 24 , 4/ 24 , - 1/ 24) 1为运算简单 ,我们也可以采用中点格式来计算时间 t n +1/ 2 和 t n +1 处的网格点 `x j +1/ 2 沿特征一维 Eu ler 方程的特征有限体积格式293线落在时间 t n层的位置 1 仅介绍时间 t n +1/ 2处的处理方法 1 度在 t n 处的线性插值形式速度 u 在时间 t n +1/ 2 可以写成速u j +1/ 2 (1 - ξ) + u j - 1/ 2ξ,(7)u η = 从而可以得到单元边界格点 x j +1/ 2 在时间 t n 处的位移为u j +1/ 2Δt / 2ξΔx (8)(Δ(Δ- u ) j +1/ 2 j - 1/ 2 从而得到单元格点 `x j +1/ 2 在时间 tn +1/ 2处的位置 ( x j +1/ 2 - ξΔx ) 1 1. 1 . 3 节点值的重构为得到高阶无震荡格式 ,本文分别采用三阶和五阶中心加权本质无震荡格式 4 25来对各变量的节点值进行重构 1 这里仅简要介绍五阶本质无震荡格式4 25 1在模版 I j - 2 , I j - 1 , I j , I j +1 , I j +2 上 引 入 四 阶 最 优 多 项 式 u o ( x ) , 并 将 其 定 义 为 u o ( x ) =5 6 x i ) i - 1 1 为了得到本质无震荡格式还需要构造 3 个二阶多项式 1 类似 u o ( x ) 定3a oi - 1( x - i = 16 义 ,首先构造 3 个一般的二阶多项式为 : 在模板3I j - 2 , I j - 1 , I j 上定义 u 1 ( x) =a 1i - 1( x -i = 16 x i )i - 1, 在模板3x i ) i - 1 , 在模板 I j - 1 , I j , I j +1 上定义 u 2 ( x) =a 2i - 1( x - I j , I j +1 , I j +2 上定义i = 16 x i ) i - 1 1u 3 ( x ) =a 3i - 1 ( x -i = 1上面定义的 4 个多项式都满足下面的关系 :∫Iu ix d x = Δx ¯uj + k ,( ) ( ) 9j + k其中 ,当 i = 0 时 , k ∈ - 2 , - 1 , 0 , 1 , 2 ; 当 i = 2 时 , k ∈ - 2 , - 1 , 0 - 1 , 0 , 1 ; 当 i = 3 时 , k ∈ 0 , 1 , 2 1 这 4 个多项式的系数由文献 4 25 ; 当 i = 1 时 , k ∈可以得到 1 为得到中心加权本质无震荡格式 ,需要引入非线性权 ,首先定义中心多项式 u c ( x ) , 而上 面定义的多项式满足关系u o ( x ) =6 Cj×u j ( x ) 1(10)j由该式可以得到五阶的中心多项式为u c ( x ) = ( u o ( x ) - C 1 u 1 ( x ) - C 2 u 2 ( x ) - C 3 u 3 ( x ) ) / C o ,(11)其中 , 6 C j = 1 , j ∈ o , 1 , 2 , 3 , 同样本文也取 C 1 = C 3 = 1/ 8 , C 2 = 1/ 2 1j由上面所得到的 4 个多项式 u k ( x ) , k 质无震荡方法引入下面关系式 :∈ c , 1 , 2 , 3 可以定义本质无震荡权 1 类似加权本 2xl = 1∫x i +1/ 2IS i =6 h 2 l - 1 ( u ( l ) 2d x ,k = c , 1 , 2 , 3 1 (12)) k k i - 1/ 2整理可以得到 IS i为kIS i = a 2 Δx 2 + ( 13/ 3) a 2 Δx 4+ O (Δx 6) , (13)k k 1 k 2 其中 , k ∈ 1 , 2 , 3 1对中心多项式 u c ( x) , 可以定义 IS i 为c IS i = a 2Δx 2 + ( (13/ 3) a 2 + ( 1/ 2) a a )Δx 4 + O (Δx 6) 1 (14)c c1 c2c1 c3294郭 彦 刘 儒 勋由式 (13) 和 (14) 可以定义本质无震荡权αk C kωk , αk =(15)= 2 ,i 6(ε + IS k ) αpp ∈ c , 1 , 2 , 3其中 , k ∈ c , 1 , 2 , 3 1 对于不同的算例以及不同网格剖分取不同的小量 ε1在上 面 所 得 到 的 本 质 无 震 荡 权 的 基 础 上 可 以 定 义 单 元 c , 1 , 2 , 3 的凸组合 R j ( x )R i ( x ) = ωc u c ( x ) + ω1 u 1 ( x ) + ω2 u 2 ( x ) + ω3 u 3 ( x ) 1详细的中心本质加权无震荡构建可以参见文献 5 1 I i 上 多 项 式 u k ( x ) , k ∈ (16) ^^从而节点值 u ( x i +1/ 2 , t n + d t ) 和 u ( x i +1/ 2 , t n + d t / 2) 可从上面中心加权本质无震荡重构的 多项式 R i ( x ) 得到^ ^ u ( x i +1/ 2 , t n + d t ) = R i ( x n +1) , u ( x i +1/ 2 , t n+ d t / 2) = R i( x n +1/ 2) 1(17) 下面给出由时间 t n 层的各变量单元平均值 ¯u n 求解下一时间层的值 ¯u n +1 的算法 1i i 算法 :给定时间 t n 层处的单元平均值 ¯u n 计算下一时间层处的值i 1) 计算多项式 u k ( x ) 的系数并计算权 ωk , k ∈ c , 1 , 2 , 3 在 t n层的值 ;2) 由式 (5) 和 (6) 分别计算单元边界 x n +1 和 x n +1/ 2在时间层 t n +1 和 t n +1/ 2 沿特征线落在时间层 t n 的位置 ;^3) 利用式 (16) 计算 u ( x i +1/ 2 , t n+ βl d t ) 并计算方程 (2) 的右端项 ; 4) 由方程 (2) 计算下一时间层单元平均值 ¯u n +1 1 i 1 . 2 Euler 方程组本节将前面所提出的格式运用到一维无粘 E uler 守恒方程组 1 首先考虑如下一维无粘守 恒 E uler 方程组5 U 5 F(18)= 0 ,5 t + 5 x其中ρu ρu 2 + pρρu1U = , F =Eu ( E + p )其中 , U 为守恒变量向量形式 , F 为无粘通量 1 分别记密度为 ρ, 速度为 u , 压力为 p , 总能量为 E 且有关系式 E = p/ (γ - 1) + ρu 2/ 2 , 其中 γ = 1 . 4 1 1. 2 . 1 一维 E uler 方程的数值格式类似与构造一维标量守恒方程的格式 ,由下式得到单元平均值在下一时间层的值 U ¯n +1: i3d t ^ ^ U ¯n +1n d x 6n n = U ¯i- N l [ F ( U ( x i +1/ 2 , t + βl d t ) ) - F ( U ( x i - 1/ 2 , t + βl d t ) ) ] 1 (19)il = 1设 V = (ρ, u , p ) T 为原始变量的向量形式 ,方程 (18) 可以改写为5 V 5 V+ A p 5 x= 0 ,(20)5 t 其中ρ 0u 0 u 1/ ρ 1A p = 0 γpu一维 Eu ler 方程的特征有限体积格式295由原始变量得到的方程 (20) 是非守恒方程 1 设 c 为声速且记 c = γp/ ρ, 则 A p 的特征值为λ1 = u - c ,λ2 = u ,λ3 = u + c1 相应与这些特征值的右特征向量 ,左特征向量和特征矩阵分别 为 1/ (2 c 2) - 1/ c21/ (2 c 2) λ1λ2 0 00 λ31 - c / ρ c 210 1 010 - ρ/ ( 2 c )0 ρ/ (2 c )c/ ρ , Λ =1R p =, L p = c 2从而方程 (20) 的特征形式可以重写为5 V 5 V L p + ΛL p5 x = 0 1(21)5 t 线性化特征矩阵可以将方程 (21) 化为ρ 1 d x 2 c d u - 2 c2 d p = 0 , l 1 ( X 0):= λ1 = u - c , d t d ρ - 1 d p = 0 , l 2 ( X 0) : d x = λ2 = u , c2 d tρ 1 d x l 3 ( X 0) : d t = λ3 = u + c1 2 c d u + 2 c2 d p = 0 , 原始变量在 X 0 处的值可以由下面的式子得到 , 其中 X ( l k ) ( k = 1 , 2 , 3) 为单元边界格点分别在 时 间 t n 和 t n +1 沿特征线 l k ( k = 1 , 2 , 3) 落在时间 t n 处的位置 1 对每条特征线来说 ,这里求解节 点落在前一时间层位置的方法类似求解标量守恒率方程中求解的方法 1 从而通过求解线性方程就可以沿特征线求解 (ρ, u , p ) 为=1 ( U ( X ( l 1) ) + U ( X ( l 3) ) + 1( P ( X ( l 3) ) - P ( X ( l 1) ) ) ) ,u ( X 0) ρc2 = 1( P ( X ( l 1) ) + P ( X ( l 3) ) + ρc ( U ( X ( l 3) ) - U ( X ( l 1) ) ) ) , p ( X 0) 2= R ( X ( l 2) ) + 1( p ( X 0) - P ( X ( l 1) ) ) ,ρ( X 0) c2 其中 , R ( x ) , U ( x ) 和 P ( x ) 分别表示为ρ, u 和 p 的近似 , 这3 个节点值都可以用中心加权本 质无震荡重构得到 1 这里的重构类似与求解标量方程的重构 R i ( x ) 求解节点值 (17) 1算法 :给定时间 t n 层处的单元平均值ρ¯n ,ρ¯u n 和 ¯E n 计算下一时间层的单元平均值ρ¯n +1,ρ¯u n +1 和 i i i i i ¯E n +1 ;1) 利用时间 t n 层的单元平均值ρ¯n,ρ¯u n 和 ¯E n 计算 ¯u n 和 ¯p n ;iiiiii2) 计算在时间 t n处多项式 R k ( x ) , U k ( x ) , P k ( x ) 的系数和权 ωk , k = c ,1 ,2 ,3 , R i ( x ) = ωc R c ( x ) + ω1 R 1 ( x ) + ω2 R 2 ( x ) + ω3 R 3 ( x ) , U i ( x ) = ωc U c ( x ) + ω1 U 1 ( x ) + ω2 U 2 ( x ) + ω3 U 3 ( x ) ,P i ( x ) = ωc P c ( x ) + ω1 P 1 ( x ) + ω2 P 2 ( x ) 3) 利 用 式 ( 5) 和 ( 6) 计 算 单 元 边 界 点 在 时 间 + ω3 P 3 ( x ) ;t n +1 和 t n +1/ 2 沿 特 征 线 落 在 t n 层 的 位 置 X ( l k ) n +1 和 X ( l k ) n +1/ 2 ( k = 1 , 2 , 3) 类似式 (5) 和 (6) , 其中式子中的特征值分别替换为λ1 ,λ2 和λ3 ;^4) 利用方程 (16) 和 (19) 的右端项可以得到单元节点值 U ( x i +1/ 2 , t n+ βl d t ) ; 5) 最后利用方程 (19) 计算下一时间层的各变量单元平均值 U ¯n +1 1 i296郭 彦 刘 儒 勋2 数 值 结 果本节主要将本文所提出的格式运用到一些一维的算例 1 首先将该格式运用到经典的线性 ,非线性 Burgers 方程和一维欧拉方程的密度扰动问题检验该格式的精度和收敛性 1例 1 求解如下线性标量方程 :u t + u x = 0 ,(22)分别给出两个初值条件 u ( x ,0) = sin (πx ) 和 u ( x ,0) = sin 4 (πx ) ,考虑以 2 为周期的边界条件 1 两种初值都计算到时间 t = 10 1 第一个初始条件下的精度和收敛性可以由表 1 得到 1 从表 1 可以看到该格式几乎可以得到理论上的三阶和五阶精度且收敛性也是比较好的 1 对于初始 = sin 4 (πx ) 可以得到类似的结论 1u t + u x = 0 , u ( x , 0) = sin (πx) 的精度和收敛性 ( t = 10)条件 u ( x , 0) 表 1e L 误差1O L 阶1e L 误差∞O L 阶∞方法n10 2 . 31 E - 01 3 . 13 E - 01206 . 58 E - 02 1 . 81 1 . 47 E - 01 1 . 10 4080 7 . 52 E - 034 . 00 E - 04 3 . 134 . 23 2 . 00 E - 021 . 25 E - 032 . 904 . 00 3 阶160 3 . 97 E - 05 3 . 33 8 . 71 E - 05 3 . 80 320 4 . 98 E - 06 3 . 008 . 61 E - 06 3 . 341020 3 . 72 E - 021 . 77 E - 03 6 . 22 E - 023 . 22 E - 034 . 39 4 . 27 405 . 71 E - 05 4 . 96 1 . 11 E - 04 4 . 86 5 阶80160 1 . 81 E - 065 . 70 E - 08 4 . 984 . 99 3 . 54 E - 061 . 09 E - 07 4 . 975 . 03 3201 . 79 E - 095 . 003 . 37 E - 095 . 02u t + ( u 2 / 2) x = 0 , u ( x , 0) = 0. 5 + sin (πx ) 的精度和收敛性 ( t = 0. 5/π) 表 2e L 误差1O L 阶1e L 误差∞O L 阶∞方法n102040 5 . 13 E - 021 . 76 E - 021 . 69 E - 03 1 . 62 E - 019 . 37 E - 021 . 08 E - 02 1 . 553 . 38 0 . 793 . 12 3 阶80160320 1 . 12 E - 041 . 06 E - 051 . 29 E - 06 3 . 913 . 403 . 048 . 18 E - 041 . 04 E - 041 . 32 E - 053 . 722 . 972 . 9810 1 . 47 E - 02 7 . 73 E - 022040801 . 32 E - 031 . 03 E - 044 . 71 E - 06 3 . 483 . 694 . 45 1 . 46 E - 021 . 43 E - 036 . 61 E - 05 2 . 413 . 354 . 445 阶1601 . 73 E - 07 4 . 762 . 12 E - 06 4 . 96 3206 . 18 E - 094 . 804 . 97 E - 085 . 42例 2 数值求解如下非线性方程 :一维 Eu ler 方程的特征有限体积格式297u 2= 0 , (23) = 0. 5/π时解仍为光滑u t + 2x其中 ,初值为 u ( x , 0) = 0. 5 + sin (πx ) , 以 2 为周期 1 当数值求解到 t 解 1 表 2 给出数值解的 L 1 和 L ∞ 误差估计 ,同时给出收敛性情况 1 从该表格中同样可以得到 理论中的三阶和五阶精度并且可以看出格式的收敛性 1 随着时间的增加数值解出现激波现 象 ,该格式也能较好地捕捉到激波 1例 3 E uler 方程的密度扰动问题11本例考虑 E uler 方程组 (18) 1 为了得到所提格式对 E uler 方程的收敛性以及精度分析 ,我们考虑初始条件为 ρ( x , 0) = 1 + 0 . 2sin (πx ) , u ( x , 0) = 1 , p ( x , 0) = 1 , 同样取以 2 为周期的 边界条件 1 E uler 方程在此初值条件下的真解为 :ρ( x , t ) = 1 + 0. 2sin (π( x - t ) ) , u ( x , t ) = 1 , p ( x , t ) = 1 1 由本文格式所得到的数值解在时间 t = 2 时密度的 L 1 和 L ∞误差从表 3 可以看到 1 由此表可以看出三阶和五阶格式都能达到理论阶 ,数值上得到格式的收敛性 1 另外通过 所得到数值信息可以得出 :本格式与本质加权无震荡和中心加权本质无震荡格式所得到的精 度阶和收敛性具有类似的效果 1表 3初始密度 ρ0 = 1 + 0. 2sin(πx ) , 初值速度压力 u 0 = p 0 = 1 的精度和收敛性e L 误差1O L 阶1e L 误差∞O L 阶∞方法n10 1. 185 E - 02 3 . 41 E - 02201 . 39 E - 03 3 . 09 4 . 23 E - 03 3 . 01 4080160 1 . 79 E - 042 . 40 E - 053 . 14 E - 06 2 . 962 . 902 . 94 4 . 02 E - 044 . 33 E - 055 . 21 E - 06 3 . 403 . 213 . 06 3 阶320 3 . 98 E - 07 2 . 986 . 41 E - 07 3 . 0210201 . 99 E - 031 . 25 E - 04 6 . 53 E - 033 . 36 E - 04 3 . 99 4 . 28 40 4 . 42 E - 06 4 . 82 1 . 10 E - 05 4 . 93 5 阶801603201 . 47 E - 074 . 75 E - 091 . 45 E - 104 . 914 . 955 . 033 . 22 E - 078 . 93 E - 092 . 55 E - 105 . 095 . 175 . 13下面 3 个数值算例用来考察该格式对激波的捕获能力 1例 4 Lax 问题12考虑 Lax 激波管问题 ,这种问题是由方程 (18) 在下面初值条件下所得到的物理现象 :(0. 445 , 0 . 698 , 3 . 528) , (0. 5 , 0 , 0 . 571) ,- 5 ≤ x ≤0 ,0 < x ≤5 1(ρ, u , p ) =本文所提格式所得到的时间 t = 1. 3 时数值解与 C WE NO 格式得到的数值解 5和解析解 的比较可以从图 1 中观察到 1 由图 1 可以看到 ,在网格数为 n = 200 时 ,该格式所算的数值结 果有一些震荡 ,但是 ,当我们继续加密网格时 ,就能明显地看到震荡在逐渐变弱 1 在引入本质 无震荡权抑制震荡的同时 ,该格式还是能较好地捕获到激波现象 1 将此数值结果与中心加权 本质无震荡的格式 C WE NO 和 UW E NO 5 得到的数值结果相比较 ,可看到在抑制震荡方面本格 式要好一些 ,能得到较为准确的数值解 1298郭 彦 刘 儒 勋(a ) 密度 ( b ) 速度 (c ) 压力( d ) 密度 (e ) 速度 (f ) 压力图 1 Lax 问题在 t = 1. 3 时的密度 、速度 、压力的数值解和解析解的比较 (网格数为 n = 200 1上排图为本文所提格式 ( G CF L = 0. 5) ,下排图为 C WE NO 格式 ( G CF L = 0. 5 )例 5 S od 问题13同样考虑 E uler 方程 (18) 1 该方程的解由 3 部分组成 1 我们考虑下面初值条件 :(1 , 0 , 1) ,( 0. 125 , 0 , 0 . 1) ,- 0. 5 ≤ x ≤0 ,0 < x ≤0. 5 1(ρ, u , p ) =在时间 t = 0. 2 处的数值解与真解的比较可以从图 2 观察到 1 通过与真解比较可以得出 该格式 ,能较好捕获激波位置 ,并且在间断处与真解也能较好的吻合 1 图 2 给出的是 n = 200 时的数值结果 ,可以看出此时还是有一点震荡 ,但是随着网格的进一步加密震荡会减弱 1 这些 数值结果与其它的高阶方法特别是本质无震荡型方法所得到的结果一致 1(a ) 密度 ( b ) 速度 (c ) 压力图 2 So d 问题在 t = 0. 2 时的密度 、速度 、压力 (网格数为 n = 200 )例 6 激波相互扰动问题2本例依然考虑方程 (18) 在下面初值下的数值模拟结果 :一维Eu ler 方程的特征有限体积格式299(3. 857 148 , 2. 629 369 ,10 .333 333) ,(1 + 0. 2sin (5 x) , 0, 1) ,- 5 ≤x ≤- 4 ,- 4 < x ≤5 1(ρ, u , p) =图3 给出在时间t= 1. 8 时数值结果与由五阶中心加权本质无震荡格式,在网格数为2 000的情形下所得到的“真解”1 同别的一些高阶格式类似,本格式也不能较好地模拟激波下游正确的流动特征1 但是,随着网格地加密数值结果也是会更接近真实解1(a) 密度( b) 速度(c) 压力图3 激波交互扰动问题在t = 1. 8 时的密度、速度、压力(网格数为n = 400 )3 总结本文运用特征线理论,结合中心加权本质无震荡插值思想,对一维的双曲守恒律方程提出一种新的有限体积格式1 通过算例可以发现所提格式能较准确的给出数值解,通过对几个经典的例子的数值求解,可以得出所提出的三阶和五阶格式在光滑区域能达到理论阶,在间断处仍然保持本质无震荡的性质1 该方法所得到的数值解既保持WE NO 格式的本质无震荡性的优点,同时,相对与已有的一些的WE NO 和C W E NO 格式,该格式又体现出特征线方法的分辨率高的优点1[ 参考文献]S h u C W , Os h e r S .Efficie n t i m p le m e n t a ti o n of es s e n tially no n2os c illat o ry s h oc k cap t u ri n g s c h e m esJ . J Comp u t Ph ys , 1988 , 77 (2) :4392471.S h u C W , O s he r S . E fficie n t i m p le m e n t a ti o n o f es s e n tially no n2os c illat o ry s h oc k cap t u ri n g s c h e m es ⅡJ . J C o mp ut Phys , 1989 , 83 (1) :32278.[1 ][ 2 ][ 3 ] J Co mp ut Phys , 1996 , 126 J ia n g G , S h u C W. E fficie n t i m p le m e n t a ti o n of w e ight e d ENO s c h e m es J .(1) :2022228.[ 4 ] Levy D , Pup o G , Rus s o G. Comp a ct ce n t r al W E NO s c h e m es f or m ulti d i m e n si o nal c o ns e r vati o n laws J . SIAM J Sci C o mp ut , 2000 , 22 (2) :6562672.Cap deville G. A ce n t r al W E NO s c h e m e f or s o lvi n g hyp e r b olic c o ns e r vati o n laws on non2unif o r mmes hes J . J C o mp ut Phys , 2008 , 227 (5) :297723014.陈荣三. 大密度和大压力比可压缩的数值计算J . 应用数学和力学, 2008 , 29 (5) :6092617.涂国华, 袁湘江, 陆利蓬.激波捕捉差分方法研究J . 应用数学和力学, 2007 , 28 (4) :4332440.HU J un , G UO S h a o2ga n g . S o l u ti o n t o Eule r equati o ns by high2res ol u ti o n u p w i n d c o mp a ct s c h e m e [ 5 ][ 6 ][ 7 ][ 8]I n te r n a t J N u m e r Met h Fl u i ds , 2008 , 56 (11) :213922150.bas e d o n fl u x sp litti n g J .[ 9 ] Xia o F , P e n g X. A c o nve x it y p res e r vi n g s c h e m e f or c o ns e r vative a d vecti o n t r a n sp ort J . J C o mp ut Ph ys , 2004 , 198 (2) :3892402.300郭 彦 刘 儒 勋Ii S , Xiao F . CI P / m ulti 2mome n t fi n it e vol u me met h o d f or Eule r equati o ns :A s e m i 2Lagra n gia n c h a r ac 2 [ 10 ] J Co mp ut Phys , 2007 , 222 (2) :8492871.t e r is t ic f or m ulati o n J .[ 11 ]Qi u J , S h u C W. He r m it e W E NO s c h e m es a n d t h eir app licati o n as li m it e r s f or Runge 2Kutt a dis c o nti n 2 uo us Gale r k i n met h o d : o ne di m e n si o nal cas e J . J Co mp ut Phys , 2004 , 193 (1) :1152135.La x P D . Wea k s o l u ti o ns o f nonli n ea r hyp e r b olic equati o ns a n d t h eir n ume r ical c o mp ut a ti o n J . Com m u n P u r e Appl Ma t h , 1954 , 7 (1) :1982232.S o d G. A s urvey of s e ve r al fi n it e diff e r e n ce met h o ds f or s y s t e m s of non 2li n ea r c o ns e r vati o n laws J . J Comp u t Ph ys , 1978 , 27 (1) :1231.[ 12 ] [ 13 ]C h a r a c t e r i s t i c 2B a s e d Fi ni t e Vol u m eS c h e m e f o r 1D E u l e r E q u a t i o n sGUO Y a n , LIU Ru 2x u n( Dep a r t m e n t of Ma t h e m a t i c s , U n i v e r s i t y o f Sci e n ce a n d Tech n o logy of Chi n a ,Hef e i 230026 , P. R. Chi n a )A b s t r a c t : A high 2or d e r fi n it e 2vol u me s c h e m e w a s p res e n t e d f or t h e o ne 2di m e n si o nal s c ala r a n d i n vis 2 ci d Eule r c o ns e r vati o n laws . The Si m p s o n ’s qua d rat u re r u le w a s us e d t o ac h ieve high 2or d e r acc u racy i n ti m e . To g et t h e p o i n t val u e of t h e Si m p s o n ’s qua d rat u re , t h e c h a r act e r is t ic t h e o ry w a s us e d t o o b 2 t a i n t h e p ositi o ns o f t he gri d p oi n ts at eac h s ub 2ti me s t ages al o ng t he c h a ract e r is t ic c u rves , a n d t h e t h ir d 2or de r a n d fif t h 2or de r ce n t ral w eight e d es s e ntially non 2os c illat o ry ( C W E NO ) rec o ns t r u cti o n w a s a dop t e d t o es t i mat e t he cell p o i n t val ues . Seve ral s t a n da r d one 2di me n si o nal e xa m p les w e r e u s e d t o ve r if y high 2or d e r acc u racy , c o nve r ge n ce a n d cap a b ilit y o f cap t u ri n g s h oc k .Ke y w o r ds : hyp e r b olic equati o n ; fi n it e vol u me met h o d ; c h a r act e r is t ic t h e o ry ; W E NO rec o ns t r u cti o n ;Runge 2Kutt a met h o d。
二维抛物方程的有限差分法
二维抛物方程的有限差分法二维抛物方程的有限差分法摘要二维抛物方程是一类有广泛应用的偏微分方程，由于大部分抛物方程都难以求得解析解，故考虑采用数值方法求解。
有限差分法是最简单又极为重要的解微分方程的数值方法。
本文介绍了二维抛物方程的有限差分法。
首先，简单介绍了抛物方程的应用背景，解抛物方程的常见数值方法，有限差分法的产生背景和发展应用。
讨论了抛物方程的有限差分法建立的基础，并介绍了有限差分方法的收敛性和稳定性。
其次，介绍了几种常用的差分格式，有古典显式格式、古典隐式格式、Crank-Nicolson隐式格式、Douglas差分格式、加权六点隐式格式、交替方向隐式格式等，重点介绍了古典显式格式和交替方向隐式格式。
进行了格式的推导，分析了格式的收敛性、稳定性。
并以热传导方程为数值算例，运用差分方法求解。
通过数值算例，得出古典显式格式计算起来较简单，但稳定性条件较苛刻；而交替方向隐式格式无条件稳定。
关键词：二维抛物方程；有限差分法；古典显式格式；交替方向隐式格式FINITE DIFFERENCE METHOD FORTWO-DIMENSIONAL PARABOLICEQUATIONAbstractTwo-dimensional parabolic equation is a widely used class of partial differential equations. Because this kind of equation is so complex, we consider numerical methods instead of obtaining analytical solutions. finite difference method is the most simple and extremely important numerical methods for differential equations. The paper introduces the finite difference method for two-dimensional parabolic equation.Firstly, this paper introduces the background and common numerical methods for Parabolic Equation, Background and development of applications. Discusses the basement for the establishment of the finite difference method for parabolic equation And describes the convergence and stability for finite difference method.Secondly, Introduces some of the more common simple differential format,for example, the classical explicit scheme, the classical implicit scheme, Crank-Nicolson implicit scheme, Douglas difference scheme, weighted six implicit scheme and the alternating direction implicit format. The paper focuses on the classical explicit scheme and the alternating direction implicit format. The paper takes discusses the derivation convergence,and stability of the format . The paper takes And the heat conduction equation for the numerical example, using the differential method to solve. Through numerical examples, the classical explicit scheme is relatively simple for calculation, with more stringent stability conditions; and alternating direction implicit scheme is unconditionally stable.Keywords:Two-dimensional Parabolic Equation; Finite-Difference Method; Eclassical Explicit Scheme; Alternating Direction Implicit Scheme1绪论1.1课题背景抛物方程是一类特殊的偏微分方程，二维抛物方程的一般形式为u Lu t ∂=∂ (1-1)其中1212((,,))((,,))(,,)(,,)(,,)u u u u u u L a x y t a x y t b x y t b x y t C x y t x x y y x y∂∂∂∂∂∂=++++∂∂∂∂∂∂ 120,0,0a a C >>≥。
一类线性双曲型偏微分方程的有限差分格式求解
一类线性双曲型偏微分方程的有限差分格式求解刘付军;龚东;王高放【摘要】针对一类二阶双曲型偏微分方程,利用有限差分法建立了显式和隐式两种差分格式.对两种差分格式进行加权平均,得到了一种新的加权平均格式,给出了新加权平均差分格式解的存在性、收敛性和稳定性分析,最后给出了数值算例验证.【期刊名称】《河南工程学院学报（自然科学版）》【年(卷),期】2014(026)003【总页数】4页(P73-76)【关键词】显式差分格式;隐式差分格式;加权平均;收敛性【作者】刘付军;龚东;王高放【作者单位】河南工程学院理学院,河南郑州451191;河南工程学院理学院,河南郑州451191;河南工程学院理学院,河南郑州451191【正文语种】中文【中图分类】O241.84有限差分法是求解微分方程的一种有效的数值解法.对于线性双曲型偏微分方程,已建立了一些典型差分格式,如显式、隐式和紧差分格式等[1-4].文献[5]利用差分格式对二阶常微分方程进行了求解.文献[6]利用一种加权平均格式对一阶双曲型偏微分方程进行了数值分析.双曲型偏微分方程是描述振动或波动现象的一类重要的偏微分方程,在实际问题中有着广泛的应用[7].针对二阶常系数线性双曲型偏微分方程建立了两种差分格式并对这两种差分格式进行加权平均,求出加权平均值,建立了一个新的差分格式,验证了该差分格式解的存在性、收敛性和稳定性.考虑二阶双曲型问题区域记作Ω={(x，t)|0≤x≤1，0<t≤T}.对区域Ω进行剖分.记xi=ih, 0≤i≤m；tk=kτ，0≤k≤n.其中,h=1/m，τ=T/n.记Ωh={xi|0≤i≤m}, Ωτ={tk|0≤k≤n}, Ωhτ=Ωh×Ωτ.称在t=tk上的节点{(xi,tk)|0≤i≤m}为第k层节点.定义Ωhτ上的网格函数用来表示,其中).此外,记n.1.1 显式差分格式在节点(xi,tk)上考虑定解问题(1),有将u(xi-1, tk)和u(xi+1, tk)分别在节点(xi, tk)处进行Taylor展开,得到其中,xi-1<ζik<xi+1.将u(xi, tk-1)和u(xi, tk+1)分别在节点(xi, tk)处进行Taylor 展开,得到其中,tk-1<ηik<tk+1.将式(3)和式(4)代入式(2)得到针对问题(1),建立如下近似差分格式：该格式称为3层5点显式差分格式.1.2 隐式差分格式针对给出如下的差商格式:将u(xi,tk-1)和u(xi,tk+1)分别在节点(xi,tk)处进行Taylor展开,得将式(7)和式(8)代入到式(2)中,可得如下近似差分格式:该格式称为3层5点的隐式差分格式.1.3 加权平均差分格式对显格式(6)和隐格式(9)进行加权平均,加权平均数用θ表示,得到将在点(xi，tk)利用Taylor展开,得到其截断误差记为将式(11)至式(14)代入式(10),可得为了便于处理,不考虑二阶项,由式(1)可得=a2,分别对x和 t求四阶偏导,得=a2，=a2，于是有=a4,将其代入式(15),可得令--+=0，则θ=, 令r=,则θ=1+, 1-θ=-, 代入式(10)可得如下差分格式:该格式称为一种新的加权平均差分格式.首先,得到该加权平均差分格式解的存在性.定理1 加权平均差分格式(17)的解是存在唯一的.证明设记r=ατ/h, 给定u0和u1,则将格式(17)进行变形,得到利用上式当中的递推关系,计算可以得到因而差分格式(17)是显式的,对于任意的r 都是唯一可解的.下面给出该加权平均差分格式解的收敛性和稳定性.定理2 加权平均差分格式(17)的解是收敛的.证明设{u(x,t)|0≤x≤1,0≤t≤T}是定解问题(1)的解,是差分格式(17)的解.记则由式(16)可知当h,τ→0时,误差e(h3+τ3)→0, 差分格式(17)的解收敛到定解问题的精确解,则差分格式(17)与相应的微分方程(1)相容,故该差分格式是收敛的且差分格式具有三阶精度.定理3 加权平均差分格式(17)的解是稳定的.证明 Lax定理指出，对于一个适当提出的线性微分方程初值问题以及它的一个满足相容性条件的差分逼近,收敛性的充分必要条件是稳定性.根据前面已经证明了的差分格式(17)的相容性和收敛性,由Lax定理即可得差分格式(17)的稳定性.应用新的加权平均差分格式(17)计算如下定解问题：该定解问题的精确解为u(x,t)=ex+t.下面给出了当步长h=1/10,τ=1/20(即步长比r=1/2)时,计算得到的部分数值结果,如表1所示.由表1可以看出，精确解与数值解的误差控制在量级范围之内,且当h，τ→0时,所得的结果越接近精确解.【相关文献】[1] 李立康.微分方程数值解法[M].上海:复旦大学出版社,2003.[2] 孙志忠.偏微分方程数值解法[M].北京:科学出版社,2012.[3] 胡建伟,汤怀民.微分方程数值方法[M].北京:科学出版社,2007.[4] 李胜坤,冯民富,李珊.Benjamin-Bona-Mahony方程的有限差分近似解[J].四川师范大学学报:自然科学版,2003,26(4):363-365.[5] 张守贵.用差分法求解二阶常微分方程初值问题[J].重庆理工大学学报:自然科学版,2012,26(8):110-112.[6] 杨韧.求解一阶线性双曲型偏微分方程组的一个差分格式[J].四川师范大学学报:自然科学版,2009,32(5):614-617.[7] 谷超豪,李大潜,陈恕行,等.数学物理方程[M].北京: 高等教育出版社，2002.。
偏微分方程的有限差分方法
二阶线性偏微分方程的一般形式为：
A 2 u B 2 u C 2 u D u E u F G u 0 x 1 2 x 1 x 2 x 2 2 x 1 x 2
对于变量 x1 和 x 2 给定的值 xˆ1 和 xˆ 2 若 4 A (x ˆ 1 ,x ˆ 2 ) C (x ˆ 1 ,x ˆ 2 ) B 2 (x ˆ 1 ,x ˆ 2 )
这里，[ u ] ij 表示 u(xi, yj )。上两式分别简记为
x p u x ijh 1 1 2x(pijx[u]i)jO (h1 2)
yp u yijh 12 2y(pij y[u]ij)O (h2 2)
则 L u x p u x y p u y q u f (x ,y ) 在 (i, j) 点被表示为
余弦是 (co,scos)。
由
u nij
u xijc
os u yijc
os
用单侧差商逼近 x方向和 y方向的导数，然后列
出边界网点上的差分方程。
（2）邻近边界的网格点 (xi , yj ) 不在上 可以采用直接转移法近似处理，即将边界
条件用于邻近边界的网格点，然后再在该点列 出差分方程。
2 用积分插值法构造差分格式 3 差分格式的稳定性和收敛性 4 差分方程求解的一些方法
— 数值积分 有限元法
— 函数插值
不同的数值微分和数值积分方法、不同的函数插值方 法，就产生了不同的有限差分法与不同的有限元法。
其它数学基础：
数理方程、数值代数、最优化理论与方法等
偏微分方程的有限差分方法
基本思想：使用离散的、只含有限个未知 数的差分方程去近似代替连续变量的微分方程 及边值条件，并将相应的差分方程解作为(初)边 值问题的近似解。
二维波动方程的一种高精度紧致差分方法
（一 ） —２） ＝ ） Ａ （ ， （ ； Ｔ２ ＋ ×
ｅ
，
Ｊ
＋Ｒ
０≤ ， ≤ 一１， Ｏ≤ ｎ≤ Ⅳ 一 １
用 ｒ表示时 间步 长 ， ｈ表 示空 间步 长 ， １给 出 了 当 ７＝ 表 －
（） ８
ｅ ＝ｏ
ｈ ［ （ Ｌ 曰
Ａ
（
３ 数值 验 证
对于式 （ ）一（ ） 令 （ ｙ ￣ ， ｉ（ ）ｉ（ｒ） （ １ ４， ， ）＝ ２ ｒｎ ｓ Ｔ ， ， ｒ ｓ ｎ ｙ
Ｙ ）＝０ｇ ， ，） ， 问题 的精确解 为 ／ ， ，）＝ ｉ（ 订 ） ，（ Ｙ ｔ ＝０ 则 ２ Ｙ ｔ ｓ ｔ （ ｎ ｓ （Ｔ）ｉ（ｒ） ｉ ＇ ｓ 叮 。数值 实验计算是用 Ｆ ￣ａ ７语言进行编程 ｎＩ ｎ ｙ Ｘ ｏ ｒｎ７
ｈ） 。 阶精度的数值 解。孙 志 忠 提 出 了求 解二 维波 动方 程 的 高精度交替方 向隐式 方法 ， 并且是无条件稳定 的。有关这方 面 最新 的一 些 工作 可参 见 文献 ［ ６—８ 。本文 在 此工 作基 础 之 ］
上， 利用 Ｒｃ ａ ｓｎ外推法进 一步 提 高计 算精 度 ， ｉ ｒｏ ｈｄ 最终 可得 到
ｏ≤ √≤ 一１
４ ｔ ｈ ，＝１时刻 ， 本文格 式在不 同网格步 长下误 差 的 、 、
范数 ， 以及与四阶 Ａ Ｉ Ｄ 格式 计 算结 果 的 比较 。Ｌ 范 数定 义 ２
厂］ｉ 面＝广——一 『
ｄ］ ＝ ｏ 。 ｅ
＼ ｅ ＝０
。 ≤ Ｍ
一
１
层的。即每一次时间推进都需要知道前 两个 时间步 的值 ， ０ 第
求解该 问题 精度为 Ｏ（ ＋ 。 的数值解 。 ｈ）