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

数値計算法
第13回
偏微分方程式と差分法
2008年度 九州工業大学工学部電気電子工学コース用講義資料
講師：趙孟佑
[email protected]
1
偏微分方程式
•  関数が二つ以上の独立変数に依存する時
–  例えば
•  時間tと距離x
•  距離xと距離y
•  距離x,距離y,距離z
•  時間tと距離x,距離y,距離z
•  普通は最大４つまで
–  その関数の独立変数の変化に対応した振るまい
を記述する。
2
偏微分方程式の例（楕円形）
•  Laplace方程式
–  電位f(x,y)の空間分布
f=0
等電位線�
4角形の中の�
電位分布�
f=10
f=0
f=0
f=0
3
偏微分方程式の例（楕円形）
•  Poisson方程式
–  電位f(x,y)の空間分布
r(x,y)は空間電荷�
等電位線�
f=0
+電荷�
f=0
-電荷�
f=0
4角形の中の�
電位分布�
f=0
4
偏微分方程式の例（放物形）
•  熱伝導方程式
–  温度T(x,t)の時間・空間分布
r:密度、C：比熱、k:熱伝導率�
T
時間経過�
x=0
x=L
5
偏微分方程式の例（双曲形）
•  波動方程式
–  電磁波による電界E(x,t)の伝搬
c:光速�
x
前進波�
t
後進波�
6
差分法
•  微分を空間（あるいは時間）の連続的な領
域でするのでなく、ある離散的な点と点の
間の傾きで近似
f(x)のxiでの接線f’(xi)
後退差分�
前進差分�
中心差分�
xi-1
xi
xi+1
7
後退差分
•  関数f(x)のxiでの解析的な微分はf’(xi)
•  これをxiとxi-1の間の直線的な傾きで近似す
る。
これを後退差分という。�
8
後退差分
•  f(xi-1)のxiの値を使ってTaylor展開で求めると、
但し、�
に代入すると、�
誤差�
9
後退差分
•  よって、後退差分を微分の近似に使うと、誤差
は
となり、刻み幅Dxの1乗に比例する。�
1次の精度と言う。�
10
前進差分
•  関数f(x)のxiでの解析的な微分はf’(xi)
•  これをxi+1とxiの間の直線的な傾きで近似す
る。
これを前進差分という。�
11
前進差分の精度
•  f(xi+1)のxiの値を使ってTaylor展開で求めると、
但し、�
に代入すると、�
誤差�
12
前進差分の精度
•  よって、前進差分の誤差も
となり、刻み幅Dxの1乗に比例する。�
13
中心差分
•  関数f(x)のxiでの解析的な微分はf’(xi)
•  これをxi+1とxi-1の間の直線的な傾きで近似す
る。
これを中心差分という。�
14
中心差分の精度
•  f(xi+1)とf(xi-1)のxiの値を使ってTaylor展開で求め
ると、
差をとって、2Dxで割ると、�
誤差�
15
中心差分の精度
•  中心差分の誤差は
となり、刻み幅Dxの2乗に比例する。�
2次の精度と言う。�
16
差分法の応用
•  微分方程式を差分法を用いて解く。
•  例
–  1次元Poisson方程式の解法
•  電磁気（静電気）で多用する。
17
Poisson方程式の物理的意味
•  Gaussの定理
電束密度�
電界強度�
18
Poisson方程式の物理的意味
•  電界は電位の傾き
マイナス(-)がつくことに注意�
f
x
fの傾き負、電界Exは正�
19
Poisson方程式の物理的意味
•  Poisson方程式
に�
を代入する�
電位の2階微分=-電荷/誘電率�
20
Poisson方程式の物理的意味
f(x=0)=V
f(x=L)=0
x
V
簡単化のため1次元で考える�
中に電荷がないと、r=0
電界Exは一定�
21
Poisson方程式の物理的意味
f(x=0)=V
f(x=L)=0
x
f
V
V
Ex
0
L
22
Poisson方程式の物理的意味
f(x=0)=V
f(x=L)=0
+ +
x
+ +
V
＋の電荷が中にあると�
電位は上に凸になる�
23
Poisson方程式の物理的意味
f(x=0)=V
f(x=L)=0
+ +
x
+ +
f
V
V
Ex
0
L
24
Poisson方程式の物理的意味
f(x=0)=V
f(x=L)=0
ー� ー�
ー� ー�
x
V
ーの電荷が中にあると�
電位は下に凸になる�
25
Poisson方程式の物理的意味
f(x=0)=V
f(x=L)=0
ー� ー�
ー� ー�
f
x
V
V
Ex
0
L
26
Poisson方程式の物理的意味
f(x=0)=V
f(x=L)=0
ー� ー�
+ +
ー� ー�
x
+ +
V
＋とーの電荷がわかれてあると�
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（ディリクレ）型の境界条件
  傾き（一階微分）を与える
  例：
  Neuman（ノイマン）型の境界条件
  値か、傾きかのどちらか一方を与えればよい。両
