差分法によるラプラス方程式の数値計算（pdf

差分法によるラプラス方程式の解法
電気情報工学科 ∗ 5 年通年コンピュータシミュレーション
2015 年 10 月 25 日 (金)
後期第 4 回目（通期第 18 回目）
概要
偏微分方程式の数値計算を考える．我々が生活する世界は時間も含めると 4 次元
空間であるため，独立変数が 1 つの常微分方程式ではなく，偏微分方程式が支配する
現象がほとんどだからである．
1 はじめに
偏微分方程式の解法は，過去に応用解析の講義などで学習している．そこでは，以下の
ような手順で偏微分方程式の解を求めたと思う．
1. 解を変数分離して，偏微分方程式を連立の常微分方程式に直す．
2. 境界条件を満たすように常微分方程式を解く．
3. 各々の常微分方程式の解の積がもとの偏微分方程式の解となる．
このようにして得られる解析解は，厳密で誤差はゼロである．でも，厳密な解析解が得ら
れるのは境界が単純な場合に限られる．通常の物理問題では，複雑な境界のもと偏微分方
程式の解の値が必要となる．このようなときに，数値計算が便利になってくる．
ここでは，ラプラス方程式を差分法というテクニックで数値計算する方法をみてみる．
実際には，図 1 の静電磁場や，図 2 の熱の問題に，このラプラス方程式は表れる．
• 図 1 は，紙面と垂直方向に無限に長い正方形の金属筒に，電線が 2 本通っている．正
方形の筒は 0V にアースされており，2 本の電線はそれぞれ，30V と-20V である．こ
の状態で，筒内部のポテンシャル (電圧) の分布を求めなさいという問題である．
• 図 2 は，紙面方向に非常に長い豆腐があり，その周りは 0 ℃の水で満たされている．
そして，その豆腐に 30 ℃と-20 ℃の金属棒が突き刺さっている状況である．豆腐内
部の温度分布を求めなさいというのが問題である．
∗
独立行政法人国立高等専門学校機構秋田工業高等専門学校電気情報工学科
1
これらは，物理的には異なる問題だが，ポテンシャルや温度が満たす方程式は同じであ
る．方程式が同じならば解は同じで，同じ計算手法が使える．これらが満たす方程式は，
ラプラス方程式と呼ばれ，以下のような形をしている．
∇2 ϕ = 0
(1)
ϕ はポテンシャル (電圧) であったり温度で，(x, y, z) の 3 つの独立変数がある．ここでは，
三次元の問題は大変なので，図 1 や図 2 のように紙面の方向 (z 方向) には一様とする．そ
のため，z 方向の微分はゼロとなるので，ラプラス方程式は，
∂2 ϕ ∂2 ϕ
+
=0
∂x2 ∂y2
(2)
のように二次元問題になる．ここではこの偏微分方程式の近似解を数値計算により求める
のが目的である．プログラムが完成すると，図 3 のような解のグラフが得られる．
0℃
-20V
-20℃
0V
30V
30℃
図 1: ポテンシャルを求める問題．
図 2: 温度を求める問題．
2 差分法による偏微分方程式の数値計算
2.1
差分方程式
二次元のラプラス方程式を数値計算で近似解を求めることを考える．まずは，解 ϕ(x, y)
をテイラー展開する． x，y 方向に微小変位 ±h があった場合のポテンシャルは，
∂ϕ
h+
∂x
∂ϕ
ϕ(x − h, y) = ϕ(x, y) − h +
∂x
ϕ(x + h, y) = ϕ(x, y) +
1 ∂2 ϕ 2
h +
2! ∂x2
1 ∂2 ϕ 2
h −
2! ∂x2
2
1 ∂3 ϕ 3
h +
3! ∂x3
1 ∂3 ϕ 3
h +
3! ∂x3
1 ∂4 ϕ 4
h + ···
4! ∂x4
1 ∂4 ϕ 4
h − ···
4! ∂x4
(3)
(4)
30
20
10
0
-10
-20
0
0.2
0.4
0.6
0.8
1 0
0.2
0.4
0.6
0.8
1
図 3: 差分法により計算された，ポテンシャルや温度のグラフ．
となる．これらの式の辺々を足し合わせると，
]
∂2 ϕ 1 [
=
ϕ(x
+
h,
y)
−
2ϕ(x,
y)
+
ϕ(x
−
h,
y)
− O(h2 )
∂x2 x,y h2
(5)
が得られる．このことから，2 階の偏導関数の値は微小変位 h の場所の関数の値を用いて，
h2 の精度で近似計算ができることが分かる．つめり，式 (5) の右辺の第 1 項を計算すれば
よいのである．同じことを y 方向についても行うと
]
∂2 ϕ 1 [
2
(6)
= 2 ϕ(x, y + h) − 2ϕ(x, y) + ϕ(x, y − h) − O(h )
2
∂y x,y h
が得られる．
これらの式 (5) と (6) を元の 2 次元ラプラス方程式 (2) に代入すれば，
ϕ(x + h, y) + ϕ(x − h, y) + ϕ(x, y + h) + ϕ(x, y − h) − 4ϕ(x, y) = 0
(7)
となる．一般に微分方程式を上に示した引き算に直したものを差分式と呼ばれており，上
の式は 2 次元ラプラス方程式の差分式である．この式を眺めると，座標 (x, y) のポテンシャ
ルの値 ϕ(x, y) は，周りの値の平均であることがわかる．
実際にこの式を数値計算する場合，例えば図 1 のポテンシャルを求める時には，図 5 の
ように格子状 1 に区切り，その交点での値を求めることになる．ここでは， x および y 方
1
この格子のことをメッシュ(mesh) と言う事もある．
3
向には等間隔 h で区切り計算を進めるが，等間隔である必要はない．そのようにしても差
分の方程式は複雑になるが，同じような計算は可能である．これまでの説明が理解できて
いれば，x と y 方向の間隔が異なっても，式 (7) に対応する差分の式が作れるはずである．
Ny 等分
h
h
座標 (0, 0)
(0,0)=U00
Nx 等分
図 4: 解くべき領域を格子に分割
数値計算をする場合，ϕ(x, y) や ϕ(x + h, y) の形は不便なので，配列を使った表現に改め
る．図 4 の左下の座標を (0,0) として，格子点でのポテンシャルを
ϕ(x, y) = ϕ(ih, jh)
= Ui j
(8)
とする．このようにすると，式 (7) は
Ui+1 j + Ui−1 j + Ui j+1 + Ui j−1 − 4Ui j = 0
(9)
となり，数値計算し易い形になる．このようにした場合の各格子点の様子を図 5 に示す．
次の節で述べる境界条件を考えないとすると，ラプラス方程式は式 (9) の連立方程式を
解くだけである．格子に領域を分割することにより，難しい偏微分方程式が連立方程式の
問題に置き換えられたわけである．
3 練習問題
• 二次元直交座標系のラプラス方程式から差分式を求めなさい．
4
h
h
h
h
ui-1 j+2
ui j+2
ui+1 j+2
ui+2 j+2
ui-2 j+1
ui-1 j+1
ui j+1
ui+1 j+1
ui+2 j+1
ui-2 j
ui-1 j
uij
ui+1 j
ui+2 j
ui-2 j-1
ui-1 j-1
ui j-1
ui+1 j-1
ui+2 j-1
ui-2 j-2
ui-1 j-2
ui j-2
ui+1 j-2
ui+2 j-2
h
h
h
h
ui-2 j+2
図 5: 差分の格子
5

Download Report