Transcript 数値計算法I - 九州工業大学

数値計算法
第11回
偏微分方程式と差分法
2006年度
九州工業大学工学部電気電子工学コース用講義資料
講師：趙孟佑
[email protected]
1
偏微分方程式
• 関数が二つ以上の独立変数に依存する時
– 例えば
•
•
•
•
•
時間tと距離x
距離xと距離y
距離x,距離y,距離z
時間tと距離x,距離y,距離z
普通は最大４つまで
– その関数の独立変数の変化に対応した振るまい
を記述する。
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 (xi1 )
f '(xi ) 
xi  xi1
これを後退差分という。
8
後退差分
• f(xi-1)のxiの値を使ってTaylor展開で求めると、
x2
x 3
f (xi1 )  f (xi )  xf (xi ) 
f (xi ) 
f (xi ) L
2!
3!
但し、 x  xi  xi1
f (xi )  f (xi1 )
に代入すると、
xi  xi1
x2
f (xi )  f (xi1 ) f (xi )  f (xi )  xf (xi )  2! f (xi ) L

xi  xi1
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次の精度と言う。
10
前進差分
• 関数f(x)のxiでの解析的な微分はf’(xi)
• これをxi+1とxiの間の直線的な傾きで近似す
る。
f (xi1 )  f (xi )
f '(xi ) 
xi1  xi
これを前進差分という。
11
前進差分の精度
• f(xi+1)のxiの値を使ってTaylor展開で求めると、
x2
x 3
f (xi1 )  f (xi )  xf (xi ) 
f (xi ) 
f (xi ) L
2!
3!
但し、 x  xi1  xi
f (xi1 )  f (xi )
に代入すると、
xi1  xi
x2
f (xi1 )  f (xi ) f (xi )  xf (xi )  2! f (xi ) L  f (xi )

xi1  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の１乗に比例する。
13
中心差分
• 関数f(x)のxiでの解析的な微分はf’(xi)
• これをxi+1とxi-1の間の直線的な傾きで近似
する。
f (xi1 )  f (xi1 ) f (xi1 )  f (xi1 )
f '(xi ) 

xi1  xi1
2x
これを中心差分という。
14
中心差分の精度
• f(xi+1)とf(xi-1)のxiの値を使ってTaylor展開で求
めると、
2
3
x
x
f (xi1 )  f (xi )  xf (xi ) 
f (xi ) 
f (xi ) L
2!
3!
x2
x 3
f (xi1 )  f (xi )  xf (xi ) 
f (xi ) 
f (xi ) L
2!
3!
差をとって、2xで割ると、
x 3
f (xi1 )  f (xi1 ) 2xf (xi )  2 3! f (xi ) L

2x
2x
誤差
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 x0
x
0
xL
 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
fi1  2fi  fi1

2
x i
x2
これはまず、i-1/2とi+1/2での一階微分を中心差分で定義して
（計算する必要はない）
f
fi  fi1

,
x i1/2
x
f
fi1  fi

x i1/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 i1/2
x i1/2
f
f

2
x i1/2 x i1/2
f

2
x i
x
33
2階微分の差分の誤差
f x2 2f x 3 3f
fi1  fi  x


L
2
3
x i 2! x i 3! x i
f x2 2f x 3 3f
fi1  fi  x


L
2
3
x i 2! x i 3! x i
より
誤差
2f
fi1  2fi  fi1 2f 2x2 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個
fi1  2fi  fi1
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
fi1  2fi  fi1  
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  fi1
e1  
ei  
x2
0
x2
0
eN 1  
2  f1
i1
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
①
②
②’
③
③’
④
bi1
bi1
bi1






bi  ci1
bi , ci  
ci , ei  ei1
ei
 
 
ai
ai
ai
bi  ci1
  2bi1
 , ci  bi1
 , ei  ei1
  bi1
 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  cidi1 
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