方は要らない。
29
1次元Poisson方程式の解法
結局、1次元のPoisson方程式�
をx=0からx=Lの間で差分法を使って解くとは�
境界x=0とx=LでのΦの値と、x=0~x=LのΔx刻みのN+1個の点�
x1,x2,x3,-----,xN,xN+1での空間電荷密度r(x1),r(x2),r(3),--,r(xN),r(xN+1)
が分かっていて、�
x2,x3,-----,xN-1,XNでの電位fを求める�
x1
0
xN
x2
DX
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 xi+1
xi+2
xiでのfの2階微分を次式で近似�
これはまず、i-1/2とi+1/2での一階微分を中心差分で定義して�
（計算する必要はない）�
32
2階微分の差分
fi+1
fi
fi-1
fi-2
xi-1 xi-1/2 xi xi+1/2 xi+1
xi-2
を�
fi+2
と�
xi+2
の中心差分で表す�
33
2階微分の差分の誤差
より�
誤差�
よって、この差分のやり方は2次の精度�
34
1次元Poisson方程式の解法
は以下の連立方程式になる�
N-1個�
i=1,i=N+1での�
式は要らない。�
Dirichletの境界条件だから�
35
1次元Poisson方程式の解法
36
1次元Poisson方程式の解法
行列に直すと、�
係数行列はN-1XN-1の対角対称行列で、�
対角要素とその前後以外は全て0
消去法で解くのが最も手軽�
37
1次元Poisson方程式の解法
ここで�
38
1次元Poisson方程式の解法
まず2行目から対角要素の下の要素を順繰りに消していく。�
①
①-②x(-2)=②’
②
②’
②’-③x②’の係数=③’
③
③’
以下同様に�
④
39
1次元Poisson方程式の解法
2行目からN-1行目まで、対角要素の下の要素を消すと、�
あとは、N-1行目から上に解いていく。�
40
1次元Poisson方程式の解法
41
1次元Poisson方程式の例
境界条件は�
刻み幅Dx=0.1m
42
1次元Poisson方程式の例
電荷の無いとき�
43
例
x
r(C/m3)
0.0��0
1
1
2
2
3
3
4
2
5
0
6
- 2
7
- 3
8
-2
9
-1
10��0
44
例
x
r(C/m3)
0.0��0
1
1
2
2
3
3
4
2
5
0
6
- 2
7
- 3
8
-2
9
-1
10��0
式は全部で9個�
45
例
x
r(C/m3)
0.0��0
1
1
2
2
3
3
4
2
5
0
6
- 2
7
- 3
8
-2
9
-1
10��0
46
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
c
phi(1)=phil
phi(n+1)=phir
c
e(1)=-dx*dx/eps0*rho(2)-phil
c
e(n-1)=-dx*dx/eps0*rho(n-1)-phir
� e(n-1)=-dx*dx/eps0*rho(n)-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
47
宿題
x=0からx=Lまでの一次元のpoisson方程式を解いて、�
x=0からx=Lまでの電位fをグラフでプロットせよ。�
L=1m, Dx=0.1m,f(x=0)=100,f(x=L)=0とし、rの分布は次ページ�
の表で与えられる。�
このmはメートルの意味�
48
宿題
x
r(C/m3)
0.0��0
0.1
1x10-9
0.2
2x10-9
0.3
(m+1)x10-9
0.4
2x10-9
0.5
0
0.6
- 2x10-9
0.7
- 3x10-9
0.8
-(m+1)x10-9
0.9
-1x10-9
1.0��0
ここでmは学籍番号の下一桁である。�
49

Download Report