计算电磁学2有限差分法

差分格式
1 ⎛ ∂ 2ϕ ⎞ ⎛ ∂ϕ ⎞ α (ϕ1 − ϕ0 ) + β (ϕ3 − ϕ0 ) = ⎜ ⎟ (α h1 − β h3 ) + ⎜ 2 ⎟ (α h12 + β h32 ) + 2! ⎝ ∂x ⎠0 ⎝ ∂x ⎠0
' n '' n ''' n
如果该关系式能够一直维持到第m阶仍能成立， 但m＋1阶不再成立，就称为m阶龙格-库塔法
γ j = f n , j = 1, 2, , m
j −1
cem@
龙格-库塔法
⎧ω1 + ω 2＝1 ⎪ 1 二阶龙格-库塔法的系数方程组 ⎨ ⎪α 2ω 2＝ β 21ω 2＝ 2 ⎩ 这是一个不定方程组，也就是方程的变量的数目 多于方程个数，因而方程的解有多个。上式如果 取α2=β21=1，ω1=ω2=1/2
⎧ K1 = hf ( xn , yn ) ⎪ ⎪ K 2 = hf ( xn + α 2 h, yn + β 21 K1 ) ⎪ ⎨ K 3 = hf ( xn + α 3 h, yn + β 31 K1 + β 32 K 2 ) ⎪ ⎪ ⎪ K r = hf ( xn + α r h, yn + β r1 K1 + β r 2 K 2 + ... + β r ,r −1 K r −1 ) ⎩ cem@
龙格-库塔法
然后，再用ω常数对它们进行线性组合
K＝ω1 K1＋ω2 K 2＋...+ωr K r
并且作为yn+1到yn的增量，得到n+1步长处的结果
yn +1 = yn + K＝yn + ω1 K1＋ω2 K 2＋...+ωr K r
这些与计算函数无关的常数如何去确定？这也是 龙格-库塔法的关键，思路是：即要计算函数的增 量，就要想到泰勒展开方法，上式中每个K值就 是在处附近的值。
⎧ yn +1 = yn + K1 / 2 + K 2 / 2 ⎪ ⎨ K1＝hf ( xn , yn ) ⎪ K ＝hf ( x + h, y + K ) n n 1 ⎩ 2
这就是梯形近似的预估计-校正法。其他阶的龙 格-库塔法可以类推
cem@
有限差分格式
二维Possion方程差分格式
cem@
1、差分与差商
用差分代替微分，是有限差分法的基本出发点。 这一点由微分原理保证的，当自变量的差分趋于 零时，差分变成微分
Δf ( x ) = f ( x + h ) − f ( x )
df Δf ( x ) = lim f ( x) = dx Δx ⎯⎯ 0 Δx →
2、欧拉近似、梯形法则和龙格-库塔 法
除了用最简洁差分去代替微分外，还可以从微分方程解 为积分形式出发，利用数值积分计算的手段处理差分格 式，通过直观图形法去分析。
⎧ dy ⎪ = f ( x, y ) ⎨ dx ⎪y ⎩ x = x 0 = y0
y ( x) = y0 + ∫ f (t , y (t )dt
⎜ ⎟ Δx ⎝ dx x +Δx dx x ⎠ dx 1 ⎡ f ( x + h) − f ( x ) f ( x ) − f ( x − h) ⎤ ≈ ⎢ − ⎥ h⎣ h h ⎦ f ( x + h) − 2 f ( x ) + f ( x − h) = h2
2
2
''
后向 差分
cem@
⎧ yn +1 ≈ yn + hf ( xn , yn ) ⎨ ⎩ xn = x0 + nh
cem@
梯形近似
改善数值分析的精度关键是如何更好地近似曲边梯形的面 积。用函数曲线的端点连线代替曲线，而不是用矩形的平 行线
hn yn+1 ≈ yn + ( f (xn, yn ) + f (xn+1, yn+1)) 2
cem@
龙格-库塔法
采用预估计-校正法的梯形近似法，只计算一次就 能得到的截断误差为O（h3）结果，但当期望一 步到达R阶精度，就发展了龙格-库塔法 选取与函数f无关又与步长数n无关的常数 α 2 , β 21 , α 3 , β 31 , β 32 ,..., α r , β r1 , β r 2 ,..., β r ,r −1及ω1 , ω2 ,...ωr 用αβ常数依次计算一系列的值
'
两式相减
2 3 ''' f ( x0 + h) − f ( x0 − h) = 2hf ( x0 ) + h f ( x0 ) + 3!
'
cem@
差分与差商
前向、后向差分截断于 h f ( x0 ) / 2! ，具有h的一阶 精度，而中心差分法截断于 2h3 f ''' ( x0 ) / 3!，具有h的 二阶精度，中心差分的精度比较高。 前向 d2 f 1 ⎛ df df ⎞ 差分 = −
2 2
h αh + βh = 0 ⇒ α = − β h
2 1 2 3
2 3 2 1
求出二阶精度精度为一阶偏导数差分格式
2 2 ⎛ ∂ϕ ⎞ α (ϕ1 − ϕ0 ) + β (ϕ3 − ϕ0 ) h3 (ϕ1 − ϕ0 ) − h1 (ϕ3 − ϕ0 ) = ⎜ ⎟ ≈ α h1 − β h3 h1h3 (h1 + h3 ) ⎝ ∂x ⎠0 ⎛ ∂ϕ ⎞ h42 (ϕ 2 − ϕ0 ) − h22 (ϕ 4 − ϕ0 ) ⎜ ⎟ ≈ cem@ h2 h4 (h2 + h4 ) ⎝ ∂y ⎠0
∂ 2ϕ ∂ 2ϕ ∇ 2ϕ ( x, y ) = 2 + 2 = F ( x, y ) ∂ x ∂ y
有限差分法的网格划分，通常采用完全有规律的 分布方式，这样可使每个离散点上得到相同形式 的差分方程，有效的提高解题速度。对能填满平 面域的三种规则网格（正方形，正三角形和正六 边形）的划分方式，经常采用的是正方形网格划 分，
'
Δf ( x ) f ( x + h ) − f ( x ) f ( x) ≈ = Δx h
'
cem@
差分与差商
前向差分
df Δf ( x) f ( x + h) − f ( x) ≈ = Δx dx h
后向差分 df Δf ( x) f ( x) − f ( x − h) ≈ = Δx dx h 中心差分
cem@
差分格式
h2
3
h3
2 0
h1
ϕ2
j+1 1
一阶偏导数差分格式
⎧⎛ ⎪⎜ ⎪⎝ ⎨ ⎪⎛ ⎜ ⎪⎝ ⎩
ϕ − ϕ0 ∂ϕ ⎞ ≈ 1 + O (h 2 ) ⎟ h1 ∂x ⎠0 ϕ − ϕ3 ∂ϕ ⎞ ≈ 0 + O (h 2 ) ⎟ ∂x ⎠0 h3
h4
ϕ3
ϕ0
ϕ4
差分与差商
对偏导数，可仿照上述方法，将表示为：
∂u u ( x + h, y, z ) − u ( x, y, z ) ≈ h ∂x
∂ u u ( x + h, y, z ) − 2u ( x, y, z ) + u ( x − h, y, z ) ≈ 2 h2 ∂x
2
cem@
df Δf ( x) f ( x) − f ( x − h) ≈ = Δx dx h
cem@
差分与差商
通过泰勒公式分析上面差分精度，在点上的一阶 导数的逼近度可由泰勒公式展开
1 2 '' f ( x0 + h) = f ( x0 ) + hf ( x0 ) + h f ( x0 ) + 2! 1 2 '' ' f ( x0 − h) = f ( x0 ) − hf ( x0 ) + h f ( x0 ) − 2!
hn yn+1 ≈ yn + ( f (xn, yn ) + f (xn+1, yn+1)) 2 ⎧ yn +1 = yn + hf ( xn , yn ) ⎪ ⇒⎨ hn ⎪ yn +1＝yn + 2 ( f ( xn , yn ) + f ( xn +1 , yn +1 ) ) ⎩
cem@
f(x)
f(xn+1) f(xn)
上式为梯形近似法， 又称改进的欧拉折线法， 改进后的误差为，具有二阶精度， 比欧拉折线法小一个数量级
xn
xn+1
x
cem@
梯形近似
可以看出方程两边都包含n＋1步的待求未知量， 一般情况下无法直接求解，为此通常采用“预估 计-校正法”：先用欧拉折线法计算出0到n＋1步的 未知值的预估计值，然后，再用梯形近似法计 算。该方法的截断误差仍然为O（h3）。
Hale Waihona Puke 龙格-库塔法一般梯形近似法不能一步计算得到所要求的精 度，需要采用多次迭代才能够得到积分精度要求
[0] ⎧ yn +1 = yn + hf ( xn , yn ) ⎪ ⎨ [ s +1] hn [ yn +1 ＝yn + ( f ( xn , yn ) + f ( xn +1 , yns+]1 ) ) , s = 0,1, 2,... ⎪ ⎩ 2 [ [ max { yns++1] − yns+]1 } < ε 1 当迭代到一定的次数时，即满足精度 n∈[0, N ] 时，就认为获得所要求的精度积分。