数値計算法I - 九州工業大学
Download
Report
Transcript 数値計算法I - 九州工業大学
数値計算法
第11回
偏微分方程式と差分法
2006年度
九州工業大学工学部電気電子工学コース用講義資料
講師:趙孟佑
[email protected]
1
偏微分方程式
• 関数が二つ以上の独立変数に依存する時
– 例えば
•
•
•
•
•
時間tと距離x
距離xと距離y
距離x,距離y,距離z
時間tと距離x,距離y,距離z
普通は最大4つまで
– その関数の独立変数の変化に対応した振るまい
を記述する。
2
偏微分方程式の例(楕円形)
• Laplace方程式
– 電位f(x,y)の空間分布
2f 2f
2 0
2
x y
f=0
等電位線
4角形の中の
電位分布
f=10
f=0
f=0
f=0
1.0
2.0
3.0
5.0
3
偏微分方程式の例(楕円形)
• Poisson方程式
– 電位f(x,y)の空間分布
(x,y)は空間電荷
2f 2f
2
2
x y
o
等電位線
f=0
+電荷
f=0
-電荷
f=0
4角形の中の
電位分布
f=0
4
偏微分方程式の例(放物形)
• 熱伝導方程式
– 温度T(x,t)の時間・空間分布
T
2T
C
2
t
x
:密度、C:比熱、:熱伝導率
T
時間経過
x=0
x=L
5
偏微分方程式の例(双曲形)
• 波動方程式
– 電磁波による電界E(x,t)の伝搬
2 E 2 2 E
c
2
t
x 2
c:光速
x
前進波
t
後進波
6
差分法
• 微分を空間(あるいは時間)の連続的な領
域でするのでなく、ある離散的な点と点の
間の傾きで近似
f(x)のxiでの接線f’(xi)
後退差分
前進差分
中心差分
xi-1
xi
xi+1
7
後退差分
• 関数f(x)のxiでの解析的な微分はf’(xi)
• これをxiとxi-1の間の直線的な傾きで近似す
る。
f (xi ) f (xi1 )
f '(xi )
xi xi1
これを後退差分という。
8
後退差分
• f(xi-1)のxiの値を使ってTaylor展開で求めると、
x2
x 3
f (xi1 ) f (xi ) xf (xi )
f (xi )
f (xi ) L
2!
3!
但し、 x xi xi1
f (xi ) f (xi1 )
に代入すると、
xi xi1
x2
f (xi ) f (xi1 ) f (xi ) f (xi ) xf (xi ) 2! f (xi ) L
xi xi1
x
2
x
xf (xi )
f (xi ) L
2!
x
x
誤差
f (xi )
f (xi ) L
2!
9
後退差分
• よって、後退差分を微分の近似に使うと、誤
差は
x
x 2
f (xi )
f (xi ) L
2!
3!
となり、刻み幅xの1乗に比例する。
1次の精度と言う。
10
前進差分
• 関数f(x)のxiでの解析的な微分はf’(xi)
• これをxi+1とxiの間の直線的な傾きで近似す
る。
f (xi1 ) f (xi )
f '(xi )
xi1 xi
これを前進差分という。
11
前進差分の精度
• f(xi+1)のxiの値を使ってTaylor展開で求めると、
x2
x 3
f (xi1 ) f (xi ) xf (xi )
f (xi )
f (xi ) L
2!
3!
但し、 x xi1 xi
f (xi1 ) f (xi )
に代入すると、
xi1 xi
x2
f (xi1 ) f (xi ) f (xi ) xf (xi ) 2! f (xi ) L f (xi )
xi1 xi
x
2
x
xf (xi )
f (xi ) L
2!
x
x
12
誤差
f (xi )
f (xi ) L
2!
前進差分の精度
• よって、前進差分の誤差も
x
x2
f (xi )
f (xi ) L
2!
3!
となり、刻み幅xの1乗に比例する。
13
中心差分
• 関数f(x)のxiでの解析的な微分はf’(xi)
• これをxi+1とxi-1の間の直線的な傾きで近似
する。
f (xi1 ) f (xi1 ) f (xi1 ) f (xi1 )
f '(xi )
xi1 xi1
2x
これを中心差分という。
14
中心差分の精度
• f(xi+1)とf(xi-1)のxiの値を使ってTaylor展開で求
めると、
2
3
x
x
f (xi1 ) f (xi ) xf (xi )
f (xi )
f (xi ) L
2!
3!
x2
x 3
f (xi1 ) f (xi ) xf (xi )
f (xi )
f (xi ) L
2!
3!
差をとって、2xで割ると、
x 3
f (xi1 ) f (xi1 ) 2xf (xi ) 2 3! f (xi ) L
2x
2x
誤差
x 2
f (xi )
f (xi ) L
3!
15
中心差分の精度
• 中心差分の誤差は
x2
x 4 (5)
f (xi )
f (xi ) L
3!
5!
となり、刻み幅xの2乗に比例する。
2次の精度と言う。
16
差分法の応用
• 微分方程式を差分法を用いて解く。
• 例
– 1次元Poisson方程式の解法
• 電磁気(静電気)で多用する。
2f
2
x
o
17
Poisson方程式の物理的意味
• Gaussの定理
r
D
D 電束密度
r
D 0 E
E 電界強度
r
r
Ey
Ez
Ex
D 0 E 0
0
0
x
y
z
18
Poisson方程式の物理的意味
• 電界は電位の傾き
E f
f
f
f
Ex , Ey , Ez
x
y
z
マイナス(-)がつくことに注意
f
x
fの傾き負、電界Exは正
19
Poisson方程式の物理的意味
• Poisson方程式
に
Ex Ey Ez
0
x
y
z
E f を代入する
2
2
2
f f f
2
2
2
x y z
0
電位の2階微分=-電荷/誘電率
20
Poisson方程式の物理的意味
f(x=0)=V
f(x=L)=0
x
V
簡単化のため1次元で考える
2f
2
x
0
中に電荷がないと、=0
f
0
2
x
2
Ex
0
x
電界Exは一定
21
Poisson方程式の物理的意味
f(x=0)=V
f(x=L)=0
x
f
V
V
V
Ex
L
Ex
0
L
22
Poisson方程式の物理的意味
f(x=0)=V
f(x=L)=0
+
+
+
x
+
V
+の電荷が中にあると
2f
0
2
x
0
電位は上に凸になる
23
Poisson方程式の物理的意味
f(x=0)=V
f(x=L)=0
+
+
f
+
x
+
V
V
Ex
0
L
24
Poisson方程式の物理的意味
f(x=0)=V
f(x=L)=0
ー ー
x
ー ー
V
ーの電荷が中にあると
2f
0
2
x
0
電位は下に凸になる
25
Poisson方程式の物理的意味
f(x=0)=V
f(x=L)=0
ー ー
x
ー ー
f
V
V
Ex
0
L
26
Poisson方程式の物理的意味
f(x=0)=V
f(x=L)=0
ー ー
+
+
ー ー
x
+
+
V
+とーの電荷がわかれてあると
2f
2
x
0
27
Poisson方程式の物理的意味
f(x=0)=V
f
f(x=L)=0
ー ー
+
+
ー ー
x
+
+
V
•誘電体の分極
•プラズマ中の電界
等で重要
V
Ex
0
L
28
楕円形微分方程式の境界条件
• Poisson, Laplace方程式を解くには、計算する領域
の境界での値についての条件が必要
値を与える
例:f(x=0)=V,f(x=L)=0
Dirichlet(ディリクレ)型の境界条件
傾き(一階微分)を与える
例:
f
f
0,
x x0
x
0
xL
Neuman(ノイマン)型の境界条件
値か、傾きかのどちらか一方を与えればよい。両
方は要らない。
29
1次元Poisson方程式の解法
結局、1次元のPoisson方程式
f
2
x
o
2
をx=0からx=Lの間で差分法を使って解くとは
境界x=0とx=LでのΦの値と、x=0~x=LのΔx刻みのN+1個の点
x1,x2,x3,-----,xN,xN+1での空間電荷密度(x1),(x2),(3),--,(xN),(xN+1)
が分かっていて、
x2,x3,-----,xN-1,XNでの電位fを求める
x1
0
xN
x2
X
xN+1
L
30
捕捉
なぜ、x0,x1,x2,x3,-----,xNでなく、 x1,x2,x3,-----,xN,xN+1とxに
番号をうつのか?
Fortran77ではx(0)という配列の添字は使えない。
よって、n+1個の数の配列を表すには
x(1),x(2),----,x(n),x(n+1)と1から始めないといけない。
ここで、x(1)=0,x(n+1)=Lとなる。
CやJAVAではx(0)という配列の添字が使えるので、
x0,x1,x2,x3,-----,xNとしてもよい。
31
2階微分の差分
fi+1
fi-2
xi-2
fi-1
fi+2
fi
xi-1 xi-1/2 xi
xi+1/2 x
i+1
xi+2
xiでのfの2階微分を次式で近似
2f
fi1 2fi fi1
2
x i
x2
これはまず、i-1/2とi+1/2での一階微分を中心差分で定義して
(計算する必要はない)
f
fi fi1
,
x i1/2
x
f
fi1 fi
x i1/2
x
32
2階微分の差分
fi+1
fi-1
fi-2
xi-2
fi+2
fi
xi-1 xi-1/2 xi
xi+1/2 x
i+1
xi+2
2f
f
f
の中心差分で表す
を
と
2
x i
x i1/2
x i1/2
f
f
2
x i1/2 x i1/2
f
2
x i
x
33
2階微分の差分の誤差
f x2 2f x 3 3f
fi1 fi x
L
2
3
x i 2! x i 3! x i
f x2 2f x 3 3f
fi1 fi x
L
2
3
x i 2! x i 3! x i
より
誤差
2f
fi1 2fi fi1 2f 2x2 4f
2
L
2
2
4
x i
x
x i
4! x i
よって、この差分のやり方は2次の精度
34
1次元Poisson方程式の解法
f
2
x
o
2
は以下の連立方程式になる
f1 2f2 f3
2
2
x
0
f2 2f3 f4
3
2
x
0
M
M M
N-1個
fi1 2fi fi1
i
2
x
0
M
M M
fN 1 2fN fN 1
N
x2
0
i=1,i=N+1での
式は要らない。
Dirichletの境界条件だから
35
1次元Poisson方程式の解法
2f2 f3
f2 2f3 f4
M
M M
fi1 2fi fi1
M
fN 1 2fN
x 2
0
x 2
0
x 2
0
2 f1
3
i
M M
x 2
0
N fN 1
36
1次元Poisson方程式の解法
x2
2 f1
0
x2
2 1 0 L 0 0 f2
1 2 1 O O 0 f 3
0
3
0 O O O O M M
M
2
0 O 1 2 1 0 fi x
i
M O O O O O M
0
M
0 0 L 0 1 2 fN
2
x
係数行列はN-1XN-1の対角対称行列で、
f
N
N 1
対角要素とその前後以外は全て0
0
行列に直すと、
消去法で解くのが最も手軽
37
1次元Poisson方程式の解法
b1
a
2
0
0
M
0
ここで
c1
b2
O
O
O
0
0
c2
O
ai
O
L
L
0
0 d1 e1
O O
0 d2 e2
O O
M M M
bi
ci
0 di ei
O O
O M M
0 aN 1 bN 1 dN 1 eN 1
ai 1
bi 2
ci 1
di fi1
e1
ei
x2
0
x2
0
eN 1
2 f1
i1
x2
0
N fN 1
38
1次元Poisson方程式の解法
まず2行目から対角要素の下の要素を順繰りに消していく。
2d1 d2
e1
d1 2d2 d3 e2
①-②x(-2)=②’
3d2 2d3
e1 2e2
d2 2d3 d4 e3
②’-③x②’の係数=③’
以下同様に
4d3 3d4
e1 2e2 3e3
d3 2d4 d5 e4
①
②
②’
③
③’
④
bi1
bi1
bi1
bi ci1
bi , ci
ci , ei ei1
ei
ai
ai
ai
bi ci1
2bi1
, ci bi1
, ei ei1
bi1
ei
39
1次元Poisson方程式の解法
2行目からN-1行目まで、対角要素の下の要素を消すと、
b1
0
0
0
M
0
c1
b2
O
O
O
0
0
c2
O
0
O
L
L
O
O
bi
O
0
0
0 d1 e1
O
0 d2 e2
O
M M M
ci 0 di ei
O O M M
0 bN 1 dN 1 eN 1
あとは、N-1行目から上に解いていく。
eN 1
dN 1
bN 1
1
dN 2
eN 2 cN 2dN 1
bN 2
40
1次元Poisson方程式の解法
eN 1
dN 1
bN 1
1
dN 2
eN 2 cN 2dN 1
bN 2
M
1
di ei cidi1
bi
M
1
d1 e1 c1d2
b1
41
1次元Poisson方程式の例
3
(C/m )
5 10-9
0 100
-5 10-9
0
0.2
0.6
0.4
0.8
1
x (m)
境界条件は f(x 0) 100(V)
f(x 1) 0
刻み幅x=0.1m
42
1次元Poisson方程式の例
120
100
f (V)
80
電荷の無いとき
60
40
20
0
0
0.2
0.4
0.6
0.8
1
x (m)
43
1次元Poisson方程式の例
implicit real*8 (a-h,o-z)
parameter(nmax=11)
dimension phi(nmax)
dimension a(nmax),b(nmax),c(nmax)
dimension d(nmax),e(nmax)
dimension rho(nmax),phi(nmax)
dimension x(nmax)
n=10
dx=0.1
eps0=8.85d-12
phil=100
phir=0
c
open(unit=10,file='rho.dat',status='old')
do i=1,n+1
read(10,*)idum,rho(i)
end do
close(10)
c
do i=1,n+1
x(i)=(i-1)*dx
end do
c
do i=1,n-1
a(i)=1.d0
b(i)=-2.d0
c(i)=1.d0
end do
phi(1)=phil
phi(n+1)=phir
c
e(1)=-dx*dx/eps0*rho(2)-phil
e(n-1)=-dx*dx/eps0*rho(n-1)-phir
do i=2,n-2
e(i)=-dx*dx/eps0*rho(i+1)
end do
c
c Solution of tri-diagnola matrix
c
do i=2,n-1
b(i)=c(i-1)-b(i-1)/a(i)*b(i)
c(i)=-b(i-1)/a(i)*c(i)
e(i)=e(i-1)-b(i-1)/a(i)*e(i)
end do
c
d(n-1)=e(n-1)/b(n-1)
c
do i=n-2,1,-1
d(i)=(e(i)-c(i)*d(i+1))/b(i)
end do
c
do i=1,n-1
phi(i+1)=d(i)
end do
c
open(unit=20,file='poisson2.lis',status='unknown')
do i=1,n+1
write(20,110)x(i),phi(i)
end do
110 format(2(1x,e12.5))
close(20)
c
stop
end
c
44