筑波大学 工学基礎学類電子量子専攻3年次 『シミュレーション工学』

1
シミュレーション工学
筑波大学
工学基礎学類 電子量子専攻 ３年次
『シミュレーション工学』
筑波大学 物理工学系
佐野 伸行
Office: 3M 405
Phone: 6479
e-mail: [email protected]
2
目次
目次
1
2
3
数値解析の基礎
1.1 数値微分 . . . . . . . . . . . . . . . . . . .
1.1.1 差分公式 . . . . . . . . . . . . . . .
1.1.2 差分公式の幾何学的意味と導出 . .
1.2 数値積分 . . . . . . . . . . . . . . . . . . .
1.2.1 積分公式 . . . . . . . . . . . . . . .
1.2.2 特異性のある積分 . . . . . . . . . .
1.3 非線形方程式の解法 . . . . . . . . . . . .
1.4 微分方程式の解法 . . . . . . . . . . . . . .
1.4.1 常微分方程式の解法：陽解法 . . .
1.4.2 数値解の収束性と誤差 . . . . . . .
1.5 物理への応用 . . . . . . . . . . . . . . . .
1.5.1 物理的背景：２原子分子の振動運動
1.5.2 数値解析の基本アルゴリズム . . .
1.6 モンテカルロ法の基礎 . . . . . . . . . . .
1.6.1 モンテカルロ法の基礎概念 . . . . .
1.6.2 重み関数を用いたモンテカルロ法 .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
デバイスモデリング入門
2.1 デバイス特性シミュレーション：ポアソン方程式 . . . . . . . .
2.1.1 MOS キャパシター: 陽解法 . . . . . . . . . . . . . . . .
2.1.2 陰解法 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.3 ３重対角行列の再帰的な解法 . . . . . . . . . . . . . . .
2.2 デバイス特性シミュレーション：ドリフト拡散法 . . . . . . . .
2.2.1 半導体デバイスの基本方程式 . . . . . . . . . . . . . . .
2.2.2 pn ダイオードへの応用とドリフト拡散法のアルゴリズム
2.2.3 ドリフト拡散法における物理モデル . . . . . . . . . . . .
2.2.4 ドリフト拡散方程式の離散化 . . . . . . . . . . . . . . .
2.2.5 ドリフト拡散方程式の線形化 . . . . . . . . . . . . . . .
参考文献
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
3
3
3
6
9
9
11
14
17
17
25
29
29
31
34
34
36
.
.
.
.
.
.
.
.
.
.
41
42
42
44
46
48
48
50
51
54
56
59
1. 数値解析の基礎
1
3
数値解析の基礎
1.1
数値微分
数値微分、数値積分、方程式の数値解法の３つは、物理的問題を計算機で取り扱うときに欠
かせない基本的な演算である。
いま、独立変数 x の任意の値に対しての関数 f (x) が与えられたとする。微分では、この関
数 f の導関数を求めることであり、積分ではこの微分の逆演算を、方程式では f がゼロになる
x の値を探す、ということである。関数 f (x) が解析的に与えられていれば、微分や積分は解析
的に求めることができる場合もある。しかし、多くの場合、微分や積分を解析的に実行するこ
とはかなりの困難を伴う。最近では、Mathematica や Maple 等の数式処理ソフトウェアを使い
こなすことにより、かなりの困難は緩和されているが、未だ完璧ではない。例えば、関数 f (x)
が非常に複雑であり、公式を用いた微分や積分等が困難な場合や、関数 f (x) が解析的に与えら
れておらず、ある離散的な値 x にのみ知られている場合などがこれにあたる。そこで、自分で
プログラミングをすることを前提とした数値的な解法が必要になる。
数値解析により解を求めるということは、殆どの場合、近似解を求めることと等価である。
つまり、連続の x について定義される関数 f (x) を離散的な値 xi での関数値 f (xi ) で近似し、方
程式を差分方程式に置き換えることになる。勿論、離散的な値 xi の間隔が狭まれば、厳密な連
続的な関数 f (x) に近づく。従って、数値解析での中心的な考え方は、関数 f をできるだけ良い
精度でこれらの演算が行いやすいように近似するということにある。
本書の最初の３章では、この３種類の基本的な演算を計算機上で実行する方法について解説
する。これらの手法は、今後の講義内容および数値解析全般を理解するうえで必要になる基礎
的素養とみなされるものである。より詳細な一般的説明は参考文献に挙げてある物理数学およ
び数値解析の専門書を参考にされたい。
1.1.1
差分公式
数値解析における１階微分
関数 f (x) が与えられたとき、いま、x = 0 における導関数 f (0) を求めるとしよう。
一般に、導関数 f (0) に対する差分公式は、関数 f (x) を原点近傍 x ≈ 0 で Taylor 展開する
ことにより導くことができる。いま、x 軸上での等間隔ごとの関数 f (x) の値がわかっているも
のとしよう。つまり、
fn = f (xn );
xn = nh (n = · · · , −2, −1, 0 + 1, +2, · · · )
の値が既知であるとする（図 1 参照）。そこで、x = 0 における導関数 f (0) を fn を用いて表
わすことを考える。
f (x) を x = 0 の近傍で Taylor 展開すれば、
f (x) = f (0) + x f (0) +
x2 x3
f (0) + f (0) + . . .
2!
3!
(1)
1. 数値解析の基礎
4
f-1
f0
-h
0
f-2
-2h
f(x)
f1
f2
h
2h
x
図 1: 関数 f (x) の離散化
うえの式において、x = ±h とすれば、
f±1 ≡ f (±h) = f (0) ± h f (0) +
h2 h3
f (0) ± f (0) + O(h4 )
2
6
ここで、O(h4 ) は h4 以上の高次の項を表わす。
２点公式
うえの式から、f0 および f±1 の２点を用いることにより、導関数 f (0) を近似することを考
える。この場合の２点公式（2-point formula）は、h について２次以上の項を無視することに
より、以下のように定義される。
f (0) =
f1 − f0
+ O(h) （前進２点公式）
h
(2)
f0 − f−1
+ O(h) （後退２点公式）
(3)
h
これらの２点公式は、次のセクションで説明するように、区間 [0, ±h] における f の線形近似
に基づいている。しかしながら、この場合の誤差は O(h) のオーダーであり、誤差が O(h2 ) で
ある以下の３点公式に比べて、精度はかなり悪い。
f (0) =
３点公式
導関数 f (0) を f1 と f−1 を使って表わすとすると、
f1 − f−1 = 2hf (0) +
である。従って、
f (0) =
h3 f (0) + O(h4 )
3
f1 − f−1 h2 − f (0) + O(h3 )
2h
6
1. 数値解析の基礎
5
２次の項はキャンセルしていることに注意せよ。もし間隔 h が十分に小さいとするならば、上
式の第２項以降は無視できるであろう。結局、x = 0 における導関数 f (0) は、以下の３点公式
（3-point formula）で与えられる。
f (0) =
f1 − f−1
+ O(h2 )
2h
(4)
一般に、間隔 h を小さくすればするほど精度は高くなる。もし関数 f (x) が二次までの多項
式で与えられるならば３点公式は厳密になる。言い換えれば、３点公式は、３点 x = ±h, 0 の
関数値をもとにした２次多項式で関数 f を近似することを意味する（セクション 1-3 参照）。
５点公式
３点公式よりもさらに精度を上げるためには、x = 0 からさらに離れた点における関数値（例
えば、f±2 ）を用いて、導関数 f (0) を近似すれば良い。
Tylor 展開の式 (1) において、x = ±2h とすれば以下の式を得る。
f±2 ≡ f (±2h) = f (0) ± 2h f (0) + 2h2 f (0) ±
4h3 f (0) + O(h4 )
3
３点公式を求めたとときと全く同様の手順をふめば、
f−2 − f2 = −4hf (0) −
8h3 f (0) + O(h5 )
3
8h3 f (0) + O(h5 )
3
であるから、以下のような５点公式（5-point formula）を得る。
8(f1 − f−1 ) = 16hf (0) +
f (0) =
f−2 − 8f−1 + 8f1 − f2
+ O(h4 )
12h
(5)
３次および４次の項がキャンセルしていることに注意せよ。
５点公式は、５点区間 [−2h, 2h] における f の４次多項式の近似に基づいており、h に関し
て３次のオーダーまで有効である。ただし、場合によって、５点公式は、３点公式に比べて各
項の打ち消し合いの起きることが多いために、幅 h が大きいときにその精度に注意が必要であ
る。このようなケタ落ちによる誤差をまるめ誤差（round-off error）という。
任意の場所 x = a における導関数 f (a) を求めるには、うえの議論における x を a だけ平行
移動すれば良い。つまり、f (x) を x = a の近傍で Taylor 展開すれば、
f (x) = f (a) + (x − a) f (a) +
(x − a)2 (x − a)3 f (a) +
f (a) + . . .
2!
3!
ここで、新しい座標 x を x = x − a とすれば、うえの展開式は f (x) の x = 0 の近傍での Taylor
展開式と一致する。即ち、任意の場所 x = a における導関数 f (a) を求めるためには、うえで
求めたすべての公式において、
f0 → f (a),
f±1 → f (a ± h),
···
1. 数値解析の基礎
6
と読み代えてやれば良い。
数値解析における高次の微分
高次の微分に対する公式も、これまでの１階微分の場合と同様にして求めることができる。
例えば、テイラー展開の式 (1) を用いれば、
f1 − 2f0 + f−1 = h2 f (0) + O(h4 )
であるから、x = 0 における二次の導関数 f (0) は、以下の公式で与えられる。
f (0) =
f−1 − 2f0 + f1
+ O(h2 )
h2
(6)
うえの公式は、一次の導関数に対する３点公式を用いて二次の導関数 f (0) を表したものと
等価である。すなわち、メッシュ間隔を半分 h/2 に想定したうえで、３点公式を２回用いるこ
とにより、二次の導関数 f (0) を得ることができる。
f (0) =
=
1 f1 − f0 f0 − f−1
−
h
h
h
f1/2 − f−1/2
h
このようなトリックは、メッシュを用いて高階の（偏あるいは常）微分方程式を解くことが必
要な理工学の各分野で用いられる。
1.1.2
差分公式の幾何学的意味と導出
これまでのセクションで求めた微分に対する差分公式は、かなり天下り的な印象を受ける。
そこで、ここでは、これらの公式の意味を理解するとともに、より系統的に差分公式を導出す
ることを試みよう。
３点公式
いま、関数 f (x) が x = 0 近傍で x についての２次の多項式で近似できたと仮定する。
f (x) = c0 + c1 x + c2 x2
ここで、c0 , c1 , c2 は定数である。x = ±h, 0 の３点での関数値がそれぞれ f±1 , f0 であるとすれ
ば、これら既知の関数値を用いて定数 c0 , c1 , c2 を決めることができる。
f0 = f (x = 0) = c0
f1 = f (x = h) = c0 + c1 h + c2 h2
f−1 = f (x = −h) = c0 − c1 h + c2 h2
1. 数値解析の基礎
7
これらを解くと、
f1 − f−1
f1 − 2f0 + f−1
, c2 =
2h
2h2
従って、x = 0 近傍で関数 f (x) は、以下のように近似される。
c0 = f0 ,
c1 =
f (x) = f0 +
f1 − 2f0 + f−1 2
f1 − f−1
x+
x
2h
2h2
微分をすることにより、以下のような１階微分と２階微分の３点公式を得る。
f (0) =
f1 − f−1
2h
f1 − 2f0 + f−1
h2
即ち、３点公式は２次曲線で f (x) を近似することに対応している。
f (0) =
５点公式
全く同様にして、高次の微分公式も導くことができる。例えば、５点公式は、f (x) を４次曲
線で近似することと等価である。そこで、
f (x) = c0 + c1 x + c2 x2 + c3 x3 + c4 x4
とおいて、この曲線が f (x = 0) = f0 , f (x = ±h) = f±1 , f (x = ±2h) = f±2 を通ること
から、






これを解いて、

2h 4h2
8h3 16h4
c1

2
3
4

h h
h
h 
  c2

2
3
4
−h h
−h
h 
 c3
2
3
4
−2h 4h −8h 16h
c4






c1
c2
c3
c4


 
 
=
 
 


 
 
=
 
 
f2 − f0
f1 − f0
f−1 − f0
f−2 − f0
f−2 −8f−1 +8f1 −f2
12h
0 −16f1 +f2
− f−2 −16f−1 +30f
24h2
+2f1 −2f2
− f−2 −2f−1
12h3
f−2 −4f−1 +6f0 −4f1 +f2
24h4












従って、例えば、１階微分の５点公式はまえに導出したものと一致する。
f (0) = c1 =
f−2 − 8f−1 + 8f1 − f2
12h
その他の高階微分の５点公式は以下のようになる。
f (0) = 2c2 = −
f−2 − 16f−1 + 30f0 − 16f1 + f2
12h2
f (0) = 6c3 = −
f−2 − 2f−1 + 2f1 − 2f2
2h3
1. 数値解析の基礎
8
f (4) (0) = 24c4 =
演習問題 1-1-1
つぎの関数
f−2 − 4f−1 + 6f0 − 4f1 + f2
h4
x
f (x) = xx
に対して、f (x = 1.5) の値を２点公式（前進および後退）、３点公式、５点公式を用いて評価す
るプログラムを作成せよ。きざみ幅 h は、h = 0.5, 0.2, 0.1, 0.05, 0.02, 0.01, 0.005, 0.002, 0.001,
0.0005, 0.0002, 0.0001, 0.00005, 0.00002, 0.00001, 0.000005, 0.000002, 0.000001 とせよ。
f (x) を１階微分することにより、f (x) を解析的に求めよ。そして、f (x = 1.5) に対する厳
密値を求めたうえで、それぞれの公式から求めた数値解との誤差を求めて公式の精度について
議論せよ。
演習問題 1-1-2
x = 0 近傍における関数 f (x) を４次曲線で近似することにより、f∓1 , f0 , f±1 , f±2 を使って、
f (0) の４点公式を導け。
（ヒント：うえの複合記号の一方をもちいることにより、公式は２通り存在する）。
1. 数値解析の基礎
1.2
9
数値積分
ある関数 f (x) の a < x < b における定積分を求めるとしよう。
a, b のあいだを偶数個 N の等間隔の格子に分割したとする。つまり、
N=
(b − a)
h
である。ここで、N は偶数。つぎに、求めたい定積分を以下のように変形する。
b
a+2h
f (x)dx =
a
a+4h
f (x)dx+
a
a+2h
f (x)dx+ · · · +
b
f (x)dx
(7)
b−2h
上式の各項を評価する際に、用いる近似法によっていくつかの積分公式が存在する。ここでは、
最も一般的な台形公式とシンプソン公式について学ぶ。
1.2.1
積分公式
台形公式（trapezoidal formula）
以下のような定積分
h
−h
f (x)dx
を考える際、区間 [−h, 0] と区間 [0, h] における関数 f を前セクションの数値微分で求めた１次
差分式（２点公式）でそれぞれおきかえてみる。
区間 [−h, 0] においては、
f (x) ≈ f (0) + f (0)x ≈ f0 +
従って、
0
−h
f (x)dx ≈
0
−h
dx f0 +
f0 − f−1
x
h
f0 − f−1
h
x = (f−1 + f0 )
h
2
区間 [0, h] においても同様にして、
h
0
f (x)dx ≈
h
0
dx f0 +
f1 − f0
h
x = (f1 + f0 )
h
2
従って、区間 [−h, h] に対しては、
h
h
f (x)dx = (f−1 + 2f0 + f1 ) + O(h3 )
2
−h
を得る。これは、数値積分の近似式として最も良く知られており、台形公式（trapezoidal formula）と呼ばれる。台形公式の幾何学的な意味を考えてみれば、何故台形公式と呼ばれるのか
容易に理解できる。
1. 数値解析の基礎
10
うえの公式における原点を順々にずらせば、式 (7) の積分を計算することができる。結局、台
形公式を用いた場合の式 (7) の積分公式は、
b
h
h
[f (a) + 2f (a + h) + f (a + 2h)] + [f (a + 2h) + 2f (a + 3h) + f (a + 4h)]
2
2
h
+ · · · + [f (b − 2h) + 2f (b − h) + f (b)]
2
h
= [f (a) + 2f (a + h) + 2f (a + 2h) + · · · + 2f (b − h) + f (b)]
(8)
2
f (x)dx =
a
で与えられる。
シンプソン公式（Simpson formula）
台形公式よりもさらに精度を上げるためには、f (x) に対する１階微分 f と２階微分 f の差
分式を用いることによって達成することでができる。１階微分 f と２階微分 f の差分式を用
いると、
f (x) ≈ f (0) + xf (0) +
f1 − 2f0 + f−1 2
x2 f1 − f−1
f (0) ≈ f0 +
x+
x
2!
2h
h2
となる。１階微分 f に対しては、台形公式の場合と異なり３点公式を用いていることに注意。そ
の結果、差分式は原点に対して対称になり、以下のように奇数次の項はキャンセルしてしまう。
h
−h
f (x)dx =
=
≈
h
h
h
x2
x3
dx f (0) + xf (0) + f (0) + f (0) + · · ·
2!
3!
−h
x2
dx f (0) + f (0) + · · ·
2!
−h
x2 f1 − 2f0 + f−1
dx f (0) +
2
h2
−h
この式を区間 [−h, h] において積分すると、以下のようなシンプソン公式を得る。
h
−h
f (x)dx =
h
(f−1 + 4f0 + f1 ) + O(h5 )
3
この公式は、台形公式と殆ど変わらないように見えるが、台形公式に比べて２次だけ正確な近
似になっている。
式 (7) の各積分をシンプソン公式を用いて計算すれば、以下のようなシンプソン公式の積分
公式が得られる。
b
f (x)dx =
a
h
[f (a) + 4f (a + h) + 2f (a + 2h) + 4f (a + 3h) + · · · + 4f (b − h) + f (b)] (9)
3
一般的注意 台形公式およびシンプソン公式に含まれる f の符号はすべて同符号であることか
ら、これらの公式は、数値計算上、極めて安定である。精度をさらにあげるためには、高次の
微分係数を用いれば良いが、この場合、公式に含まれる f の符号が正にも負にもなる。従って、
1. 数値解析の基礎
11
差分公式の次数を上げて精度をあげるよりも、むしろ、低次の差分公式を用いて、刻み幅 h を
小さくするほうが得策である。
1.2.2
特異性のある積分
これまでの積分公式は、被積分関数が積分範囲において正則であり、積分の発散等の問題が
起こらないことを暗に仮定していた。しかし、被積分関数に特異点が存在していても、実際に
は、有限の積分値を与えるものも少なくない。そこで、被積分関数に特異性がある積分のいく
つかの技術的方法について述べる。
積分変数の変換
以下のような積分を考える。
1
I1 =
0
g(x)(1 − x2 )−1/2 dx =
1
0
f (x)dx
ここで、g(x) は x ∈ [0, 1] で正則であるとする。f (x = 1) = ∞ となることから、これまでの積
分公式を使うことができない。そこで、t = (1 − x)1/2 として、変数変換をすると、
1
I1 =
0
g(1 − t2 )2(2 − t2 )−1/2 dt
となり、t ∈ [0, 1] において被積分関数は正則になり、これまでの積分公式を使うことができる。
積分区間の分割
つぎの積分
1
I2 =
0
f (x)dx
において、x = 0 近傍で被積分関数が f (x) → Cx−1/2 であったとする。このとき、f (x = 0) = ∞
となることから、やはり積分公式を使うことができない。そこで、以下のように積分区間を分
割する。
1
h
I2 =
0
f (x)dx +
f (x)dx
h
ここで、h は積分区間に比べて十分に小さいとする。第１項は、f (x) の x = 0 近傍での漸近形
から、以下のように評価することができる。
h
0
f (x)dx ≈
h
0
Cx−1/2 dx = 2Cx1/2
h
0
= 2Ch1/2
第２項は正則な被積分関数であるので、これまでの積分公式を用いて評価することができる。
例題
以下に具体例をあげて、特異性のある積分を実行してみよう。
1
I=
0
1
dxf (x) =
0
1
dx √
1 − x2
1. 数値解析の基礎
12
を考える。この厳密解は、
1
0
dx √
1
π
= [arcsin x]10 =
2
2
1−x
となる。
x = 1 で f (x) は発散するから、x = 1 近傍で積分範囲を分割する。
1
I=
0
1−h
dxf (x) =
1
dxf (x) +
0
1−h
dxf (x)
第１項は被積分関数が有限であるから、積分公式を用いて求めることができる。シンプソン公
式を用いれば、
1−h
0
dxf (x) =
h
[f (0) + 4f (h) + 2f (2h) + · · · + 4f (1 − 2h) + f (1 − h)]
3
第２項を評価するために、x = 1− 近傍での f (x) の漸近形を求める。
1
f (x) = √
2
1
1
+
1+x 1−x
1/2
1
1
→√ √
2 1−x
従って、t = 1 − x と変数変換して、
1
1−h
dxf (x) ≈
h
0
√
1
dt √ t−1/2 = 2h1/2
2
最終的に、求める積分は、
I=
√
h
[f (0) + 4f (h) + 2f (2h) + · · · + 4f (1 − 2h) + f (1 − h)] + 2h1/2
3
で与えられる。刻み幅 h が小さくなるにつれて、最後の被積分関数の発散に対応する項の寄与
が緩やかに小さくなることに注意せよ。
演習問題 1-2-1
台形公式およびシンプソン公式を用いて、次の定積分 I を求めるプログラムを作成し、評価
せよ。
1
I=
0
2
e−x dx
尚、格子の数 N は、N = 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048 とすること。
I の厳密解は、誤差関数 erf(x) を用いて以下のように与えられる。
√
π
erf(1) = 0.7468241328124270
I=
2
うえで求めた解の精度（誤差）について論ぜよ。
演習問題 1-2-2
1. 数値解析の基礎
13
つぎの定積分 I を求めよ。
1
I=
1
dxf (x) =
0
0
dxx−2/3 (1 − x)−1/3
なお、積分は以下のように分割して求めよ。
1
0
h
dxf (x) =
0
1−h
dxf (x) +
1
dxf (x) +
h
1−h
dxf (x)
正則な被積分関数に対してどのような積分公式を用いたか明記すること。
この積分は解析的に解くことができ、その解は以下のように与えられる。
1
0
2π
dxx−2/3 (1 − x)−1/3 = √
3
刻み幅 h を変化させることにより、数値解と厳密解との誤差について論ぜよ。
1. 数値解析の基礎
1.3
14
非線形方程式の解法
関数 f (x) が複雑な関数型をもつとき、一般に、方程式 f (x) = 0 の解析的な解を求めること
は不可能である。そこで、本章では、関数 f (x) が任意の x に対して評価が可能であると仮定
したうえで、方程式 f (x0 ) = 0 の満たす数値解 x0 を求めるのに用いられる最も基本的な方法に
ついて紹介する。
２分法（直接法）
方程式の数値解を求めるためのもっとも確実かつ単純な方法は、２分法と呼ばれるものであ
る。２分法のアルゴリズムは以下のように纏められる。
1. 解の最初の試行値が、正確な解 x0 よりも小さかったとする。この場合、次の試行値を少し
大きい方へ進める。
2. もし f の符号が変われば、次の歩幅を半分にして x を逆方向に進める。
3. 最初に設定した精度に到達するまでうえのステップを繰り返すことにより、数値解を求
める。
うえのアルゴリズムから明らかなように、歩幅が求める数値解の精度に対応する。最初の試
行値が大きすぎた場合、解が存在せずに求める事ができない場合もある。この場合は、ある程
度の探索の後に、探索の方向を（うえの場合は、x が小さい方向に）変えることにより、必ず
解を求めることができる。
２分法のアルゴリズムは単純で、しかも、かなり高い確率で解を求めることができることか
ら、最初に近似解を求める際に多用される。ただし、歩幅の大きさと解近傍での関数の変化の
激しさによってはいくつかの解を見のがしてしまうこともあり、注意が必要である。
ニュートン法（Newton-Raphson method）
２分法よりも解を求める効率性を高めたものにニュートン法と呼ばれるものがある。ただし、
関数 f の導関数 f (x) が任意の x について解析的に計算できることが必要である。
いま、求めたい方程式 f (x0 ) = 0 の解 x0 の近くで、関数 f が線形であると仮定する。そのう
えで、x = xi での導関数 f (xi ) が解析的に計算することができると仮定する。一方、差分近似
を用いれば、この導関数 f (xi ) は、
f (xi ) ≈
f (xi+1 ) − f (xi )
xi+1 − xi
と近似できる。常に、f (xi+1 ) = 0 と仮定すれば（図 2 参照）、即ち、導関数の差分式に用いら
れる x 座標の一方（ここでは、xi+1 ）は常に f (x) = 0 との交点であると仮定すれば、求める
ニュートン法による公式は以下のように与えられる。
xi+1 = xi −
f (xi )
f (xi )
(10)
1. 数値解析の基礎
15
具体的には、解の最初の試行（推定）値を x0 として、うえの公式を用いて x1 を求める。こ
れを順次くり返すことにより、速い収束性のもとで f (x) = 0 の数値解を求めることができる。
通常、計算機の平方根に対する組み込み演算は、このアルゴリズムを用いている。
f(x)
x0
xi+1
xi
図 2: ニュートン法の幾何学的な意味。
ただし、解近傍で導関数がゼロに近付く場合、つまり解近傍で関数 f (x) が極値をもつ場合
は、公式からわかるように発散してしまい、解を求めることが困難になる。
セカント法（Secant method）
ニュートン法では、導関数 f (x) が解析的に計算できることが前提であった。そこで、ニュー
トン法に含まれる導関数を、差分近似
f (xi ) ≈
f (xi ) − f (xi−1 )
xi − xi−1
でおきかえることにより、導関数に対する前提条件を避けたものをセカント（Secant）法とい
う。うえの差分近似をニュートン法に対する公式に用いれば、
xi+1 = xi − f (xi )
xi − xi−1
f (xi ) − f (xi−1 )
(11)
となる。ただし、この公式を用いるためには、２点 xi , xi−1 での関数値が必要になる。従って、
最初の試行値としては x0 および x1 を設定する必要がある。
一般に、最初の試行（推定）値の設定が真の解に十分近くない場合は、ニュートン法と同様
にセカント法は収束性がかなり悪い。導関数を２点間で線形近似しているのであるから、この
２点が十分に近い必要があるのは、容易に推測できる。つまり、この２点で関数が複雑な形状
（導関数が符号も含めて大きく変わる性質）をもっていれば、殆ど収束させるのは不可能であ
ろう。そこで、まず２分法を使って近似解を求め、つぎに、ニュートン法やセカント法によっ
て精度をあげていくという方法が一般的である。
1. 数値解析の基礎
16
演習問題 1-3-1
２分法、ニュートン法およびセカント法を用いて、次の方程式を解くプログラムを作成し、数
値解を求めよ。
f (x) = x5 + 3x − 5 = 0
尚、２分法およびニュートン法に対する初期値は x0 = 0.5、２分法における最初の刻み幅は
∆x = 0.25 とせよ。セカント法に対しては、初期値は x0 = 0.5, x1 = 1 とせよ。反復数 n ごと
に得られる数値解を書き出すこと。２分法により求まる解の精度は 10−6 として、それぞれの
方法の収束する早さを比べよ。
さらに、初期値をいろいろ変えることにより、それぞれの方法の収束性について論ぜよ。
f (x0 ) = 0 の有効数字 16 桁までの厳密解は以下のように与えられる。
x0 = 1.108575508057246
ちなみに、y = f (x) をプロットすると図 3 のようになる。
30
20
10
-1
-0.5
0.5
1
1.5
-10
図 3: y = x5 + 3x − 5 の曲線。
2
1. 数値解析の基礎
1.4
1.4.1
17
微分方程式の解法
常微分方程式の解法：陽解法
物理学における基本法則の多くは、微分方程式で定式化される場合が多い。従って、物理的
現象をモデル化するためには、微分方程式の数値解を求める作業が重要になる。本章では常微
分方程式の数値解を求める代表的な方法について学習する。
常微分方程式は、数個の１階の常微分方程式を連立させることにより解くことができる場合
が多い。例えば、質量 m の粒子に力 F (z) が働いている場合の１次元運動に対するニュートン
の運動方程式は、以下のような２階の常微分方程式で与えられる。
m
d2 z
= F (z)
dt2
いま、運動量 p(t) を
dz
dt
と定義すれば、ニュートンの運動方程式は次の２元の連立１階微分方程式に帰する。
p
dz
=
dt
m
dp
= F (z)
dt
つまり、多くの微分方程式が連立した１階の微分方程式に変形することができることから、
適当な初期条件（例えば、y(x = 0) = y0 ）のもとで以下のような１階微分方程式の数値解を求
めることができれば良い。
dy
= f (x, y)
dx
p(t) = m
基本的なアルゴリズム
一般に、１階微分方程式の数値解を求める基本的なアルゴリズムは簡単である。
初期条件 y(x = 0) = y0 に対して、以下の微分方程式
dy
= f (x, y)
(12)
dx
を数値的に求めることを考えよう。つまり、ある x（例えば、x = 1）における解 y(x) を求め
ることを考える。
1. 区間 [0, 1] を長さ h = 1/N の区間に等分割する。
2. 与えられた微分方程式を差分化して yn と yn+1 等を関係づける漸化式を作る。
3. x = 0 から１段階づつ積分する。
このときの漸化式に含まれる yn , yn+1 , yn+1 , . . . の数や yn+1 に対する次数によって、計算精
度が変わる。
オイラー法（Euler method）
最低次のオイラー法は、yn , yn+1 の漸化式に基づく。
1. 数値解析の基礎
18
式 (12) における左辺を、１章で導いた前進差分でおきかえると
yn+1 − yn
= f (xn , yn )
h
これを変形すれば、
yn+1 = yn + h f (xn , yn ) + O(h2 )
(13)
のような漸化式となる。これが、与えられた微分方程式に対する最低次数でのオイラー法によ
る漸化式である。与えられた初期条件から、うえの漸化式を用いて順々に yn+1 を求めること
により、求めたい y(x) での値を得ることができる。
漸化式によって積分する各段階での誤差は O(h2 ) であるが、区間全域にわたってこの誤差が
積算されることから、全体としての誤差は N O(h2 ) ≈ O(h) となる。従って、誤差を抑えるた
めには、刻み幅 h を小さくする必要がある（ただし、それに応じて計算時間は長くなる）。し
かし、計算精度が低いことから、積分区間が大きくなるとすぐにこの方法は破たんする。
テイラー級数を用いたオイラー法（Tylor method）
うえの方法の精度をあげるために、通常、高次のオイラー法を用いる。yn+1 を高次までテイ
ラー展開すれば、
yn+1
dy h2 d2 y = y(xn + h) = yn + h
+
+···
dx n
2 dx2 n
1
= yn + hf (xn , yn ) + h2 yn + O(h3 )
2
となる。ここで、与えられた微分方程式 yn = f (xn , yn ) を用いた。また、
yn =
∂f
∂f dyn
∂f
df (xn , yn )
∂f
=
+
=
+ f (xn , yn )
dx
∂x ∂y dx
∂x
∂y
であるから、テイラー級数を用いたオイラー法（テイラー法）は
1
∂f ∂f yn+1 = yn + hf (xn , yn ) + h2
+ f (xn , yn )
+ O(h3 )
2
∂x n
∂y n
(14)
となる。これは、先ほどのオイラー法 [式 (13)] に比べて１次だけ高精度な公式になっている。
与えられた微分方程式の右辺 f (xn , yn ) が簡単に微分できるときは、式 (14) は非常に有用で
ある。テイラー展開のさらに高次の項まで残せば、さらに精度をあげることができるが、実際
にはあまり使われない。
多段階法（Adams method）
微分方程式
dy
= f (x, y)
dx
の形式解は、
yn+1 = yn +
xn+1
xn
dxf (x, y)
1. 数値解析の基礎
19
で与えられる。この第２項の積分を xn における値 hf (xn , yn ) で近似したものがオイラー法で
あった。多段階法では、同様に第２項の積分を xn における値で近似をするが、その際、被積分
関数をいくつかの x 点での関数 f (x, y) から補間する。
二段階法（Two-step method）
二段階法では、関数 f (x, y) を x ∈ [xn−1 , xn ] での値を用いて線形補間する。つまり、x =
xn−1 , xn での関数 f (x, y) をそれぞれ fn−1 , fn として、x ∈ [xn , xn+1 ] での関数 f (x, y) を以下
のように１次関数で近似する。
fn − fn−1
(x − xn−1 ) + fn+1
xn − xn−1
x − xn−1
x − xn
=
fn −
fn−1
h
h
従って、二段階法での漸化式は以下のように与えられる。
xn+1
x − xn
x − xn−1
fn −
fn−1
yn+1 = yn +
dx
h
h
xn
1
3
= yn + h fn − fn−1
2
2
f (x, y) =
(15)
三段階法（Three-step method）
三段階法では、さらに１つ前の xn−2 での fn−2 を用いて、x ∈ [xn , xn+1 ] での関数 f (x, y) を
２次関数で近似する。
x ∈ [xn , xn+1 ] での関数 f (x, y) をまえと同様に以下のように仮定する。
f (x, y) = c0 + c1 + +c2 x2
この関数が fn , fn−1 , fn−2 を通るという条件から、定数 c0 , c1 , c2 をきめる決めることができる。
その結果、
−3fn−2 + 4fn−1 − fn
fn−2 − 2fn−1 + fn
(x − xn−2 ) +
(x − xn−2 )2
2h
2h2
となる。従って、三段階法での漸化式は以下のように計算することができる。
f (x, y) = fn−2 +
yn+1 = yn +
xn+1
dxf (x, y)
(16)
xn
演習問題 1-4-1
三段階法における上式（式 16）の積分を実行することにより、公式を導け。公式は、yn , h, fn , fn−1 , fn−2
のみを使って表すこと。
陰解法（Implicit method）
オイラー法では、既に分かっている yn を用いて、漸化式から順次 yn+1 を求めるという方法
であり、このような方法を陽解法（Explicit method）という。一方、本セクションで紹介する
陰解法は、未知の yn+1 を含んだ漸化式を用いて微分方程式を解く方法である。
1. 数値解析の基礎
20
いま、２つの格子点 xn , xx+1 の中間点を xx+1/2 とする。微分方程式 (12) における左辺を、
点 xx+1/2 を含めた微分の３点公式を用いて置き換えると、
yn+1 − yn
h
= f (xn+1/2 , yn+1/2 )
≈
1
[f (xn , yn ) + f (xn+1 , yn+1 )]
2
となる。つまり、陰解法での漸化式は以下のようになる。
yn+1 = yn +
h
[f (xn , yn ) + f (xn+1 , yn+1 )] + O(h3 )
2
(17)
勿論、yn+1 は未知であるから、この漸化式を評価するためにはニュートン法等を用いる必要が
あり、オイラー法に比べて１ステップ余計なプログラミングが必要になる。
しかし、f (x, y) = g(x)y のように、微分方程式 (12) の右辺が y について線形である場合は、
うえの漸化式は、以下のように簡単になる。
yn+1 = yn
1 + h2 g(xn )
1 − h2 g(xn+1 )
(18)
この場合、右辺はすべて簡単に求めることができるから、オイラー法と殆ど変わらない手間で
高精度な解析が可能になる。
演習問題 1-4-2
境界条件 y(0) = 1 のもとで、微分方程式
dy
= −xy
dx
を、オイラー法、テイラー法、三段階法、陰解法を用いて解くプログラムを作成せよ。そして、
y(1) および y(3) を求めよ。刻み幅としては、h = 0.5, 0.2, 0.1, 0.05, 0.02, 0.01, 0.005, 0.002,
0.001 として、各方法の精度を調べよ。
うえの微分方程式の厳密解は
2
y = e−x /2
で与えられる。
ルンゲクッタ法（Runge-Kutta method）
微分方程式を解くアルゴリズムには、うえで説明したオイラー法や陰解法の他に数多く存在
する。しかし、もっとも一般的なものはルンゲ・クッタ法と呼ばれるものであろう。ルンゲ・
クッタ法には、さまざまな次数（精度）のものがあるが、ここではもっとも基本的な２次およ
び高次のルンゲ・クッタ法について説明する。
２次のルンゲ・クッタ法
まず、２次のルンゲ・クッタ法を導くために、微分方程式
dy
= f (x, y)
dx
1. 数値解析の基礎
21
を１段階だけ積分したものを考えよう。
yn+1 = yn +
xn+1
f (x, y)dx
xn
積分に含まれる f (x, y) を積分区間の中央で展開すると、
yn+1 = yn + hf (xn+1/2 , yn+1/2 ) + O(h3 )
ここで、誤差はテイラー展開したときの２次の項から生じることに注意せよ。１次の項は奇関
数であることから消える。これを評価するためには、yn+1/2 が必要になる。そこで、yn+1/2 を
展開して、
yn+1/2 = y(xn +
h
h dy ) = = yn +
+ ···
2
2 dx n
h
≈ yn + f (xn , yn )
2
k
= yn +
2
ここで、k = hf (xn , yn ) である。従って、２次のルンゲ・クッタ法は以下のような２段階の評
価を繰り返すことにより、求めたい解 y(x) を得ることができる。
まず、yn+1/2 を評価するために、
k = hf (xn , yn )
(19)
を求める。次に、この yn+1/2 と xn+1/2 を用いて、以下のように yn+1 を求める。
yn+1 = yn + hf (xn + h/2, yn + k/2) + O(h3 )
(20)
この方法は、テーラー級数を用いたオイラー法と同程度の精度をもつが、（オイラー法で必要
となった）関数 f が解析的に微分可能であるといった条件がいらない。そのかわり、各格子点
において、f を２回評価する手間が必要になる。
高次のルンゲ・クッタ法
高次のルンゲ・クッタ法は、積分の部分をシンプソン公式等の高精度なものにおきかえるこ
とにより導くことができる。
積分の部分をシンプソン公式を用いて書き換えると、
yn+1 = yn +
h
f (xn , yn ) + 4f (xn+1/2 , yn+1/2 ) + f (xn+1 , yn+1 ) + O(h5 )
6
となる。そこで、３次のルンゲ・クッタ法は以下の３段階の評価からなる。
k1 = hf (xn , yn )
k2 = hf (xn + h/2, yn + k1 /2)
1. 数値解析の基礎
22
k3 = hf (xn + h, yn − k1 + 2k2 )
1
(k1 + 4k2 + k3 ) + O(h4 )
6
４次のルンゲ・クッタ法では、シンプソン公式のなかの 4f (xn+1/2 , yn+1/2 ) を２項に分ける
ことにより、以下の４段階の評価からなる。
yn+1 = yn +
k1 = hf (xn , yn )
k2 = hf (xn + h/2, yn + k1 /2)
k3 = hf (xn + h/2, yn + k2 /2)
k4 = hf (xn + h, yn + k3 )
1
(k1 + 2k2 + 2k3 + k4 ) + O(h5 )
6
４次のルンゲ・クッタ法は安定性が非常に高い方法と考えられており、ルンゲ・クッタ法のな
かでも最も一般的である。
yn+1 = yn +
演習問題 1-4-3
つぎの微分方程式
dy
= −xy,
y(0) = 1
dx
を、２次、３次および４次のルンゲ・クッタ法を用いて解くプログラムを作成せよ。そして、
y(1) および y(3) を求めよ。刻み幅として、h = 0.5, 0.2, 0.1, 0.05, 0.02, 0.01, 0.005, 0.002,
0.001 として、厳密解と比較することにより、精度について調べよ。
厳密解は、
2
y = e−x /2
である。
1. 数値解析の基礎
23
２階の線形微分方程式
本章では、物理の重要な微分方程式として頻繁に顔を出す２階の線形微分方程式について考
察する。ここで考察する微分方程式は以下のようなものである。
d2 y
+ k(x)2 y = S(x)
dx2
(21)
このようなかたちで表されるものとしては、静電ポテンシャルに対するポアソン方程式
∇2 Φ = −4πρ
がある。（たびたびそうであるように）ポテンシャル Φ と電荷密度 ρ が球対称（つまり r のみ
の関数）であるならば、ポアソン方程式は、
1 d 2 dΦ
r
= −4πρ(r)
r2 dr
dr
と書くことができる。そこで、
φ(r)
r
とすれば、ポアソン方程式は、式 (21) と同様な２階の線形微分方程式になる。
Φ(r) =
d2 φ(r)
= −4πrρ(r)
dr2
(22)
また、中心力ポテンシャル V (r) のもとで運動する粒子に対する Schrödinger 方程式は、粒子
の波動関数を
R(r)
Yl,m (θ, ϕ)
Ψ(r) =
r
のように表したとすると、
d2 R(r)
+ k(r)2 R(r) = 0
(23)
dr2
のようになり、やはり、式 (21) と同様な２階の線形微分方程式になる。ここで、Yl,m (θ, ϕ) は
球面調和関数、E は粒子のエネルギー、m は粒子の質量である。k 2 (r) は、
2m
l(l + 1)
k(r) = 2 E −
− V (r)
2mr2
h̄
2
で与えられる。
これらの微分方程式は、前章の微分方程式と基本的に同じものである。前章では初期条件の
み与えられた条件のもとで微分方程式を解いたが、本章では離れた２点での変数の値が拘束条
件として与えられているところが異なっている。このような問題を（２点）境界値問題という。
ヌメロフの方法（Numerov method）
２階の線形微分方程式
d2 y
= −k(x)2 y + S(x)
dx2
1. 数値解析の基礎
24
を解くもっとも簡単で効率的な方法としてヌメロフの方法と呼ばれるものがある。この方法で
は、左辺の２階微分を３点差分公式でおきかえて近似する。つまり、
yn±1 = yn ± hyn +
h2 h3 h4 (4)
y ± yn + yn + O(h5 )
2 n
6
24
を用いて、２階微分の３点差分公式を書くと、
h2 (4)
yn+1 − 2yn + yn−1
=
y
+
y + O(h4 )
n
h2
12 n
となる。従って、
yn+1 − 2yn + yn−1 h2 (4)
d2 y
=
− yn + O(h4 )
dx2
h2
12
yn の４階微分は、もとの微分方程式を用いることにより、
yn(4) =
d2
2
(−k
y
+
S)
2
dx
x=x
n
(k 2 y)n+1 − 2(k 2 y)n + (k 2 y)n−1 Sn+1 − 2Sn + Sn−1
≈ −
+
+ O(h2 )
h2
h2
となるから、最終的に、微分方程式 (21) は以下のように変形される。
h2 2
5h2 2
k yn +
1 + kn+1
yn+1 − 2 1 −
12
12 n
=
h2 2
1 + kn−1
yn−1
12
h2
(Sn+1 + 10Sn + Sn−1 ) + O(h6 ) (24)
12
微分方程式が与えられていれば、各格子点での kn2 と Sn は解析的に簡単に評価できるから、
ヌメロフの方法は簡単で効率の良い方法である。式 (24) は、yn+1 と yn−1 に関する線形方程式
になっているから、初期条件として、y0 , y1 が与えられれば、順々に yn を評価することができ
る。この公式は、４次のルンゲクッタ法よりも精度が良い。
演習問題 1-4-4
次のような境界条件のもとで微分方程式
d2 y
= −4πy,
dx2
y(0) = 1,
y (0) = 0
を、ヌメロフの方法を用いて解いて、y(1), y(3) および y(5) を求めよ。最初の漸化式を評価す
るためには、y1 を知る必要があるが、y1 = y(x = h) を x = 0 のまわりでテイラー展開して、与
えられている境界条件を用いることにより y1 を求めよ。刻み幅としては、h = 0.1, 0.01, 0.001,
0.0001 とせよ。
厳密解
√
y = cos(2 πx)
と比較して、得られた解の精度について調べよ。
1. 数値解析の基礎
1.4.2
25
数値解の収束性と誤差
１階微分方程式：解の収束性
つぎのような微分方程式を単純なオイラー法を用いて解くことを考えよう。
dy
= −y
dx
ここで、初期条件は y(x = 0) = 1 とする。明らかに、厳密解は y = e−x で与えられる。うえの
式の左辺を３点公式を用いて差分化すると、
yn+1 − yn−1
= −yn
2h
となるから、次のような漸化式を解くことになる。
yn+1 = yn−1 − 2hyn
(25)
ここで、初期条件は以下のようなものを用いればよい。
y0 = y(x = 0) = 1
dy(0) h2 d2 y(0)
h2
+
+
.
.
.
≈
1
−
h
+
dx
2 dx2
2
実際にうえの初期条件のもとで漸化式を評価すると、大きい x で y(x) が振動してその絶対
値は発散することがわかる。つまり、十分大きいで x で数値解は収束しない。
そこで、漸化式 (25) を解析的に解いてみよう。いま、yn = Arn（A, n は定数）とすると、特
性方程式は以下のようになる。
y1 = y(h) = y(0) + h
Arn+1 = Arn−1 − 2hArn
これを解くと、
1
r = −h ± (1 + h2 ) 2
そこで、
1
1
r+ = −h + (1 + h2 ) 2 ≈ 1 − h,
r− = −h − (1 + h2 ) 2 ≈ −(1 + h)
とおけば、漸化式 (25) の一般解は、
n
n
yn = yn+ + yn− = Ar+
+ Br−
≈ A(1 − h)n + B(−1)n (1 + h)n
と書ける。
十分大きな x = nh では、
+
y (x) =
yn+
x
=A 1−
n
y − (x) = yn− = B(−1)n 1 +
n
x
n
n
→ Ae−x
→ (−1)n Bex
1. 数値解析の基礎
26
となる。１番目の解が与えられた微分方程式の一般解に対応する。また、２番目の解 y − (x) は
発散解であり、B = 0 である限り、常に発散する。計算機上では確実に B = 0 とすることは困
難であり、計算機イプシロンにより、B = 0 の解が必ず含まれてしまう。その結果、大きい x
では数値解が発散してしまい、解が不安定になる。
演習問題 1-4-5
つぎの微分方程式
dy
= −y
dx
を、初期条件 y(x = 0) = 1 のもとでオイラー法を用いて解くプログラムを作成せよ。差分化し
た漸化式を解くときの初期条件は本文の記述を参考にせよ。
十分に大きい x で、数値解が発散することを確かめよ。
２階微分方程式：自明ではない境界条件と数値誤差
境界値問題を解く際、その境界条件が解析的に解くときは自明であっても、数値的に解く場
合は自明ででない場合も多い。そのような例について説明しよう。
つぎのようなポアソン方程式を極座標系で考える。
1 d
dΦ
∇ Φ = 2
r2
r dr
dr
= −4πρ(r)
2
=
1 d2
(rΦ)
r dr2
ここで、Φ(r) = ϕ(r)/r とすれば、
d2 ϕ(r)
= −4πrρ(r)
dr2
となり、ヌメロフ法を使うことができる。電荷密度 ρ(r) が以下のように与えられたとする。
ρ(r) =
1 −r
e
8π
このとき全空間での総電荷 Q は、
Q=
d3 rρ(r) =
∞
0
dr4πr2
1 −r
e =1
8π
となる。そして、うえのポアソン方程式はこの電荷分布に対して解析に解くことができ、その
解は、
1
ϕ(r) = 1 − (r + 2) e−r
2
で与えられる。ここで、境界条件としては、
ϕ(r = 0) = 0,
ϕ(r = ∞) = 1
とした。これらの物理的意味について考えてみよ。
1. 数値解析の基礎
27
さて、ヌメロフ法を使ってこのポアソン方程式を漸化式に書きかえると、
ϕn+1 − 2ϕn + ϕn−1 =
h2
(Sn+1 + 10Sn + Sn−1 )
12
ここで、
1
Sn = −4πrn ρ(rn ) = − rn e−rn
2
従って、ϕ0 と ϕ1 が必要になる。解析解を求める際に用いた ϕ∞ は、この場合、自明ではない
ことに注意せよ。うえの境界条件の一つは、
ϕ0 = ϕ(r = 0) = 0
ところが、ϕ1 は解が求まらなければわからない。これが、自明ではない境界条件の例である。
そこで、物理的な考察から ϕ1 を求める。与えられた電荷分布のもとでクーロンの法則を用
いると
ρ(r)
Φ(r ) = d3 r |r − r|
r = 0 として、変数変換をすれば、
Φ(0) =
ρ(r)
= 4π
d r
r
3
∞
0
∞
drrρ(r) =
0
dr
1
re−r
=
2
2
つまり、r ≈ 0 近傍で、Φ(r ≈ 0) = 1/2 となる。従って、
ϕ1 ≈ +hΦ(0) ≈ +
h
2
となる。
うえで与えられたポアソン方程式においては、実際に数値解を求めるとその解の誤差が線形
で増大する。これは、計算機に含まれるまるめ誤差により、ϕ0 , ϕ1 の境界条件に必ず誤差が含
まれるためである。そして、この誤差が n が増大すると伴に線形に増大する。これは、ポアソ
ン方程式に特有の性質である。
例えば、
d2 ϕ(r)
= −4πrρ(r)
dr2
の一般解は、右辺をゼロとおいて解いたものであり、
ϕ(r) = c0 + c1 r
となる。つまり、このポアソン方程式における解は、
ϕ(r) = [c0 + c1 r] + ϕsp (r)
(26)
ここで、ϕs p(r) は、うえで求めた右辺がゼロでないときの特殊解である。勿論、境界条件から
c0 = c1 = 0 でなくてはならない。しかしながら、いわゆる計算機イプシロンによって、数値解
1. 数値解析の基礎
28
は常に式 (26) で与えられる。その結果、厳密解 ϕsp (r) のうえに線形で増加する誤差項 c0 + c1 r
が上乗せされる。
しかしながら、このときの誤差は線形で増加することがわかっているから、誤差を取り除く
のは簡単である。最初の数ステップ（例えば、10 ステップ）から誤差項の定数 c1 を以下のよ
うに見積もることができる。
ϕ11 − ϕ1
c1 =
10h
従って、正しい数値解は
ϕn → ϕn − c1 · (nh)
とすれば良い。
あるいは以下のように、誤差を取り除くこともできる。この問題の場合、r → ∞ の境界条件
Φ(r) → 1/r から、ϕ(r = ∞) = 1 となる。従って、十分に大きい r （つまり n）で
ϕ(r = ∞) = ϕn = 1
として、n が減る方向にヌメロフの漸化式を用いれば良い。そうすれば、線形の誤差項は n が
減少するにつれて小さくなるから問題はない。
演習問題 1-4-6
電荷分布
1 −r
e
8π
に対して、つぎのポアソン方程式を、ヌメロフの方法を用いて r = 0 から r = 50 まで解け。
ρ(r) =
d2 ϕ(r)
= −4πrρ(r)
dr2
ここで、境界条件は、
ϕ0 = 0,
ϕ1 =
h
2
とし、刻み幅は h = 0.1 とせよ。
厳密解は本文で与えられているように、
1
ϕ(r) = 1 − (r + 2)e−r
2
である。
境界条件に、5%の誤差を含ませたとき、
ϕ0 = 0,
ϕ1 =
数値解の誤差が線形に増大することを確かめよ。
h
· 0.95
2
1. 数値解析の基礎
1.5
29
物理への応用
本章では、これまで学習した数値計算の手法を具体的な物理の問題に適用することによって、
実際に、数値計算がどのように応用されるかを体験する。
まず、研究課題についての物理的な背景と意味を説明したうえで、プログラミングを行うた
めのアルゴリズムを詳細に示す。そのうえで、各自が作成したプログラムを用いて、いくつか
の興味深い問題に適用し、研究課題の物理的内容のより深い理解を試みる。
1.5.1
物理的背景：２原子分子の振動運動
本年度の研究課題は、量子力学の創設期において大きな役割を果たしたボーアの半古典的量
子化条件を用いて、解析的に解くことが不可能なポテンシャルのもとでの量子化されたエネル
ギー準位を求めるものとする。
具体的な問題設定としては、２個の原子核が電子によって結合された H2 のような２原子分
子の振動運動を考察する。２原子分子の振動を引き起こすポテンシャル場は、原子核および電
子の存在する場所によって主にクーロン力によって確定される。しかし、電子は原子核に比べ
てずっと軽いために原子核の変位に対して即座に応答することができる（このような仮定をボ
ルン近似あるいは断熱近似という）。従って、ある静的なポテンシャル場のもとで、原子核の相
対距離の振動が起きる（つまり、電子の運動はポテンシャルのなかに既に繰り込まれていると
考える）。
一般に、このような状況のもとでのポテンシャル V (r) は、原子核の相対距離が非常に小さ
いときは原子核間のクーロン力とパウリの排他律による強い斥力が働き、長距離においては弱
い引力になる。このようなポテンシャルとして良く用いられるものとしては、以下のようなレ
ナード・ジョーンズ型あるいは 6-12 型ポテンシャルと呼ばれるものがある。
12
V (r) = 4V0
a
r
−
6 a
r
図 4 に示してあるように、このポテンシャルは、rmin = 21/6 a ≈ 1.12a において極小値をもつ。
いま（相対）振動のエネルギー E を −V0 < E < 0 とすれば、原子核は距離 r が rin < r < rout
の間の相対的な振動運動となる（図 4 参照）。つまり、原子核の相対距離の運動は、内側と外側
の転回点 rmin , rmax の間を振動する。ただし、ポテンシャルは放物線型ではないため、調和振
動とは限らないことに注意せよ。エネルギーの保存則より、相対運動の運動エネルギーとポテ
ンシャル・エネルギーの和
p2
+ V (r)
E=
2m
は、一定に保たれる。その結果、位相空間（座標 r と運動量 p の空間）で相対運動は閉じた軌
道を描く。この軌道の式は、
p = ± 2m(E − V (r))
で与えられる。
1. 数値解析の基礎
30
V(r)/Vo
1
0.5
rmin
rin
1
1.2
rout
1.4
1.6
1.8
2
r/a
-0.5
-1
図 4: レナード・ジョーンズ型ポテンシャル
この軌道に対して、ボーアの量子化条件を用いる。いま、与えられたエネルギー E に対する
作用 S を考える。
S(E) = p(r)dr
作用 S は位相空間における閉じた軌道の面積に等しい。ボーアの量子化条件によれば、この作
用Sが
1
S = (n + )h
2
の条件を満たす場合にのみ運動が許される。ここで、h はプランク定数であり、n は自然数、
n = 1, 2, . . .。上式の 1/2 は、通常、ボーアの量子化条件において省かれるが、量子力学による
結果（不確定性原理から生じる）と一致させるためにここではつけ加えている。その結果、運
動の許されるエネルギーはとびとびの離散化されたエネルギー En になり、つぎの式で決める
ことができる。
rout 1
S(En ) = 2
2m(En − V (r))dr = n +
h
(27)
2
rin
これは積分方程式であり、積分に対しては台形公式やシンプソン公式等、解 En を求めるため
には２分法やニュートン法等を用いて数値的に解くことができる。これが、ボーアの量子化条
件を用いて離散化されたエネルギー En を求める方法である。
参考までに、量子力学的には、この相対距離の振動運動は以下のような１次元 Schrödinger
方程式で表される。
h̄2 d2
−
+ V (r) ψn = En ψn
2m dr2
ここで、m は２個の原子核の換算質量である（m = M/2）。振動エネルギー E が −V0 < E < 0 の
ときは、運動は束縛状態の波動関数 ψn で規定され、エネルギーは離散的になる。この Schrödinger
方程式は、ヌメロフの方法等を用いて解くことができる。
1. 数値解析の基礎
1.5.2
31
数値解析の基本アルゴリズム
それでは、式 (27) を数値的に解くアルゴリズムについて具体的に考察しよう。まず、式 (27)
に含まれるポテンシャルとしてレナード・ジョーンズ型を仮定する。計算を簡単にするために
すべての量を無次元化する。そこで、以下のように無次元化したエネルギー、距離等を定義す
れば、
1/2
E
r
2ma2 V0
ε=
, x= , γ=
V0
a
h̄2
無次元化したレナード・ジョーンズ・ポテンシャル v(x) は、
1
1
v(x) = 4
−
x12 x6
と表される。図 4 は、実は、このポテンシャルをプロットしたものである。
量子化されたエネルギー εn を求める式 (27) は、
s(εn ) = 2γ
xout 1
εn − v(x)dx = n +
2π
2
xin
(28)
となる。パラメータ γ は、系を特徴付けている物質パラメータに対応する。γ が大きい場合、
つまり h̄ が小さい場合は、古典系に近い系を表す。γ に含まれる a, V0 は、それぞれ、分子の回
転エネルギーに対する慣性モーメントと分子をばらばらにするのに必要な解離エネルギーから
実験的に求めることができる。H2 分子は γ = 21.7 である。また、重い原子核からなる酸素分
子 O2 分子の場合は γ = 150 である。
— ステップ１ —
無次元化されたレナード・ジョーンズ・ポテンシャル v(x) において、あるエネルギー（εn = −0.75
とする）に対する転回点 xin , xout （図 4 参照）を２分法により求めるプログラムを作成せよ。
無次元化したエネルギーは εn = −0.75、精度は δx = 10−8 、座標の刻み幅 δx = xmin /10、
転回点を探すための x の初期値としては極小点 x = xmin とせよ。
———————
— ステップ２ —
うえで求めた転回点 xin , xout およびシンプソン公式を用いて、無次元化した作用 s(εn ) を求め
るための式 (28) の積分を計算するプログラムを作成し、ステップ１で作成したプログラムと結
合せよ。
物質パラメータは γ = 100 とし、積分を行うための分割数は N = 100 とせよ。
———————
この時点において、各自の作成したプログラムは、無次元化したエネルギー εn を与えれば、
式 (28) の作用 s(εn ) を評価することができるようになっているはずである。この作用 s(εn ) を
用いれば、式 (28) からエネルギー εn の量子準位 n を求めることができる。即ち、
n=
s(εn )
− 0.5
2π
1. 数値解析の基礎
32
により、量子数 n が求まる。
— ステップ３ —
物質パラメータを γ = 100 として、エネルギーが εn = −0.0005 に対して、エネルギー準位 n
を求めよ。
これは、振動エネルギーが連続的になる領域までの量子化されたエネルギー準位の最大数（た
だし、この数値計算で仮定されている精度、つまり、δx = 0.0005 の範囲以内で求めることが
できる準位数）に対応する。
———————
さて、ステップ３で求めた量子準位の最大数 n までのエネルギー準位 εn をセカント法を用
いて順々に求めよう。
そのためには、未知数 εn に対する以下の方程式
f (εn ) = s(εn ) − 2π(n + 0.5) = 0
を解いて、エネルギー準位 εn をセカント法により求めればよい。
例えば、基底状態 n = 1 のエネルギー準位は、方程式
f (ε1 ) = s(ε1 ) − 2π(1 + 0.5) = 0
の解であることから、エネルギー準位 ε1 が決まる。同様に、n = 2, 3, . . . , n に対するエネル
ギー準位 εn を求めることができる。
作成したプログラムが正しく動作するかどうか検証するためには、解析的に解くことができ
るポテンシャルが解析できるように、作成したプログラムを修正し、数値解と解析解を比較す
るということが良く行われる。例えば、レナード・ジョーンズ・ポテンシャルをポテンシャル
の極小点（xmin ）の近傍でテイラー展開し、放物線型ポテンシャルで近似するということも良
く行われる。また、この研究課題では、半古典的手法により量子化された束縛エネルギー準位
を求めた。勿論、方程式を射撃法等の方法を用いて、直接、解くことも可能である。力のある
者は挑戦してみると良い。
問題 S-1
セカント法を用いて、最大量子準位数 n までのエネルギー準位 εn を求めるプログラムを完成
させ、それぞれの量子準位におけるエネルギー準位を求めよ。尚、それぞれの準位における転
回点（xmin , xmax ）、エネルギーが収束するためのセカント法での反復数も同時に示すこと。
物質パラメータをいろいろと変化させることにより、量子効果の重要性を認識することがで
きる。γ = 20 − 200 において、エネルギー準位の変化を調べ、量子効果の重要性と関連づけて
議論せよ。
ヒント
1. 数値解析の基礎
33
基底状態のエネルギー ε1 を求めるためのエネルギーの初期値は
ε11 = −1 + 0.001,
ε21 = ε11 + 0.001
とし、そのときの関数値は
f (ε11 ) = 0 − 2π0.5 = −π
とせよ。また、セカント法でのエネルギーの刻み幅は、δε = 0.001 とせよ。
つぎのエネルギー ε2 を求めるには、初期値を
ε12 = ε1 ,
ε21 = ε12 + δε
としたうえで、
f (ε12 ) = f (ε1 ) − 2π1
とすれば良い。
1. 数値解析の基礎
1.6
34
モンテカルロ法の基礎
近年、計算機の高速化に伴って、大規模な数値シミュレーションがさまざま科学工学分野で行
われている。なかでも、モンテカルロ法（Monte Carlo method）と分子動力学法（Molecular
Dynamics method）は、大規模シミュレーションで使われる手法のなかでも、従来の差分方程
式を用いたものと伴に双璧をなすと言っても過言ではない。そこで、このセクションでは、モ
ンテカルロ法（Monte Carlo method）のごく基礎的な考え方と手法について説明する。
モンテカルロ法は、一般に、多重積分においてその威力を発揮する。例えば、統計力学など
では系の自由度に応じた多重積分の実行が必要になる。例えば、自由度が 10 という非常に小さ
な系について、30 重積分を台形公式等を用いて実行することは殆ど不可能である。各次元に対
して、積分区間を 10 点に差分化して台形公式あるいはシンプソン公式を用いた場合、全部で
1030 個の点について計算を行う必要がある。1 秒間に 108 回の演算が可能な計算機を用いたと
しても、計算時間は 1022 秒もかかる。宇宙の年齢が約 1020 秒程度であることを考えれば、こ
の計算がどれほど困難なものであるかは容易に想像できるであろう。しかし、モンテカルロ法
を用いることにより、簡単に解を見つけることができる場合も多い。
1.6.1
モンテカルロ法の基礎概念
モンテカルロ法の基本的な考え方は非常に単純である。つぎのような１次元積分をモンテカ
ルロ法で求めるには、N 個の乱数 xi ∈ [0, 1] を発生させ、そのときの被積分関数の値 f (xi ) か
ら平均値を求める。
1
N
1 I=
dxf (x) ≈
f (xi ) ≡ f (x)
(29)
N i=1
0
後述するように、この平均値 f (x)
は、N が大きくなるにつれて、求めたい積分値 I に近づく。
積分区間が任意の場合
積分区間が任意の場合は、積分変数 x を適当に変数変換してやればよい。x = (x − a)/(b − a)
とすれば、
b
I=
a
dxf (x) = (b − a)
≈ (b − a)
1
0
dx f (x )
N
1 f (xi )
N i=1
≡ (b − a) f (x )
ここで、新たなランダム変数 x は、区間 [a, b] でランダムに選ばれることに注意せよ。
多重積分の場合
多重積分に対するモンテカルロ法も１次元の場合と全く同様である。
1. 数値解析の基礎
35
つぎのような２重積分の場合は、２つの変数 x1 , x2 をランダムに選び、そのときの被積分関
数値から平均を求めて積分値とする。
1
I=
1
dx1
0
0
dx2 f (x1 , x2 ) ≈
N
1 f (xi1 , xi2 ) ≡ f (x1 , x2 )
N i=1
モンテカルロ法における誤差
モンテカルロ法の積分値の誤差は、統計論の中心極限定理から以下のように考察できる。
積分値を N 倍したものを ξ とすれば、fi ≡ f (xi ) をランダム変数とし、
ξ =N ·I =
N
fi = f1 + f2 + · · · + fN
i=1
従って、ξ の平均値は
ξ
= N I
= f1 + f2 + · · · + fN = N f 中心極限定理によれば、ξ の分散 σξ は
σξ2 = ∆ξ 2 = (∆f1 + · · · + ∆fN )2 = N ∆f 2 ≡ N σf2
となる。一方、N は固定された定数であるから、
σξ2 = N 2 σI2
従って、積分値 I の誤差は、
σξ
1 σf
σI
=
=√
I
ξ
N f つまり、積分値 I の誤差は、
√
1. 乱数の発生回数 N に対して、1/ N に比例して減少する。
2. 被積分関数 f (x) の標準偏差（変化の大きさ）に比例する。
従って、精度を高めるためには、試行回数（乱数のサンプル数）を増やすという当然のことと
もに、被積分関数の変化をできるだけ少なくする必要がある。この後者を行うために用いられ
る方法が、重み関数を用いたモンテカルロ法（importance sampling）である。
演習問題 1-6-1
つぎの多重積分を単純なモンテカルロ法により求めよ。
I =
0
π
8
dx1
0
π
8
dx2
0
π
8
dx3
0
π
8
dx4
π
8
0
dx5
0
π
8
dx6
× sin (x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 )
サンプル数は、N = 101 , 102 , 103 , 104 , 105 とせよ。
0
π
8
dx7
0
π
8
dx8
1. 数値解析の基礎
36
この積分の厳密解は
π
3π
π
+ 56 sin − 112 sin
8
4
8
≈ 5.371873 × 10−4 Iexact = 70 − 16 sin
である。
2 = I 2 − I
2 を求めて、log σ と log N をプロットし、
それぞれの N について分散 σN
N
1
σN ∝ √
N
であることを示せ。
1.6.2
重み関数を用いたモンテカルロ法
重み関数を用いたモンテカルロ法の基本的な考え方は、積分変数を適当に変換して、被積分
関数の変化をできるだけ小さく抑えるということにある。そこで、重み関数 w(x) を用いて式
（29）を評価することを考えよう。ここで、重い関数は、つぎのように規格化されているとする。
1
0
dxw(x) = 1
式（29）をつぎのように変形する。
1
I=
0
1
dxf (x) =
0
dxw(x)
f (x)
w(x)
dy/dx = w(x) として変数変換を行うと、
x
y(x) =
0
dxw(x)
となり、w(x) の規格化条件から、
y(x = 0) = 0, y(x = 1) = 1
である。従って、
1
I=
0
N
1 f (x(y))
f (x(yi ))
f (x(y))
≈
=
dy
w(x(y))
N i=1 w(x(yi ))
w(x(y))
つまり、ランダム変数を x から y に変換し、乱数により y ∈ [0, 1] を選んだうえで、x(y) を評
価し、f (x(y))/w(x(y))
を求めればよい。積分区間で被積分関数 f (x(y))/w(x(y))
が滑らか
に変化するように重み関数 w(x) を選ぶことにより、効率良く積分値を求めることができる。
重み関数を用いたモンテカルロ法とは、一様ランダムなもともとの変数 x 用いる方法（単純
なモンテカルロ法）から、ある特定な分布（重み関数）をもった乱数を用いて評価することを
意味する。
1. 数値解析の基礎
37
演習問題 1-6-2
つぎの積分を、単純なモンテカルロ法と重み関数を用いたモンテカルロ法により求めよ。
1
I=
0
π
dx
=
1 + x2
4
重み関数として、
1
w(x) = (4 − 2x)
3
とせよ。サンプル数を N = 10, 50, 100, 500, 1000, 5000, 10000, 50000 として、log σN と log N
をプロットせよ。ここで、
2
σN
= (IN − IN )2
である。
うえの重み関数が何故うえの計算において有効なのか説明せよ。
特定の分布をもつ乱数の生成
重み関数を用いたモンテカルロ法の鍵は、適当な重み関数に応じた乱数を生成することにあ
る。つまり、重み関数が、以下に示すように、w(x)dx がランダム変数 x が x − x + dx に含ま
れる確率になる。そこで、ある特定の分布をもつ乱数を如何にして生成するか、ということが
問題になる。
変換法（Transformation method）
これまでの例題と同様に、ある特定の分布をもつ乱数が解析的に求めることができる場合を
変換法（Transformation method）という。
例として、被積分関数に減衰する項が含まれている場合を考えよう。このとき、減衰項を重
み関数 w(x) として選ぶ。
∞
I =
0
∞
dxf (x) =
1
=
0
0
dxe−x g(x)
dyg(y)
ここで、
dy
= c1 e−x
dx
定数 c1 は、重み関数に対する規格化条件
y(x = 0) = 0,
y(x = ∞) = 1
により決めることができる。従って、c1 = 1 となり、
y(x) = 1 − e−x
つまり、確率密度が P (x) = e−x の分布をもつ乱数 x は、
x(y) = − log(1 − y)
1. 数値解析の基礎
38
によって生成される。ここで、y は [0, 1] の乱数である。
別の例として、極座標系における多重積分を考えよう。
I=
d3 rf (x, y, z) =
∞
0
drr2
π
0
2π
dθ sin θ
0
dφf (r, θ, φ)
この場合、角度に関する積分は、つぎのような重み関数を用いた乱数を生成して評価する。
dyθ
= c1 sin θ
dθ
dyφ
= c2
dφ
勿論、定数 c1 , c2 は、yθ , yφ に対する規格化条件から決まる。これらを解いて、
yθ =
1
(1 − cos θ)
2
⇒
θ = cos−1 (1 − 2yθ )
φ
⇒ φ = 2πyφ
2π
となる。ここで、yθ , yφ は、[0, 1] の一様乱数である。
yφ =
棄却法（Rejection method）
解析的にある確率密に従った乱数が求まらないことも多い。例えば、演習問題 1-6-2 におけ
る重み関数として、
6
x2
w(x) =
1−
5
2
とすると、
6
x2
y(x) = x x −
5
6
となり、x についての３次関数になる。この場合、x = x(y) というかたちに変形することがで
きない。このようなときは、以下のようなノイマンの棄却法（Rejection method) を用いる。
重み関数（確率密度）w(x) が与えられたとする。この関数の x ∈ [0, 1] での最大値を wmax
とする。
1. [0, 1] の２つの乱数 xi , ηi を発生させる。
2. w(xi )/wmax と ηi を比較する。
3. w(xi )/wmax > ηi のときは、この xi をランダム変数として採用する。
4. w(xi )/wmax < ηi のときは、この xi をランダム変数として棄却する。
このようにして、重み関数（確率密度）w(x) に対応したランダム変数を生成することができる。
正規分布をした乱数の生成
平均がゼロで、分散が σ = 1 であるような正規（Gauss）分布は
x2
1
w(x) = √ e− 2
2π
1. 数値解析の基礎
39
で与えられる。このような確率密度をもつランダム変数 x を生成することは、理工学分野での
数値解析において非常に多い。しかしながら、この確率密度について変換法を用いて新たなラ
ンダム変数 y の関数として、x = x(y) を解析的に求めることはできない（誤差関数の逆関数が
必要になる）。勿論、うえの棄却法を用いれば正規（Gauss）分布に対応する乱数を生成するこ
とは可能である。しかしながら、幸いにも、正規（Gauss）分布になるランダム変数の生成に
は、いくつか簡単なアルゴリズムが考え出されている。
中心極限定理を用いた以下の方法は、特に簡単で有効である。si ∈ [0, 1] をランダム変数と
して、
ξ = s1 + s2 + · · · + sN
とする。まえと同様の議論により、新たなランダム変数 ξ の平均値 ξ
と分散 σξ2 = ∆ξ 2 は、
それぞれつぎのように与えられる。
ξ
= N s
σξ2 = N ∆s2
中心極限定理によれば、N が大きい場合、ランダム変数 ξ に対する確率密度 Pξ (ξ) は、つぎの
正規分布に近づく。
1
Pξ (ξ) = e
2πσξ2
−
(ξ−ξ)2
2σ 2
ξ
(ξ−N s)
−
1
2
2N σs
=
e
2πN σs2
2
(30)
ランダム変数 s は、[0, 1] の一様乱数であるから、s に対する確率密度 Ps (s) は
Ps (s) = 1
で与えられる。従って、
s
=
σs2
= ∆s
2
2
= s
1
0
dssPs (s) =
2
− s
=
1
0
1
2
dss2 Ps (s) − s
2 =
1
12
そこで、N = 12 とすれば、式（30）から、ランダム変数 ξ に対する確率密度 Pξ (ξ) はつぎのよ
うな正規分布になる。
(ξ−6)2
1
Pξ (ξ) = √ e− 2
2π
つまり、１２個の乱数を生成し、その和 ξ から 6 を引いたものが、平均がゼロで分散が 1 の正
規分布に対応する乱数になる。
具体的なアルゴリズム（プログラム）は以下のようになる。
xi = 0
for i = 1, 12
xi = xi + rand(i)
next i
xi = xi - 6
1. 数値解析の基礎
40
平均が ξ
で分散が σ の正規分布
Pξ (ξ) = √
1
2πσ 2
e−
(ξ−ξ)2
2σ 2
に対応する乱数にするには、うえのアルゴリズムで求めた ξ に、σ をかけたうえで、ξ
を加え
ればよい。
他にもいくつかアルゴリズムが存在するが、それらについてはモンテカルロ法の文献を参照
してもらいたい。
演習問題 1-6-3
うえで説明した正規分布をした乱数を生成するアルゴリズムを用いて、ランダム変数 x を 100000
個生成し、横軸を x、縦軸をそれぞれの x に対する頻度とするヒストグラムを作成せよ。
求めたヒストグラムのうえに正規分布をプロットすることにより、正しくランダム変数 x が
生成されていることを検証せよ。
2. デバイスモデリング入門
2
41
デバイスモデリング入門
本講義の後半では、これまで学習した数値計算の手法をいよいよ具体的な半導体デバイス特
性シミュレーションに応用するため、
「デバイスモデリング」の基礎（の基礎！）について学習
する。
「デバイスモデリング」（デバイス・シミュレーション）は、実際の半導体デバイスの設計
においてデバイス特性の予測を行い、設計指針を与えることを大きな目的としている。デバイ
ス・シミュレーションには、現在のところ、大きく分けて
(i) 連続的描像のもとでの流体デバイスシミュレーション
(ii) 離散的描像のもとでの粒子デバイスシミュレーション
前者の代表的なものがドリフト拡散法であり、また後者の代表的なものがモンテカルロ法と呼
ばれるものである。これまでのデバイス・シミュレーションは、ドリフト拡散法が中心である
（現在形であることに注意！）。
しかしながら、近年の半導体デバイスの微細化（特に、Si を基盤にした既存の Si-MOSFET
におきては、研究所レベルでゲート長が既に十数 nm にまで達している）に伴って、デバイス
の電気的特性を決める電子（或いは、正孔）の輸送現象はドリフト拡散法では解析することが
困難なほどに物理的複雑さを深めている。例えば、空間的な不均一性に伴った非局所な効果や
電子の波動性に伴った量子効果、高濃度にドープされることにより生じる多体効果などである。
その結果、デバイス特性を精度良く解析するためには、上述の複雑な物理機構をふまえたモン
テカルロ法や（新たな？）シミュレーションが必要不可欠になってきている。従って、現在の
デバイス・シミュレーション研究では、デバイス内での電子輸送の物理現象に基づいた物理モ
デルの構築が大きな比重を占めている。
勿論、これらのデバイスモデリングの最先端まで本稿で議論するのは（時間的にも、必要と
される物理的背景の多さからも）不可能である。そこで、本章では、デバイス・シミュレーショ
ンとして最も汎用性が高く、またデバイス設計の現場で実際に日常的に使用されているドリフ
ト拡散法について学習する。ドリフト拡散法は、デバイス基本方程式と称されるポアソン方程
式と電流連続式を離散化して自己無撞着に解く手法であり、その手法自体にさまざまな数値計
算上のトリックがちりばめられている。そこで、単純ではあるが物理的内容が豊富で半導体デ
バイスの基礎である pn ダイオードを題材として、デバイス特性シミュレーションを実際に行
いながら、ドリフト拡散法の手法を学ぶ。また同時に、デバイスシミュレーションを行うこと
によって初めて得られる物理的内容の理解の深みというものを是非実感してもらいたい。
デバイスシミュレーションについての解説（種類とその基礎付け）と最近の進展については、
拙著の以下の記事を参考にされたい。
“デバイスシミュレーションとその物理”
会、2002)
応用物理学会誌 7 月号, pp906-910 (応用物理学
2. デバイスモデリング入門
2.1
2.1.1
42
デバイス特性シミュレーション：ポアソン方程式
MOS キャパシター: 陽解法
平衡状態でのデバイス特性解析では、ポアソン方程式が重要な役割を果たす。このような非
常に単純な系においても、デバイス数値解析における困難さが出現する。これは、半導体デバ
イス内でのキャリア濃度等がポテンシャルに対して指数関数的に数桁にわたって変化すること
による。そこで、MOS キャパシターにおける蓄積層、空乏層、反転層の形成に伴った空間電荷
や静電ポテンシャルを、ポアソン方程式を離散化することにより求めてみる。
Si の基板は一様にドープされた p 型であるとして、Si 基板の深いところ（つまり酸化膜界面
から離れたところ）での静電ポテンシャルをゼロにとる。また原点を Si 基板と酸化膜との界面
にとれば、ポアソン方程式は、
q d2 ψ(x)
+
−
=
−
p(x)
−
n(x)
+
N
(x)
−
N
(x)s
a
d
dx2
ε0
(31)
で与えられる。ここで、q は素電荷で、ε0 は Si の誘電率である。また p(x), n(x), Nd+ (x), Na− (x)
は、それぞれ界面から距離 x での正孔濃度、電子濃度、イオン化したドナー不純物濃度、イオン
化したアクセプター不純物濃度である。上述のように、静電ポテンシャル ψ(x) は、ψi (x)−ψi (∞)
で定義され、ψi (x) は真性ポテンシャルで、真性エネルギーレベル Ei を用いて、
ψi (x) = −
Ei
q
のように定義される。静電ポテンシャル ψ(x) を用いて、正孔および電子の濃度は
qψ(x)
q
(ψB − ψ(x)) = Na e− kT
p(x) = ni exp
kT
(32)
n2 qψ(x)
q
(ψ(x) − ψB ) = i e kT
n(x) = ni exp
kT
Na
(33)
となる。ただし、縮退の効果は無視した。ここで、ni は真性キャリア濃度、ψB はフェルミポ
テンシャルと真性ポテンシャルとの差であり、
kT
Na
ln
ψB =
q
ni
である。基板の深いところでの電荷中性条件から、
Nd+ (∞) − Na− (∞) ≈ −Na +
n2i
Na
であるから、ポアソン方程式は以下のように書ける。
qψ(x)
qψ(x)
q
n2i
d2 ψ(x)
− kT
kT
=
−
N
e
−
1
−
e
−1
a
dx2
ε0
Na
左辺を３点公式を用いて離散化すれば、
ψn+1 = 2ψn − ψn−1 − h2
ni qψn
qni Na − qψn
e kT − 1 −
e kT − 1
ε0 ni
Na
(34)
2. デバイスモデリング入門
43
0.5
Electrostatic Potential (V)
0.4
0.3
0.2
0.1
Qs ( x 10-7 C/cm2)
0
2
3
4
5
implicit
-0.1
-0.2
-0.3
0
10
20
30
40
Depth from interface x (nm)
図 5: ポアソン方程式を離散化したものを数値積分して得られる Si-MOS キャパシター基板
内での静電ポテンシャル。用いた境界条件は ψs = 0.4 V で、パラメータとして空間電荷密度
Qs = 5, 4, 3, 2 × 10−7 C/cm2 とした場合を示す。
となり、これまで学習した数値解析の手法を用いて解くことができそうである。実際、本課題
での境界条件は、
ψ(x = 0) = ψs , ψ(x = ∞) = 0
となる。ここで、ψs は酸化膜界面での表面ポテンシャルである。また、２番目の境界条件は電
荷中性条件により界面から離れたところでポテンシャルが一定であることを表す。しかしなが
ら、実際に静電ポテンシャルを求めたうえでなければ ψ(x = ∞) = 0 を用いることができない
（射撃法を用いて原理的には解くことが可能ではある）。そこで、界面でのポテンシャルの１階
微分
Qs
dψ = −Es = −
dx x=0
ε0
をパラメータとして、
ψ0 = ψs ,
ψ1 ≈ ψs − h ·
Qs
ε0
を境界値として積分し、十分に x が大きいところで静電ポテンシャル ψ(x) がゼロになる解を
見つければ良い。
実際には、このような指針で静電ポテンシャルの数値解を求めることは非常に困難であるこ
とがわかる。図 5 にいくつかの初期（条件）のもとで求めた静電ポテンシャル ψ を示す。ここ
では、表面ポテンシャルを ψs = 0.4 V と固定したうえで、空間電荷密度 Qs をパラメータとし
て変化させて無限遠でポテンシャルがゼロになる解を求めようとしたものである。また、基板
2. デバイスモデリング入門
44
のアクセプター不純物濃度は Na = 1018 cm−3 とした。
演習問題 2-1-1
うえで説明した方法を用いて、Si-MOS キャパシターでの基板内での静電ポテンシャル ψ を求
めるプログラムを作成し、図 5 の結果を自ら確めよ。p-type の Si-MOS で基板のアクセプター
不純物濃度は Na = 1018 cm−3 とせよ。
2.1.2
陰解法
前述したように、デバイス内での物理量が指数関数的に依存することから、離散化したポア
ソン方程式を単純に積分することは困難である。一般にデバイス・シミュレーションにおいて
は、このような陽解法に代わって陰解法を用いて解く場合が多い。
離散化したポアソン方程式をもう一度書き下すと、
qni Na − qψn
ni qψn
ψn+1 − 2ψn + ψn−1
=−
e kT − 1 −
e kT − 1
h2
ε0 ni
Na
(35)
境界条件は、
ψ(x = 0) = ψ1 = ψs ,
ψ(x = L) = ψN +1 = 0
ここで、L は基板の深さ方向の大きさで空乏層幅よりも十分に大きくとる。また、N = L/h と
した。うえのポアソン方程式は静電ポテンシャル ψ について非線形（指数関数の部分）であり、
これが前述の陽解法での解析の困難さに繋がっている。そこで、静電ポテンシャル ψ について
以下のように線形化する。
ψn = ψn0 + δψn
として指数関数の部分を展開すると、線形化したポアソン方程式は
δψn−1 + σn δψn + δψn+1 = rn
となる。ここで、
q 2 ni
σn = −2 − h
ε0 kT
qni
0
0
rn = −ψn−1
+ 2ψn0 − ψn+1
− h2
ε0
2
(36)
ni qψn0
Na − qψn0
e kT +
e kT
ni
Na
0
qψn
ni
Na − qψn0
kT
kT
e
−1 −
e
−1
ni
Na
である。
演習問題 2-1-2
うえの線形化したポアソン方程式 (36) を導け。
うえのポアソン方程式 (36) は静電ポテンシャル ψ に対して線形であるから、以下のような
行列のかたちにまとめることができる。
Sδψ = R
(37)
2. デバイスモデリング入門
45
1 0-5
Charge Density (c/cm2)
MOS-C
Na = 1018 1/cm3
1 0-6
1 0-7
-0.5
0
0.5
1
1.5
Surface Potential (V)
図 6: 陰解法によって求めた電荷密度 Qs の絶対値の表面ポテンシャル ψs 依存性。
ここで、




S=


σ2 1
1 σ3
.
..
0
...
0 0
1 0
..
.

0

0 

,






δψ = 


0 1 σN
δψ2
δψ3
..
.
δψN




,






R=


r2
r3
..
.







rN
境界条件から、
δψ1 = δψN +1 = 0
であることに注意せよ。従って、S の逆行列を求めて R にかければ、各場所における静電ポテ
ンシャル ψ を求めることができる。しかしながら、２次元或いは３次元シミュレーションにな
ると行列 S はかなり大きくなり、直接逆行列を求めるのは現実的ではない。うえのような行列
S を３重対角行列と呼び、後述するような効率的に解くアルゴリズムが存在している。
実際にうえのポアソン方程式を解くには、最初にデバイス内での静電ポテンシャル ψn0 の形
状を適当に決めておき、行列 S および R を計算したうえで行列方程式を解く。得られた静電ポ
テンシャルの差分 δψn を最初の計算に用いた静電ポテンシャル ψn0 に加えて、改めて行列方程
式を解く。静電ポテンシャルの差分 δψn が静電ポテンシャル ψn0 に比べて十分に小さくなるま
でこの操作を（iterative に）くり返す。十分に収束したところで、正しい静電ポテンシャル ψn
を求めることができる。図 6 は、このような陰解法によって求めた電荷密度 Qs の絶対値の表
面ポテンシャル ψs 依存性を示す。表面ポテンシャル ψs の増大に応じて、蓄積層、空乏層およ
び反転層の形成に伴った電荷密度 Qs の変化が読み取れる。解析的に導出が可能な教科書で見
られる通常の空乏近似では、蓄積層あるいは反転層におけるキャリア濃度の増大による電荷密
度 Qs の変化を導くことができないことに注意せよ。
2. デバイスモデリング入門
2.1.3
46
３重対角行列の再帰的な解法
３重対角行列の逆行列化は、以下のような再帰的な方法（recursive method）によって解か
れる。
３重対角行列を含む式 (37) のような行列方程式は、一般的に、以下のような行列で書き表わ
される。












B2 C2 0
0
0
A3 B3 C3
0
0
0 A4 B4
C4
0
..
..
..
..
.
.
.
.
0
...
AN −1 BN −1 CN −1
0
...
0
AN
BN
 
δy2
δy3
δy4
..
.









 δy
N −1














F2
F3
F4
..
.
 
 
 
 
 
=
 
 
 
  F
  N −1
δyN












FN
δyn を中心とした方程式は、
An δyn−1 + Bn δyn + Cn δyn+1 = Fn ,
2≤n≤N
(38)
と書くことができる。x = 0 および x = L での境界条件から、
δy1 = δyN +1 = 0
である。そこで、いま式（38）が以下のような変数２個（δyn と δyn+1 ）の方程式に変換でき
たと仮定する（これは境界条件から可能であると予想される）。
Bn δyn + Cn δyn+1 = Fn
これから、
δyn =
(39)
Fn
Cn
−
δyn+1
Bn
Bn
(40)
となる。これを式（38）に代入すると、
C
Bn − An n−1
Bn−1
δyn + Cn δyn+1 = Fn − An
Fn−1
Bn−1
これが式（39）と同等であるためには、以下のような漸化式が成立しなくてはならない。
Bn = Bn − An
Cn−1
,
Bn−1
Cn = Cn ,
Fn = Fn − An
Fn−1
Bn−1
(41)
つまり、n に対する Bn , Cn , Fn がわかれば、次の次数の Bn+1
, Cn+1
, Fn+1
を計算することがで
きる。式（38）において n = 2 とすれば
B2 δy2 + C2 δy3 = F2
となるから、式（39）と比較して、
B2 = B2 ,
C2 = C2 ,
F2 = F2
2. デバイスモデリング入門
47
となる。従って、漸化式（41）を用いて 3 ≤ n ≤ N に対する Bn , Cn , Fn を計算することがで
きる。
すべての Bn , Cn , Fn がわかれば、式（40）と境界条件 δyN +1 = 0 を用いて δyn を 2 ≤ n ≤ N
について求めることができる。
これが、再帰的な方法（recursive method）と呼ばれるものである。まともに行列の逆行列
を求めて解く場合に比べて、必要とされる配列の大きさが大幅に小さくなっていることに注意
して欲しい。このような方法を用いて、現在のデバイス・シミュレーションでは百万 × 百万程
度の行列計算を行うことも珍しくない。
演習問題 2-1-3
本文で説明した iterative な方法により、静電ポテンシャル δψ についての行列式を解くプログ
ラムを作成せよ。再帰的な行列計算プログラムは、以下のホームページ
http://hermes.esys.tsukuba.ac.jp/ sano/cprog
にポストされているもの（tridag.c）を参考にせよ。デバイスの構造は演習問題 2-1-1 と同じに
する。また、デバイスの深さ方向の大きさは L = 100 nm とせよ。
作成したプログラムを用いて、表面ポテンシャルが ψs = 0.4 V のときの静電ポテンシャル ψ
のデバイスの深さ方向の形状を求めよ（図 5 参照）。さらに、さまざまな表面ポテンシャルに対
しての電荷密度 Qs を求め、図 6 のようになることを確めよ。
次に、半導体デバイスの教科書に出ている空乏近似のもとで電荷密度 Qs についての解析式
を導け。そして、この解析式を用いて電荷密度 Qs を求めて数値解析の結果と比べよ。
2. デバイスモデリング入門
2.2
2.2.1
48
デバイス特性シミュレーション：ドリフト拡散法
半導体デバイスの基本方程式
前章では、平衡状態におけるデバイス特性の数値解析について議論したが、デバイスが実際
にオンオフ等の動作する状態は、非線形性の非常に強い非平衡状態である。従って、ボルツマ
ン輸送方程式を解析的に検討するのは一般には困難である。そこで、デバイス特性解析には数
値シミュレーションが大きな役割を果たす。本章では、デバイス・シミュレーション（ドリフ
ト拡散法）で基本となるいわゆるデバイス方程式について議論する。
半古典的な範疇において（現在の極微細デバイスではこの仮定が大きく揺らいでいるが）、非
平衡な電子の運動（輸送）はボルツマン輸送方程式と呼ばれる方程式に基づく。ボルツマン輸
送方程式についてここでは深く言及しないが、これは物理的運動学（kinetics）の基本方程式で
あり、古典的リュービル方程式からある近似のもとで導出される。しかしながら、ボルツマン
輸送方程式自体は複雑な偏微分積分方程式であり、ごく単純なモデルを除いては解析的に解く
ことはできない。ボルツマン輸送方程式は、モンテカルロ法と呼ばれる方法によりデバイス・
シミュレーションに組み込まれる。一方、ボルツマン輸送方程式に対するモーメントの輸送方
程式を基本方程式とする方法が、本セクションで述べるドリフト拡散法である。
ドリフト拡散法は、キャリアの輸送を規定する輸送方程式とデバイス内でのポテンシャル形
状を規定するポアソン方程式から成る。ポアソン方程式は、前セクションでの静的なデバイス
特性解析に用いたものと同一である。１次元近似のもとで、
d2 ψ(x, t)
q +
−
=
−
p(x,
t)
−
n(x,
t)
+
N
(x)
−
N
(x)
a
d
dx2
ε0
(42)
キャリア（電子および正孔）の輸送方程式は、ボルツマン輸送方程式の第０次モーメントであ
る電流連続式を意味する。１次元近似での正孔および電子の電流連続式は
1 ∂Jp (x, t)
∂p(x, t)
=−
+ Gp (x, t) − Rp (x, t)
∂t
q
∂x
(43)
1 ∂Jn (x, t)
∂n(x, t)
=
+ Gn (x, t) − Rn (x, t)
∂t
q
∂x
(44)
となる。ここで、p(x, t), n(x, t) および Jp (x, t), Jn (x, t) は正孔および電子の数密度と電流密度
である。G(x, t), R(x, t) は単位体積、単位時間当たりのキャリアの生成および消滅率である。
正孔および電子の電流密度はドリフト電流と拡散電流から成り、以下のように書ける。
Jp (x, t) = −qp(x, t)µp
∂ψ(x, t)
∂p(x, t)
− qDp
∂x
∂x
∂ψ(x, t)
∂n(x, t)
+ qDn
∂x
∂x
ここで、µp , µn は正孔および電子の移動度、Dp , Dn は拡散係数である。これらは、アインシュ
タインの関係式を用いて、
q
q
µp =
D p , µn =
Dn
kT
kT
と表される（これを揺動散逸定理とも言う）。
Jn (x, t) = −qn(x, t)µn
2. デバイスモデリング入門
49
移動度 µ、キャリアの生成率 G(x, t) および消滅率 R(x, t) は対象とする半導体デバイスの物
質パラメータであり、後述するように、さまざまな物理モデルが存在している。拡散係数 D は
アインシュタインの関係式を用いて移動度から求められる。従って、ドリフト拡散法における
基本方程式 (42)-(44) は、
p(x, t), n(x, t), ψ(x, t)
つまり、正孔および電子の数密度と静電ポテンシャルを未知変数とする偏微分方程式になってい
る。ドリフト拡散法ではこれらの方程式を離散化して数値シミュレーションを行うことになる。
一方、正孔および電子に対する quasi-Fermi ポテンシャルを以下のように
φp (x, t) = ψ(x, t) +
p(x, t)
kT
ln
q
ni
φn (x, t) = ψ(x, t) −
n(x, t)
kT
ln
q
ni
とそれぞれ定義すると、正孔および電子に対する電流密度の式は、
Jp (x, t) = −qp(x, t)µp
∂φp (x, t)
∂x
∂φn (x, t)
∂x
のように、ドリフト電流と拡散電流をまとめて記述することができる。そこで、正孔および電
子に対する quasi-Fermi ポテンシャルと静電ポテンシャル
Jn (x, t) = −qn(x, t)µn
φp (x, t),
φn (x, t),
ψ(x, t)
を未知変数としてドリフト拡散方程式を解く場合もある。
quasi-Fermi ポテンシャルの考え方は、非平衡状態でのキャリア輸送を考えるうえで非常に
重要なものである。例えば、平衡状態で電流が流れていない平衡状態 Jp = Jn = 0 では、正孔
および電子に対する quasi-Fermi ポテンシャルは
φp = φn = const.
となり、通常の Fermi ポテンシャル（或いは Fermi エネルギー）に一致する。
φp = φn = ψf = −
Ef
q
一方、非平衡状態においては、準平衡状態にあるコンタクト領域を除いて、正孔および電子の
quasi-Fermi ポテンシャルは一致しない。このとき、キャリアの局所的な数密度は
n(x, t)p(x, t) =
となる。
n2i exp
q
[φp (x, t) − φn (x, t)]
kT
2. デバイスモデリング入門
2.2.2
50
pn ダイオードへの応用とドリフト拡散法のアルゴリズム
quasi-Fermi ポテンシャル（擬フェルミ準位）
pn 接合のようなダイオード構造に対して、うえの quasi-Fermi ポテンシャルの考え方を応用
するのは興味深い。特にダイオードの内部で、電子および正孔に対する quasi-Fermi ポテンシャ
ルがどのように変化するかを理解しておくことが望ましい。
pn ダイオードのコンタクトの両端での quasi-Fermi ポテンシャルの差
φp (0) − φn (L) = Vapp
が、ダイオードにかかる印加電圧 Vapp に対応する。ここで、L はダイオードの長さである。
演習問題 2-2-1
pn ダイオードの内部における電子および正孔に対する quasi-Fermi ポテンシャルがどのよう
に場所の関数として変化するか、考察せよ。順バイアスと逆バイアスの場合で区別して述べる
こと。
境界条件
ドリフト拡散シミュレーションにおける境界条件を明確にするために、具体的なデバイス構造と
して１次元 pn 接合を考えよう。p− および n− 領域におけるコンタクトをそれぞれ x = 0, x = L
とする。これらのコンタクトでは準熱平衡状態にあると仮定する。これらの領域での電荷中性
条件から、
Nd+ (0) − Na− (0) + p(0) − n(0) = 0,
p(0)n(0) = n2i ,
Nd+ (L) − Na− (L) + p(L) − n(L) = 0
(45)
p(L)n(L) = n2i
(46)
これらから、コンタクトでの正孔と電子の境界条件は、


$
2 

Γ(0)
2ni
p(0) = −
1+ 1+
,
2 
Γ(0) 

Γ(L) 
1+
n(L) =
2 
$
1+
2ni
Γ(L)
n(0) =
n2i
p(0)
(47)
p(L) =
n2i
n(L)
(48)

2 

,
ここで、Γ(x) = Nd+ (x) − Na− (x) である。
また、p− 領域でのコンタクトの電位を Vapp 、n− 領域でのコンタクトの電位を 0 とする。
quasi-Fermi ポテンシャルの定義から、静電ポテンシャル ψ に対する境界条件は以下のように
なる。
p(0)
n(L)
kT
kT
ln
ln
ψ(0) = φp (0) −
, ψ(L) = φn (L) +
q
ni
q
ni
ここで、
Vapp = φp (0) − φn (L),
φp (0) = φn (0),
φn (L) = φp (L)
2. デバイスモデリング入門
51
である。一般には、静電ポテンシャル ψ の原点を φn (L) に取ることにより以下のような境界条
件を用いることが多い。
ψ(0) = {φp (0) − φn (L)} −
ψ(L) =
p(0)
p(0)
kT
kT
ln
ln
= Vapp −
q
ni
q
ni
n(L)
kT
ln
q
ni
ドリフト拡散法のアルゴリズム
ドリフト拡散法の定常状態のもとでの簡単化した流れ（アルゴリズム）は、以下のようにま
とめることができる。
(1) デバイス構造とモデルパラメータの確定
(2) デバイス内部の変数（ψ(x, t), p(x, t), n(x, t)）の初期設定
(3) 平衡状態（Vapp = 0）における計算の実行
(4) 以下の計算を適当な精度のもとで、最終の印加電圧（Vapp = Vf inal ）まで ∆V ご
とに繰り返す
(i) ポアソン方程式を解く
(ii) 電子に対する電流連続方程式を解く
(iii) 正孔に対する電流連続方程式を解く
(iv) 電圧を上げる: Vapp + ∆V → Vapp
(5) 計算結果の出力
ここで、(4) における３つのドリフト拡散方程式は、それぞれ強い非線形性をもつことから、
方程式を線形化したうえで、前セクションでポアソン方程式を解いたような陰解法を用いて解
く。方程式の線形化や差分化の詳細については後述する。
2.2.3
ドリフト拡散法における物理モデル
移動度モデル
半導体でのキャリア輸送現象を記述するのはドリフト拡散においては移動度である。線形応
答の範囲では、正孔および電電子の移動度はドリフト速度を用いて
vp = µp E,
vn = µn E
のように与えられ、定数である。しかし、デバイス動作のような非線形な領域にまでこのよう
に定義された移動度を拡張すると、移動度は非常に複雑な電界等の依存性をもつ。キャリアの
移動度は半導体中でのフォノン散乱、不純物散乱、表面ラフネス散乱等によって規定されるから
2. デバイスモデリング入門
52
µp
6.3 × 1016
0.76
495
47.7
1.95 × 104
1
N0 (cm3 )
α
µmax (cm2 /V sec)
µmin (cm2 /V sec)
E0 (V/cm)
β
µn
8.5 × 1016
0.72
1330
65
8.0 × 103
2
表 1: Si における移動度モデルパラメーター
であり、このような散乱仮定は電界に対して強い非線形な依存性をもつ。バルク形状で室温程
度での最も支配的な散乱機構はフォノン散乱であり、あまり強くない電界領域では、Si での電
子移動度は 1400 cm2 /V sec に達する。しかしながら、不純物を多くドープした系や MOSFET
構造のような場合は、不純物散乱や表面ラフネス散乱の寄与が大きくなり、Si での電子移動度
は 400 cm2 /V sec 程度にまで抑制される。このような不純物濃度や表面ラフネス等の影響を考
慮に入れながら正しい移動度を再現することはかなり困難をともなう。そこで、正孔や電子に
対する移動度モデルとして、これまでに非常に多くの定式化がなされている。なかでも、山口
氏の提案した山口モデルは世界的にも最も広く用いられている移動度モデルの一つである。
ここでは、Si に対してよく用いられる移動度モデルを紹介する。不純物濃度を N 、電界を E
として、移動度 µ は、

µ=
µmax − µmin
1+
N
N0
α

1
+ µmin  1+
E
E0
β 1/β
(49)
式に含まれる Si に対するパラメータを表 1 に与える。
生成再結合モデル：Shockley-Read-Hall model
キャリアの生成再結合過程として最も基本的なものが Shockley-Read-Hall model と呼ばれ
るものである。これは、半導体のバンドギャップ中に存在するエネルギー準位にキャリアが捕
獲、あるいはエネルギー準位から励起してキャリア数を変動させる過程である。バンドギャッ
プ中に存在するエネルギー準位は、欠陥や不純物等により生じる。詳細は半導体デバイスの教
科書（例えば Sze の教科書「Semiconductor Devices:Physics and Technology」p. 49）に譲る
として、非平衡状態かつ定常状態における単位時間、単位体積当たりの生成再結合率 U (x) は、
U (x) = Up (x) = Un (x) = −
p(x)n(x) − n2i
τp [n(x) + nt (x)] + τn [p(x) + pt (x)]
(50)
で与えられる。ここで、pt (x), nt (x) はトラップ準位にフェルミ準位があったときの正孔および
電子の数密度、τp , τn は正孔および電子の捕獲されて消滅するまでの平均寿命を表す。キャリア
の数密度 p(x)n(x) が真性キャリア濃度の２乗 n2i よりも大きいときは U (x) は負となり、キャリ
アの捕獲再結合を表す。また逆のときは、キャリアの生成を表す。正孔および電子の捕獲され
2. デバイスモデリング入門
53
て消滅するまでの平均寿命が等しい場合は、うえの生成再結合率 U (x) はトラップ準位が真性
エネルギー準位に等しい（Et = Ei ）ときに最大になる。つまり、Shockley-Read-Hall model
では、トラップ準位がミッドギャップにある場合は正孔および電子の生成再結合中心として働
き、トラップ準位が浅い場合はただの捕獲中心として働くことを意味する。ドリフト拡散方程
式では、式 (43),(44) における生成消滅項 G(x) − R(x) をうえの U (x) で置き代える。
オージェ再結合モデル
キャリアの再結合過程としてもうひとつ基本的なものがオージェ再結合過程である。これは、
トラップ準位を経ずに電子正孔の再結合が起こる過程である。電子正孔の再結合の結果として、
バンドギャップに相当するエネルギーがキャリアに与えられ、高エネルギー・キャリアが生成
される。単位時間単位体積におけるオージェ再結合率は、以下のようにモデル化されることが
多い。
Uaug (x) = Up (x) = Un (x) = −r{n(x)2 p(x) + p(x)2 n(x)}
(51)
ここで、Si においては r = 2 × 10−31 cm6 /s である。上式から明らかなように、オージェ再結
合過程ではキャリア密度が大きいときに最も顕著であり、キャリア数が多いパワーデバイスに
おいて特に重要になる。
衝突イオン化モデル
電界が非常に強く衝突イオン化過程が無視できない状況に於いては、Shockley-Read-Hall
model に加えて衝突イオン化項を生成項として考慮する必要がある。衝突イオン化過程は、高
エネルギー化した電子（正孔）が束縛状態にある充満帯の電子（正孔）を価電子帯に励起し、正
孔と電子（電子と正孔）を生成する過程である。電場が十分に大きい場合は、衝突イオン化に
よって新たにキャリアが生成され、その生成されたキャリアがさらに衝突イオン化によって新た
なキャリアを生成する、といった過程が連続的に起こる（avalanche 過程）。その結果、キャリ
ア数が急激に増大し、半導体デバイスの破壊に繋がる。ドリフト拡散法のなかでの衝突イオン
化のモデル化はいまでも大きな研究課題の一つであり、世界各国で活発に研究が行われている。
最も単純で原始的なモデルが、局所電界モデルと呼ばれるものである。正孔と電子に対する
イオン化係数（単位長さ当たりの衝突イオン化を起こす数）を αp , αn とすれば、単位時間と単
位体積当たりのキャリア生成率 G(x, t) は、
G(x, t) = Gp = Gn =
1
[αn |Jn | + αp |Jp |]
q
と表される。局所電界モデルでのイオン化係数 αp , αn は、
αp or αn = A exp −
b
|E(x, t)|m
(52)
(53)
となる。ここで、E(x, t) は局所的電界強度である。含まれるパラメータ A, b は、Si において表
2 のように与えられる。このイオン化モデルは一般に広く用いられているが、局所的な平衡状
態を仮定しており電界の変化が急峻で非局所な効果が重要な不均一な微細構造の半導体ではイ
オン化大きめに見積もられてしまうので注意が必要である。
2. デバイスモデリング入門
54
電 子
正 孔
A (1/cm)
3.80 × 106
2.25 × 107
b (V/cm)
1.75 × 106
3.26 × 106
m
1
1
表 2: Si における衝突イオン化パラメーター
2.2.4
ドリフト拡散方程式の離散化
上述のデバイス方程式を離散化して解くドリフト拡散シミュレーションの具体例として、pn
接合における定常状態での DC 特性解析を試みよう。
解くべきデバイス基本方程式は、（定常状態であるから）電流連続式における時間依存性を
ゼロにおいた方程式となる。
Jp (x) = −qp(x)µp
∂ψ(x)
∂p(x)
− qDp
∂x
∂x
(54)
Jn (x) = −qn(x)µn
∂ψ(x)
∂n(x)
+ qDn
∂x
∂x
(55)
1 ∂Jp (x)
− G(x) + U (x) = 0
q ∂x
(56)
1 ∂Jn (x)
+ G(x) − U (x) = 0
q ∂x
(57)
q d2 ψ(x)
+
−
=
−
p(x)
−
n(x)
+
N
(x)
−
N
(x)
a
d
dx2
ε0
(58)
これらの方程式を離散化するために、想定している pn 接合の長さを L、領域を N 個の等間
隔メッシュに区切ったとする。ただし、実際のデバイス・シミュレーションにおいては、キャリ
ア濃度や静電ポテンシャル等の変化はデバイス内の構造に大きく依存しており、これらの量が
急峻に変化する領域（例えば接合領域）では十分に小さいメッシュを用いる必要がある。従っ
て、不均一なメッシュを用いるほうがデバイス・シミュレーションではより一般的である。い
まの簡単化した系において、p− および n− 領域でのコンタクトは、
x1 = 0,
xN +1 = L
となる。さらに、補助的（auxiliary）なメッシュ（添字を m で表す）がそれぞれの上述のメッ
シュ点の中間にあると想定する（図 7 参照）。未知変数 p(x), n(x), ψ(x) は上述のメッシュ上で
計算され、これらの微分（例えば、電流密度）は補助的（auxiliary）なメッシュ上で計算する
ことにする。
単純に電流密度の式を離散化して電流連続式に代入すると、前セクションの MOS-C で問題
となったような数値計算上の不安定性が生じる。これは、各メッシュごとの変化量が静電ポテ
ンシャルとキャリア濃度では大幅にスケールが異なることによる。つまり、キャリア濃度は指
数関数的に静電ポテンシャルに依存する。そこで、各メッシュ内においては、変化がキャリア
濃度のそれに比べて小さいと考えられる静電ポテンシャル、移動度（拡散係数）、電流密度は
2. デバイスモデリング入門
55
p-contact
n-contact
auxiliry mesh
x=0
1
1
2
x=L
M-2
2
N-2
M-1
N-1
M
N
N+1
main mesh
図 7: main メッシュと auxiliry メッシュ。
メッシュ内においては一定であると仮定する。この方法を Scharfetter-Gummel の方法と呼び、
ドリフト拡散シミュレーションのその後の発展の契機となったものである。メッシュ内で積分
すれば、正孔の電流密度に対する式は以下のようになる。
qEx
p(x) = p(0)e kT +
これを書き換えれば、
Jp = −
qEx
Jp 1 − e kT ,
qµp E
qµp E
1−e
qEh
kT
p(0)e
qEh
kT
with E = −
dψ
dx
− p(h)
となる。ここでメッシュ間隔を h とした。これは離散化した以下の式と同等である。
Jp (m) = −
=
qµp (m) ψ(n) − ψ(n + 1)
qµp (m) ψ(n) − ψ(n + 1)
p(n) +
p(n + 1)
qEh
h
h
e
−1
1 − e kT
− qEh
kT
q
[λp1 (m)p(n) + λp2 (m)p(n + 1)]
h
ここで、
λp1 (m) = µp (m)
(59)
ψ(n) − ψ(n + 1)
q
1 − e− kT [ψ(n)−ψ(n+1)]
λp2 (m) = µp (m)
ψ(n) − ψ(n + 1)
q
1 − e kT [ψ(n)−ψ(n+1)]
全く同様にして、電子に対する電流密度の式も離散化することができる。その結果だけ書
くと、
q
Jn (m) = [λn1 (m)n(n) + λn2 (m)n(n + 1)]
(60)
h
ここで、
λn1 (m) = µn (m)
λn2 (m) = µn (m)
ψ(n) − ψ(n + 1)
q
1 − e kT [ψ(n)−ψ(n+1)]
ψ(n) − ψ(n + 1)
q
1 − e− kT [ψ(n)−ψ(n+1)]
うえで与えられた電流密度の式を電流連続式
1 Jp (m) − Jp (m − 1)
− G(n) + U (n) = 0
q
h
(61)
2. デバイスモデリング入門
56
1 Jn (m) − Jn (m − 1)
+ G(n) − U (n) = 0
q
h
に代入すれば、定常状態に対する離散化された電流連続方程式を得る。
一方、ポアソン方程式は以下のようになる。
q
γ1 (n)ψ(n − 1) + γ2 (n)ψ(n) + γ3 (n)ψ(n + 1) = − [Nd (n) − Na (n) + p(n) − n(n)]
ε0
(62)
(63)
ここで、（等間隔メッシュを用いているので）
γ1 (n) = γ3 (n) =
1
,
h2
γ2 (n) = −
2
h2
である。
2.2.5
ドリフト拡散方程式の線形化
うえで与えた離散化された基本方程式は強い非線形性をもつ方程式系であることから、陽解
法で解くことは収束性において困難が伴う。そこで、含まれる各未知変数 p(n), n(n), ψ(n) に
ついて、各方程式を線形化する。
電流連続式は以下のようになる。
1 Jp0 (m) − Jp0 (m − 1)
1 δJp (m) − δJp (m − 1)
− δG(n) + δU (n) = −
+ G0 (n) − U 0 (n) (64)
q
h
q
h
1 δJn (m) − δJn (m − 1)
1 J 0 (m) − Jn0 (m − 1)
+ δG(n) − δU (n) = − n
− G0 (n) + U 0 (n) (65)
q
h
q
h
δJp (m) は変分、Jp0 (m) 等はその変分を除いた部分であり、数値シミュレーションにおいては、
iterative に方程式を解くときの１回まえの値を表す（変分 δJp (m) は各 iteration での変化量に
対応する）。また、ポアソン方程式はつぎのように書かれる。
q2 0
q2 0
n (n) −
p (n) δψ(n) + γ3 (n)δψ(n + 1) =
γ1 (n)δψ(n − 1) + γ2 (n) −
ε0 kT
ε0 kT
−γ1 (n)ψ 0 (n − 1) − γ2 (n)ψ 0 (n) − γ3 (n)ψ 0 (n + 1) −
q
Nd (n) − Na (n) + p0 (n) − n0 (n) (66)
ε
ここで、正孔および電子の密度に対する変分を
q
q 0
n (n)δψ(n)
δp(n) = − p0 (n)δψ(n), , δn(n) =
kT
kT
とした。
正孔に対する電流密度 Jp (m) は補助（auxiliry）メッシュで定義されるから、Jp (m) の変分
は以下のように与えられる。
Jp (m) = Jp0 (m) + δJp (m)
∂Jp0 (m)
∂Jp0 (m)
δp(n) +
δp(n + 1)
∂p(n)
∂p(n + 1)
∂Jp0 (m)
∂Jp0 (m)
δψ(n) +
δψ(n + 1)
+
∂ψ(n)
∂ψ(n + 1)
≈ Jp0 (m) +
2. デバイスモデリング入門
57
同様に、電子に対する電流密度 Jn (m) は
Jn (m) = Jn0 (m) + δJn (m)
∂Jn0 (m)
∂J 0 (m)
δn(n) +
δn(n + 1)
≈ Jn0 (m) + n
∂n(n)
∂n(n + 1)
∂Jn0 (m)
∂J 0 (m)
δψ(n) +
δψ(n + 1)
+ n
∂ψ(n)
∂ψ(n + 1)
となる。
また、生成再結合項 U (x) は、Shockley-Reed-Hall モデルでは正孔と電子の数密度のみの関
数であるから、
U (n) = U 0 (n) + δU (n)
∂U 0 (n)
∂U 0 (n)
δp(n) +
δn(n)
≈ U 0 (n) +
∂p(n)
∂n(n)
となる。衝突イオン化による項はすべての未知変数に依存することから、単純な局所電界モデ
ルを用いたとしても結構複雑になる。結果は以下のようになる。
G(n) = G0 (n) + δG(n)
∂G0 (n)
∂G0 (n)
∂G0 (n)
δp(n − 1) +
δp(n) +
δp(n + 1)
≈ G0 (n) +
∂p(n − 1)
∂p(n)
∂p(n + 1)
∂G0 (n)
∂G0 (n)
∂G0 (n)
δn(n − 1) +
δn(n) +
δn(n + 1)
+
∂n(n − 1)
∂n(n)
∂n(n + 1)
∂G0 (n)
∂G0 (n)
∂G0 (n)
δψ(n − 1) +
δψ(n) +
δψ(n + 1)
+
∂ψ(n − 1)
∂ψ(n)
∂ψ(n + 1)
尚、うえの各変分に含まれる偏微分係数は、それぞれの定義式から容易に導くことができる。
例えば、正孔に対する電流密度 Jp (m) の正孔密度 p(n) に対する偏微分は、
∂Jp0 (m)
q
= λp1 (m)
∂p(n)
h
となる。
これらの変分の表式を線形化した離散方程式に代入すれば、ドリフト拡散法における基本方
程式となる。これらの方程式は、各未知変数についての線形方程式であり、まえと同様に３重
対角行列で表される。従って、再帰的 (recursive）な方法を用いて効率良く解くことができる。
線形化された電流連続式もポアソン方程式も以下のような線形代数方程式に書けている。
A(n)δy(n − 1) + B(n)δy(n) + C(n)δy(n + 1) = F (n),
2≤n≤N
ここで、
δy(n) = δψ(n) or δp(n) or δn(n)
である。そこで、まず ψ(n) についての線形化されたポアソン方程式をうえのような３重対角
行列の方程式に対応させ、正孔および電子密度が一定という条件のもとで解く。つぎに、正孔
2. デバイスモデリング入門
58
の電流連続式を静電ポテンシャルと電子密度が一定という条件のもとで、正孔密度 p(n) につい
ての３重対角行列の方程式を解く。そして、電子の電流連続式を静電ポテンシャルと正孔密度
が一定という条件のもとで、電子密度 n(n) についての３重対角行列の方程式として解く。これ
を各変数の変分が十分に小さくなるまで繰り返し、iterative に解く。このようにして線形化さ
れた離散方程式を順次解いて ψ(n), p(n), n(n) を求める方法が、ドリフト拡散法である。
演習問題 2-2-2
正孔および電子に対する線形化された電流連続式を、各変数の偏微分係数まで含めて導け。
演習問題 2-2-3
シミュレーション工学のホームページにポストされている pn 接合をシミュレートするための
ドリフト拡散シミュレータのプログラムを各自ダウンロードし、プログラムを解読せよ。
このプログラムを用いて、順方向と逆方向での電流電圧特性を求めよ。半導体デバイスの教
科書にのっている電流特性についての解析式の結果と比較検討せよ。
デバイス（pn 接合）の構造パラメータや物性パラメータを変化させて電流電圧特性を解析し
てみるのも興味深い。特に、pn 接合における電流電圧特性が少数キャリアの再結合によって決
まるのか、或いは少数キャリアの拡散輸送によって決まるのか、検討せよ。
3. 参考文献
3
59
参考文献
数値解析
1. R. L. Press, et. al, Numerical Recipes, 2-nd ed. (Cambridge, New York, 1992).
2. S. E. Koonin, Computational Physics (Addison-Wesley, New York, 1986).
3. M. M. Woolfson and G. J. Pert, Introduction to Computer Simulation (Oxford Univ.
Press, New York, 1999).
4. R. L. Burden, J. D. Faires, and A. C. Reynolds, Numerical Analyses, 2-nd ed. (Prindle,
Weber, and Schmidt, Boston, 1981).
5. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions (Dover, New
York, 1964).
6. G. Arfken and H. J. Weber, Mathematical Methods for Physicists, 4-th ed. (Academic,
New York, 1995).
デバイス物理およびデバイス・シミュレーション
1. Y. Taur and T. H. Ning, Fundamentals of Modern VLSI Devices (Cambridge, New York,
1998).
2. S. M. Sze, Semiconductor Devices: Physics and Technology (Wiley, New York, 1985);
S. M. Sze（南日 他訳）、半導体デバイス（産業図書）.
3. S. Serberherr, Analysis and Simulation of Semiconductor Devices (Springer, New York,
1984).
4. C. M. Snowden（浜口、谷口訳）, 半導体デバイスのモデリング (現代工学社).
5. 富沢一隆, 半導体デバイス・シミュレーション (コロナ社).
1/2 ページ
/*****************************************************************************/
/*
シミュレーション工学 課題１
*/
/*
サンプルプログラム
*/
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#define N 21
/* 刻み幅の設定の総数 */
double func(double x);
/*****************************************************************************/
/*
main関数
*/
/*
変数
用途
*/
/*
h[]
刻み幅が入る
*/
/*
fprime[]
各刻み幅に対応した２点公式による結果を保存
*/
/*****************************************************************************/
int main()
{
double x, h[N], fprime[N];
int i, j;
2/2 ページ
/*
*/
/*
2.C言語において,
*/
/*
a/b*cは、 a/bを計算してから、その結果にcを掛けることを意味する。決し
*/
/*
て、a/bを計算してからcで割るのではないことに注意せよ。後者を計算する
*/
/*
時は、
*/
/*
a/b/c あるいは a/(b*c)
*/
/*
とする必要がある。
*/
/*
*/
/*
3.結果をprintf()で出力する際に、
*/
/*
printf("%.10e¥n", fx);
*/
/*
^^^
*/
/*
とすれば、小数点以下10桁まで表示できる。%fなどについても同様である。
*/
/*
*/
/*****************************************************************************/
x = 1.5;
/* 刻み幅の設定 */
/* h[0] = 0.5, h[1] = 0.2, h[3] = 0.1 .... */
j = 0;
for (i = 0; i < N; i = i + 3) {
h[i+2] = 0.1/pow(10.0, j);
h[i+1] = 2.0*h[i+2];
h[i]
= 5.0*h[i+2];
j = j + 1;
}
/* ２点公式(前進) */
for (i = 0; i < N; i++)
fprime[i] = (func(x+h[i]) - func(x))/h[i];
/* ２点公式(前進)による計算結果を出力する */
printf("%-9s¥t%-12s¥n", "h", "fprime");
for (i = 0; i < N; i++)
printf("%.7f¥t%.10f¥n", h[i], fprime[i]);
return 0;
}
/*****************************************************************************/
/*
関数f(x) = x^(x^x)
*/
/*****************************************************************************/
double func(double x)
{
return pow(x, pow(x,x));
}
/*****************************************************************************/
/*
プログラムにおける注意点
*/
/*
*/
/*
1.冪乗の計算には,Cのライブラリ関数pow()を使用した。関数プロトタイプは */
/*
double pow(double x, double y);
*/
/*
例えば、f(x) = x^yを計算して、その結果を変数fxに格納するには、
*/
/*
fx = pow(x, y);
*/
/*
とすればよい。なお、Cの数学関数ライブラリを使用する際には、必ずプログ */
/*
ラムのヘッダにおいてmath.hをインクルードすること。また、コンパイル
*/
/*
には、以下のように -lmをつけてコンパイルする必要がある。
*/
/*
% cc -o report1-1 report101.c -lm
*/
http://hermes.esys.tsukuba.ac.jp/~sano/cprog/ndiff.c
http://hermes.esys.tsukuba.ac.jp/~sano/cprog/ndiff.c
simp.txt
/*****************************************************************************/
/*
シミュレーション工学
*/
/*
Sample program (Simpson Method)
*/
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
simp.txt
/*****************************************************************************/
/*
External function
*/
/*****************************************************************************/
double gaussf(double x)
{
double gx;
#define DATA 10 /* number of mesh */
/*#define M_PI 3.1415926535897932385 */
gx = exp(-x*x);
return gx;
double gaussf(double x);
/*****************************************************************************/
/*
main
*/
/*****************************************************************************/
int main()
{
int i, j, n[DATA];
double exactf, h[DATA];
double fx, x, sum, simpson[DATA];
double a, b, fac;
n[0] = 4;
for(i = 1; i < DATA; i++) {
n[i] = 2*n[i-1];
}
exactf = sqrt(M_PI)*erf(1.0)/2.0;
printf("¥n");
printf(" exact of integral: %.15f¥n", exactf);
a = 0.0;
b = 1.0;
for(i =0; i < DATA; i++){
h[i] = (b-a)/n[i];
sum = gaussf(a);
fac = 2.0;
for(j = 1; j < n[i]; j++) {
x = j*h[i]+a;
fx = gaussf(x);
if (fac == 2.0) {
fac = 4.0;
}
else {
fac = 2.0;
}
sum = sum+fac*fx;
}
sum = sum+gaussf(b);
sum = sum*h[i]/3.0;
simpson[i] = sum;
}
}
exact
n
4
8
16
32
64
128
256
512
1024
2048
of integral: 0.746824132812427
h
simpson
0.2500000
0.7468553797909873
0.1250000
0.7468261205274666
0.0625000
0.7468242574357303
0.0312500
0.7468241406069850
0.0156250
0.7468241332996727
0.0078125
0.7468241328428814
0.0039062
0.7468241328143308
0.0019531
0.7468241328125457
0.0009766
0.7468241328124338
0.0004883
0.7468241328124282
n
4
8
16
32
64
128
256
512
1024
2048
h
0.2500000
0.1250000
0.0625000
0.0312500
0.0156250
0.0078125
0.0039062
0.0019531
0.0009766
0.0004883
error (simpson)
-0.0000312469785603
-0.0000019877150396
-0.0000001246233033
-0.0000000077945580
-0.0000000004872457
-0.0000000000304544
-0.0000000000019038
-0.0000000000001187
-0.0000000000000068
-0.0000000000000012
printf("%6s¥t%-9s¥t%-18s¥n", "n", "h", "simpson");
for(i = 0; i < DATA; i++) {
printf("%6d¥t%.7f¥t%.16f¥n", n[i], h[i], simpson[i]);
}
printf("¥n");
printf("%6s¥t%-9s¥t%-18s¥n", "n", "h", "error (simpson)");
for(i = 0; i < DATA; i++) {
printf("%6d¥t%.7f¥t%.16f¥n", n[i], h[i], exactf-simpson[i]);
}
return 0;
}
ページ(1)
ページ(2)
1/2 ページ
/*****************************************************************************/
/*
シミュレーション工学
*/
/*
Sample program (２分法)
*/
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
2/2 ページ
}
#define N 50 /* Max number of iterations */
#define exactx 1.108575508057246 /* Exact value of x */
double func(double x);
double dfunc(double x);
/*****************************************************************************/
/*
main
*/
/*****************************************************************************/
int main()
{
double x, tolx, dx, fold, fnew;
double ans[N], error[N];
int i, j;
for (i = 0; i < N; i++) {
ans[i]
= 0.0;
error[i]
= 0.0;
}
x = -1.0;
tolx = 1.0e-6;
dx = 0.5;
i = 0;
while (fabs(dx) > tolx) {
fold = func(x);
x = x+dx;
ans[i] = x;
error[i] = exactx-ans[i];
fnew = func(x);
if(fnew*fold < 0.0)
dx = -dx/2.0;
i = i+1;
if(i > N) {
printf("Does not converge...¥n");
exit(1);
}
}
printf("%s%5d¥n","
number of iterations: ",i);
printf("%5s¥t%16s¥t%16s¥n", "it","x","error");
for (j = 0; j < i; j++) {
printf("%5d¥t%.14f¥t%.14f¥n", j, ans[j], error[j]);
}
return 0;
}
/*****************************************************************************/
/*
External functions
*/
/*****************************************************************************/
double func(double x)
{
return pow(x,5.0)+3.0*x-5.0;
}
double dfunc(double x)
{
return 5.0*pow(x,4.)+3.0;
http://hermes.esys.tsukuba.ac.jp/~sano/cprog/nibun.c
http://hermes.esys.tsukuba.ac.jp/~sano/cprog/nibun.c
newton.c
/*****************************************************************************/
/*
シミュレーション工学
*/
/*
Sample program (Newton Method)
*/
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#define N 20 /* Max number of iterations */
#define exactx 1.108575508057246 /* Exact value of x */
double func(double x);
double dfunc(double x);
/*****************************************************************************/
/*
main
*/
/*****************************************************************************/
int main()
{
double x, tolx, dx, fx, dfx;
double ans[N], error[N];
int i, j;
for (i = 0; i < N; i++) {
ans[i]
= 0.0;
error[i] = 0.0;
}
x = 10.0;
tolx = 1.0e-6;
dx = 0.5;
i = 0;
while (fabs(dx) > tolx) {
fx = func(x);
dfx = dfunc(x);
dx = -fx/dfx;
x = x+dx;
ans[i] = x;
error[i] = exactx-ans[i];
i = i+1;
if(i > N) {
printf("Does not converge...¥n");
exit(1);
}
}
printf("%s%5d¥n"," number of iterations: ",i);
printf("%5s¥t%16s¥t%16s¥n", "it","x","error");
for (j = 0; j < i; j++) {
printf("%5d¥t%.14f¥t%.14f¥n", j+1, ans[j], error[j]);
}
return 0;
}
/*****************************************************************************/
/*
External functions
*/
/*****************************************************************************/
double func(double x)
{
return pow(x,5.0)+3.0*x-5.0;
}
double dfunc(double x)
{
return 5.0*pow(x,4.)+3.0;
return pow(x, pow(x,x));
}
ページ(1)
euler.c
/*****************************************************************************/
/*
Simulation Analysis
*/
/*
Sample program (Euler's Method)
*/
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#define N 30 /* Number of steps */
double func(double x, double y);
/*****************************************************************************/
/*
main
*/
/*****************************************************************************/
int main()
{
double x, x0, x1, xmin, xmax, h, error;
double y, y0, y1;
int i, j, k;
xmin = 0.0;
xmax = 6.0;
h = fabs(xmax-xmin)/N;
x0 = xmin;
y0 = 1.0;
x1 = x0+h;
y1 = 1.0-h+h*h/2.0-h*h*h/6.0;
y1 = 1.0-h+h*h/2.0;
printf("%s%5d¥n"," number of steps ",N);
for (i = 2; i <= N; i++) {
x = x1+h;
y = y0+2.0*h*func(x1,y1);
error = y-exp(-x);
printf("%5d¥t%10f%12.6f%12.6f%12.6f¥n",i, x, exp(-x),y, error);
x0 = x1;
y0 = y1;
x1 = x;
y1 = y;
}
return 0;
}
/*****************************************************************************/
/*
External functions
*/
/*****************************************************************************/
double func(double x, double y)
{
return -y;
}
ページ(1)
1/2 ページ
2/2 ページ
/*****************************************************************************/
/*
Bohr's Quantization Rule
*/
/*
Sample program (Direct Method)
*/
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
int i, ns;
#define N 50 /* Max number of iterations */
#define m 1
/* Energy quantum level */
#define Pi 3.1415926535897932385
/* Pi */
x = xmin;
fac = 2.0;
sum = func(en,x);
for (i = 1; i < ns; i++) {
x = x+dx;
if (fac == 2.0) {
fac = 4.0;
}
else {
fac = 2.0;
}
sum = sum+fac*func(en,x);
}
x = xmax;
sum = (sum+func(en,x))*dx/3.0;
xmin
xmax
ns =
dx =
double dsimp(double en);
double func(double en, double x);
/*****************************************************************************/
/*
main
*/
/*****************************************************************************/
int main()
{
double en, tole, de, fold, fnew;
double exacte;
double ans[N], error[N];
int i, j;
for (i = 0; i < N; i++) {
ans[i]
= 0.0;
error[i]
= 0.0;
}
en = 0.1;
tole = 1.0e-6;
de = 0.5;
i = 0;
fold = dsimp(en)-m;
while (fabs(de) > tole) {
en = en+de;
ans[i] = en;
exacte = m*m;
error[i] = exacte-ans[i];
fnew = dsimp(en)-m;
if(fnew*fold < 0.0)
de = -de/2.0;
i = i+1;
fold = fnew;
if(i > N) {
printf("Does not converge...¥n");
exit(1);
}
}
/*
= 0.0;
= 1.0;
100;
fabs(xmax-xmin)/ns;
printf("%s%5d¥t%.14f¥n","
number of mesh: ",i,sum); */
return sum;
}
/*****************************************************************************/
/*
External functions (Integrand)
*/
/*****************************************************************************/
double func(double en, double x)
{
return sqrt(en);
}
printf("%s%5d¥n","
number of iterations: ",i);
printf("%5s¥t%15s¥t%16s¥n", "it","en","error");
for (j = 0; j < i; j++) {
printf("%5d¥t%.14f¥t%.14f¥n", j, ans[j], error[j]);
}
return 0;
}
/*****************************************************************************/
/*
External functions
(Simpson's method)
*/
/*****************************************************************************/
double dsimp(double en)
{
double x, xmin, xmax, dx;
double fac, sum;
http://hermes.esys.tsukuba.ac.jp/~sano/cprog/bohr.c
http://hermes.esys.tsukuba.ac.jp/~sano/cprog/bohr.c
bohros.c
/*****************************************************************************/
/*
Bohr's Quantization Rule
*/
/*
Sample program (Harmonic Oscillator)
*/
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
bohros.c
double fac, sum;
int i, ns;
xmin = 0.0;
xmax = sqrt(en);
ns = 100;
dx = fabs(xmax-xmin)/ns;
#define N 50
/* Max number of iterations */
#define m 0
/* Energy quantum level m = 0,1,2,... */
#define Pi 3.1415926535897932385 /* pi */
x = xmin;
fac = 2.0;
sum = func(en,x);
for (i = 1; i < ns; i++) {
x = x+dx;
if (fac == 2.0) {
fac = 4.0;
}
else {
fac = 2.0;
}
sum = sum+fac*func(en,x);
}
x = xmax;
sum = (sum+func(en,x))*dx/3.0;
double dsimp(double en);
double func(double en, double x);
/*****************************************************************************/
/*
main
*/
/*****************************************************************************/
int main()
{
double en, tole, de, fold, fnew;
double exacte;
double ans[N], error[N];
int i, j;
for (i = 0; i < N; i++) {
ans[i]
= 0.0;
error[i] = 0.0;
}
en = 0.1;
tole = 1.0e-6;
de = 0.1;
i = 0;
fold = dsimp(en)-m-0.5;
while (fabs(de) > tole) {
en = en+de;
ans[i] = en;
exacte = m+0.5;
error[i] = exacte-ans[i];
fnew = dsimp(en)-m-0.5;
if(fnew*fold < 0.0)
de = -de/2.0;
i = i+1;
fold = fnew;
if(i > N) {
printf("Does not converge...¥n");
exit(1);
/*return 0;*/
}
}
/* printf("%s%5d¥t%.14f¥t%.6f¥n"," number of mesh: ",i,sum,en);*/
return sum;
}
/*****************************************************************************/
/*
External functions (Integrand)
*/
/*****************************************************************************/
double func(double en, double x)
{
double fx;
fx = en-x*x;
if (fx < 0.0)
fx = 0.0;
fx = 4.0/Pi*sqrt(fx);
return fx;
}
printf("%s%5d¥n"," number of iterations: ",i);
printf("%5s¥t%15s¥t%16s¥n", "it","en","error");
for (j = 0; j < i; j++) {
printf("%5d¥t%.14f¥t%.14f¥n", j, ans[j], error[j]);
}
/* ans[i]=func(Pi,0.25);*/
/*printf("¥n%5d¥t%.14f¥t%.14f¥n",i,ans[i],Pi);*/
return 0;
}
/*****************************************************************************/
/*
External functions (Simpson's method)
*/
/*****************************************************************************/
double dsimp(double en)
{
double x, xmin, xmax, dx;
ページ(1)
ページ(2)
moscd.c
/*****************************************************************************/
/*
Simulation Analysis : MOS-Capacitor
*/
/*
Sample program (Euler's Method)
*/
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define
#define
#define
#define
#define
moscd.c
return fk;
}
temp 300.0
/* Temperature */
q 1.60210e-19 /* elementary charge */
bk 1.38054e-23 /* Boltzmann constant */
eps0 8.854185e-12 /* vacuum dielectric constant */
pi 3.1415926535897932385 /* pi */
double func(double ac, double y);
/*****************************************************************************/
/*
main (MKS-units)
*/
/*****************************************************************************/
int main()
{
double cni, eps, ac, fct, vt, Qs, h;
double psis, psi0, psi1, psi2;
double x, y1;
int nmax, i;
cni = 1.05e16;
eps = 11.9*eps0;
ac = 1.e24;
vt = bk*temp/q;
ac = ac/cni;
fct = q*cni/eps;
psis = 0.4;
Qs = 3.4743e-3;
//Qs = 2.0e-3;
h = 1.0e-9;
nmax = (int)(50.0e-9/h+1.e-5);
printf("%s%5d¥t%.25e¥n"," number of steps ",nmax,50.0e-9/h);
printf("%s¥t%.10e¥n"," charge density Qs (C/cm^2) ",Qs/1.0e4);
psi0 = psis;
psi1 = psis-h*Qs/eps;
printf("¥n%5s%14s¥t%13s¥n", "n","x (nm)","psi (V)");
for (i = 2; i <= nmax; i++) {
x = i*h;
y1 = psi1/vt;
psi2 = 2.0*psi1-psi0-h*h*fct*func(ac,y1);
if(fabs(psi2) > 1.0e10) {
printf(" divergent...¥n");
exit(1);
return 0;
}
printf("%5d¥t%10.3f¥t%15e¥n", i, x/1.0e-9, psi2);
psi0 = psi1;
psi1 = psi2;
}
return 0;
}
/*****************************************************************************/
/*
External functions
*/
/*****************************************************************************/
double func(double ac, double y)
{
double fk;
fk = ac*(exp(-y)-1.0)-1.0/ac*(exp(y)-1.0);
ページ(1)
ページ(2)
poisson.c
/*****************************************************************************/
/*
Simulation Analysis : MOS-Capacitor
*/
/*
Sample program (Implicit Method)
*/
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define
#define
#define
#define
#define
#define
poisson.c
psi[i] = psi[i]+up[i-2];
/*printf("%5d¥t%15e¥n", i, up[i-2]);*/
if (umax < fabs(up[i-2])) {
umax = fabs(up[i-2]);
}
}
deltp = umax;
if (it >= Nmax) {
printf("Does not converge...¥n");
exit(1);
}
temp 300.0
/* Temperature */
q 1.60210e-19 /* elementary charge */
bk 1.38054e-23 /* Boltzmann constant */
eps0 8.854185e-12 /* vacuum dielectric constant */
pi 3.1415926535897932385 /* pi */
Nmax 50
/* Max number of iteration */
double func(double ac, double y);
void tridag(double *a,double *b,double *c,double *u,double *f,int n);
/*****************************************************************************/
/*
main (MKS-units)
*/
/*****************************************************************************/
int main()
{
double cni, eps, ac, fct, vt, h, tolx;
double psis, deltp, Qs;
double xmax, x, y1, umax;
double psi[1000], ap[1000], bp[1000], cp[1000], up[1000], rp[1000];
int nmax, i, it;
psis = 0.4;
cni = 1.05e16;
eps = 11.9*eps0;
ac = 1.e24;
xmax = 100.0e-9;
tolx = 1.0e-5;
vt = bk*temp/q;
ac = ac/cni;
h = 1.0e-9;
nmax = xmax/h+1.e-5;
if (nmax > 1000) {
printf("enlarge the size of arrays...¥n");
exit(1);
}
fct = h*h*q*cni/eps;
}
/*summary of calculation results */
Qs = -eps*(psi[2]-psi[1])/h;
printf("¥n%5d¥t%s¥t%.10e¥n¥n",
it," charge density Qs (C/cm^2) ",Qs/1.0e4);
for (i = 1; i <= nmax+1; i++) {
x = (i-1)*h;
printf("%5d¥t%10.3f¥t%15e¥n", i, x/1.0e-9, psi[i]);
}
return 0;
}
/*****************************************************************************/
/*
External functions
*/
/*****************************************************************************/
double func(double ac, double y)
{
double fk;
fk = ac*(exp(-y)-1.0)-1.0/ac*(exp(y)-1.0);
return fk;
}
/*****************************************************************************/
void tridag(double *a,double *b,double *c,double *u,double *f,int n)
{
double *gam, bet;
int j;
/* memory size for gam */
gam = (double *)malloc(sizeof(double)*n);
if (b[0] == 0.0) {
fprintf(stderr,"determinant in tridag() 1¥n");
exit(1);
}
printf("%s%5d¥t%.25e¥n"," number of steps ",nmax,xmax/h);
/*** initial configuration of Potential ***/
for (i = 1; i <= nmax+1; i++) {
psi[i] = psis;
if (i > nmax/2) psi[i] = 0.0;
/*printf("%5d¥t%10.3f¥t%15e¥n", i, (i-1)*h/1.0e-9, psi[i]);*/
}
it = 0;
deltp = 1.0e20;
while (deltp > tolx) {
for (i = 2; i <= nmax; i++) {
y1 = psi[i]/vt;
ap[i-2] = 1.0;
bp[i-2] = -2.0-fct/vt*(ac*exp(-y1)+1.0/ac*exp(y1));
cp[i-2] = 1.0;
rp[i-2] = -psi[i-1]+2.0*psi[i]-psi[i+1]-fct*func(ac,y1);
}
it = it+1;
tridag(ap, bp, cp, up, rp, nmax-1);
umax = 1.0e-20;
for (i = 2; i <= nmax; i++) {
ページ(1)
bet = b[0];
u[0] = f[0]/bet;
for (j = 1; j < n; j++) {
gam[j] = c[j-1]/bet;
bet = b[j] - a[j]*gam[j];
if (bet == 0.0) {
fprintf(stderr,"determinant in tridag() 2 ...¥n");
exit(1);
}
u[j] = (f[j] - a[j]*u[j-1])/bet;
}
for (j = n-2; j >= 0; j--) {
u[j] = u[j] - gam[j+1]*u[j+1];
}
free(gam);
}
ページ(2)
simdd.c
/*****************************************************************************/
/*
Dift-Diffusion Simulator (1-D)
*/
/*
(Decoupled Method)
*/
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define
#define
#define
#define
#define
#define
ECHARGE 1.602189e-19
/* 電子の素電荷 [c]
*/
EMASS 9.109534e-31
/* 電子の質量
[kg] */
BOLTZ 1.380662e-23
/* ボルツマン定数 [J/K] */
PI
3.141593e+00
/* π
*/
EPS0 8.854188e-12
/* 真空の誘電率 [F/m] */
NM 10000
/* 格子点の最大数
*/
int ddmain(double vmax, double vdlt);
void tdmtrx(int ndm, double *at, double *bt, double *ct, double *ft, int ll2);
void tridag(double *a, double *b, double *c, double *r, double *u, int n);
int poisson(double *at,double *bt,double *ct,double *ft,double *ut);
int econt(double *at,double *bt,double *ct,double *ft,double *ut);
int hcont(double *at,double *bt,double *ct,double *ft,double *ut);
void clear(double *a, int n);
void recomb(int ndm, int ll, double *rcmb);
void mobility(int ndm, int ll, double *doud);
void setting(double *acc,double *don,double *volt,double *cgm,double *elec,double *hole);
void smpmob(int ndm, int ll, double *doud);
double brn(double x);
void output(double *, double *, double *);
int ll, ll2, lnd, lpd, ltot, itmax, ibmax;
double dltex;
double rmass, eps, cdin, taup, taun, vt, temp;
double vmax, vdlt, vold;
double volt[NM], elec[NM], hole[NM];
double don[NM], acc[NM], cgm[NM];
double ecurrt[NM], hcurrt[NM];
double fl[3];
double dl1[NM], dl2[NM], dudn[NM];
double dp1[NM], dp2[NM], dudp[NM];
double gam1[NM], gam2[NM], gam3[NM], hsp[NM], hprm[NM];
/*****************************************************************************/
/*
drift diffusion main routine
*/
/*****************************************************************************/
int main()
{
extern double rmass, eps, cdin, taup, taun;
extern double don[NM], acc[NM], cgm[NM];
extern double volt[NM], elec[NM], hole[NM];
extern double vmax, vdlt;
/* Siのパラメータ */
rmass = 0.328*EMASS; /* 有効質量
*/
eps = 11.9*EPS0;
/* 誘電率
*/
cdin = 1.105e10*1e6; /* 真性キャリア密度 [m^3] */
taup = 10.0e-9;
/* 正孔のライフタイム [s] */
taun = 10.0e-9;
/* 電子のライフタイム [s] */
/* 格子点、不純物プロファイル、温度、初期解の設定 */
setting(acc, don, cgm, volt, elec, hole);
/* 許容誤差、最大反復回数の設定 */
itmax = 1000;
fl[0] = 1.0e-5; /* ポアソン方程式 */
fl[1] = 1.0e-5; /* 電子電流連続の式 */
fl[2] = 1.0e-5; /* 正孔電流連続の式 */
/* 印加電圧の設定 */
ページ(1)
simdd.c
vmax = 0.5; /* 最大電圧 */
vdlt = 0.1; /* 刻み幅 */
/* ドリフト拡散メインルーチン */
ddmain(vmax, vdlt);
/* シミュレーション結果の解析 */
output(volt,elec,hole);
return 0;
}
/*****************************************************************************/
/*
格子点、メッシュサイズ、不純物プロファイルの設定
*/
/*****************************************************************************/
void setting(double *acc,double *don,double *cgm,double *volt,double *elec,double *hole)
{
extern int ll, ll2, lpd, lnd, ltot;
extern double hsp[NM], hprm[NM], gam1[NM],gam2[NM],gam3[NM];
extern double temp, vold;
double pdope, ndope;
int i;
printf("Device structure pn-diode¥n");
/* 格子点の数の設定 */
lpd = 100;
/* p領域 */
lnd = 100;
/* n領域 */
ltot = lpd + lnd;
ll = ltot + 1; /* 全格子点数 */
ll2 = ll - 2;
/* 境界の点を除いた格子点数 */
/* メッシュ幅 [m] */
dltex = 10.0e-9;
/* 温度 */
temp = 300.0;
/* 一様不純物密度 */
pdope = 1.0e16;
/* [cm^-3] */
ndope = 1.0e16;
/* [cm^-3] */
/* ドーピングプロファイルの設定 */
pdope = pdope*1.0e6;
ndope = ndope*1.0e6;
clear(don,NM); clear(acc,NM);
for (i = 0; i < ll; i++) {
if (lpd <= i && i < ll) {
/* n-region */
don[i] = ndope;
acc[i] = 0.0;
}
else {
/* p-region*/
don[i] = 0.0;
acc[i] = pdope;
}
cgm[i] = don[i] - acc[i];
}
/* 雑多な計算 */
vt = ECHARGE/(BOLTZ*temp);
for (i = 0; i < ll; i++) {
hsp[i] = dltex;
hprm[i] = (dltex+dltex)/2.0;
}
clear(gam1,NM);
clear(gam2,NM);
clear(gam3,NM);
for (i = 1; i < ll-1; i++) {
ページ(2)
simdd.c
simdd.c
hcurrt[ib] = hcr = ECHARGE/hprm[lj]
*(dp1[lj]*hole[lj] + dp2[lj]*hole[lj+1]);
sumcur = ecr + hcr;
built = volt[ll-1] - volt[0];
/* printf("built-in: %9.5f ¥n", volt[ll-1]-volt[0]); */
printf("%5d %9.3f %10.5f %15.5E %15.5E %15.5E¥n",ib+1,ib*vdlt, built,
ecr/1.e4, hcr/1.e4, (ecr+hcr)/1.e4);
vnew = vnew + vdlt;
/* ends loop! */
gam1[i] = 1.0/hsp[i-1]/hprm[i];
gam3[i] = 1.0/hsp[i]/hprm[i];
gam2[i] = -gam1[i] - gam3[i];
}
/* 初期解の設定 */
for (i = 0; i < ll; i++) {
/* 電位 */
volt[i] = 0.0;
if (cgm[i] < 0.0)
volt[i] = log(-cdin/cgm[i])/vt;
if(cgm[i] > 0.0)
volt[i] = log(cgm[i]/cdin)/vt;
vold = volt[0];
/* 電子 */
elec[i] = don[i];
/* 正孔 */
hole[i] = acc[i];
}
hole[0] = -cgm[0]*(1.0+sqrt(1.0+
(2.0*cdin/cgm[0])*(2.0*cdin/cgm[0])) )/2.0;
elec[0] = cdin*cdin/hole[0];
elec[ll-1] = cgm[ll-1]*(1.0+sqrt(1.0+
(2.0*cdin/cgm[ll-1])*(2.0*cdin/cgm[ll-1])) )/2.0;
hole[ll-1] = cdin*cdin/elec[ll-1];
}
return 0;
}
/*****************************************************************************/
/*
ポアソン方程式を解く
*/
/*****************************************************************************/
int poisson(double *at,double *bt,double *ct,double *ft,double *ut)
{
extern int itmax;
extern double volt[NM], elec[NM], hole[NM], vt, fl[3];
double error, max_error;
int it, lj, ndm;
}
/****************************************************************************/
/*
ドリフト拡散メインルーチン
*/
/****************************************************************************/
int ddmain(double vmax, double vdlt)
{
extern int ll, ibmax;
extern double vold;
extern double volt[NM], elec[NM], hole[NM], hprm[NM];
int ib, it, lj, itrt[3];
double built, vnew, sumcur, ecr, hcr;
double at[NM], bt[NM], ct[NM], ft[NM], ut[NM]; /* 三重対角行列 */
/* 印加電圧の設定 */
if (fabs(vdlt) > fabs(vmax)) {
if(vmax != 0.0) {
printf("Too large V-step (vdlt)¥n");
exit(1);
}
}
if (vmax < 0.0)
vdlt = -vdlt;
ibmax = (int)(fabs(vmax/vdlt+1.0e-5)) + 1;
printf("bias step: %d¥n", ibmax);
vnew = vold+0.0;
for (ib = 0; ib < ibmax; ib++) {
/* loop for bias */
volt[0] = vnew;
for (it = 0; it < itmax; it++) {
itrt[0] = poisson(at, bt, ct, ft, ut);
itrt[1] = econt(at, bt, ct, ft, ut);
itrt[2] = hcont(at, bt, ct, ft, ut);
if (itrt[0] == 1 && itrt[1] == 1 && itrt[2] == 1)
goto converged;
}
fprintf(stderr, "TOO MANY IT: NOT CONVERGED¥n");
exit(1);
converged:
lj = ll/2;
ecurrt[ib] = ecr = ECHARGE/hprm[lj]
*(dl1[lj]*elec[lj] + dl2[lj]*elec[lj+1]);
ndm = 0; /* poisson equation */
for (it = 1; it < itmax; it++) {
tdmtrx(ndm, at, bt, ct, ft, ll2);
tridag(at,bt,ct,ft,ut,ll2);
max_error = -1.0;
for (lj = 1; lj < ll-1; lj++) {
/* new potential */
volt[lj] = volt[lj] + ut[lj-1];
elec[lj] = elec[lj]*exp(vt*ut[lj-1]);
hole[lj] = hole[lj]*exp(-vt*ut[lj-1]);
/* search max relative error */
if (fabs(volt[lj]) == 0.0) {
fprintf(stderr, "???¥n");
exit(1);
}
else {
error = fabs(ut[lj-1]/volt[lj]);
if (error > max_error)
max_error = error;
/* end error searcn */
}
}
/* printf("dmax %d %e¥n", it, max_error); */
if (max_error < fl[0]) {
return it;
}
}
/* ここまできたら収束せず */
fprintf (stderr,"TOO MANY IT; NOT CONVERGED IN POISSON SOLVER¥n");
exit(1);
return 0;
}
/*****************************************************************************/
/*
電子電流連続の式を解く
*/
/*****************************************************************************/
int econt(double *at,double *bt,double *ct,double *ft,double *ut)
{
extern int itmax, ll2;
extern double elec[NM];
double error, max_error;
int it, lj, ndm;
ndm = 1;
for (it = 1; it < itmax; it++) {
ページ(3)
ページ(4)
simdd.c
tdmtrx(ndm, at, bt, ct, ft, ll2);
tridag(at,bt,ct,ft,ut,ll2);
max_error = -1.0;
for (lj = 1; lj < ll-1; lj++) {
/* renew lectron density */
elec[lj] = elec[lj] + ut[lj-1];
/* search max relative error */
if (fabs(elec[lj]) == 0.0) {
fprintf(stderr, "???¥n");
exit(1);
}
else {
error = fabs(ut[lj-1]/elec[lj]);
if (error > max_error)
max_error = error;
/* end error searcn */
}
}
/* printf("dmax %d %e¥n", it, max_error); */
if (max_error < fl[1]) {
return it;
}
}
/* ここまできたら収束せず */
fprintf (stderr,"TOO MANY IT; NOT CONVERGED IN ECONT SOLVER¥n");
exit(1);
return 0;
}
/*****************************************************************************/
/*
正孔電流連続の式を解く
*/
/*****************************************************************************/
int hcont(double *at,double *bt,double *ct,double *ft,double *ut)
{
extern int itmax;
extern double hole[NM];
double error, max_error;
int it, lj, ndm;
ndm = 2;
for (it = 1; it < itmax; it++) {
tdmtrx(ndm, at, bt, ct, ft, ll2);
tridag(at,bt,ct,ft,ut,ll2);
max_error = -1.0;
for (lj = 1; lj < ll-1; lj++) {
/* new potential, eletron density, hole density... */
hole[lj] = hole[lj] + ut[lj-1];
/* search max relative error */
if (fabs(hole[lj]) == 0.0) {
fprintf(stderr, "???¥n");
exit(1);
}
else {
error = fabs(ut[lj-1]/hole[lj]);
if (error > max_error)
max_error = error;
/* end error searcn */
}
}
/* printf("dmax %d %e¥n", it, max_error); */
if (max_error < fl[2]) {
return it;
}
}
/* ここまできたら収束せず */
fprintf (stderr,"TOO MANY IT: NOT CONVERGED IN HCONT SOLVER¥n");
exit(1);
return 0;
ページ(5)
simdd.c
}
/*****************************************************************************/
/*
結果を出力
*/
/*****************************************************************************/
void output(double *volt,double *elec,double *hole)
{
extern double dltex, dl1[NM],dl2[NM], hprm[NM];
extern double ecurrt[NM], hcurrt[NM];
extern double don[NM], acc[NM], cgm[NM];
extern int ibmax;
FILE *fp0, *fp1, *fp2, *fp3;
double cjn, cjp, dvl_el, dvl_hl, field[NM];
double qro, sarea, cr1, cr2, ecr, hcr;
int ib, lj;
fp0
fp1
fp2
fp3
=
=
=
=
fopen("cc17.txt",
fopen("cc18.txt",
fopen("cc19.txt",
fopen("cc20.txt",
"w");
"w");
"w");
"w");
for (lj = 0; lj < ll; lj++) {
qro = don[lj]-acc[lj]+hole[lj]-elec[lj];
fprintf(fp0,"%8.1f %14.5e %14.5e %14.5e %14.5e ¥n", lj*dltex*1.0e9,
qro/1.e6, hole[lj]/1.e6, elec[lj]/1.e6, cgm[lj]/1.e6);
}
clear(field, NM);
for (lj = 0; lj < ll; lj++) {
cjn = cjp = 0.0;
if (lj == 0)
field[lj] = 0.0;
if (lj != 0 && lj < ll-1) {
field[lj] = (volt[lj-1] - volt[lj+1])/hsp[lj]/2.0;
cjn = ECHARGE/hprm[lj]*(dl1[lj]*elec[lj]+dl2[lj]*elec[lj+1]);
cjp = ECHARGE/hprm[lj]*(dp1[lj]*hole[lj]+dp2[lj]*hole[lj+1]);
}
if (lj == ll-1)
field[lj] = field[ll-1];
fprintf(fp1, "%8.1f %9.5f ", lj*dltex*1.0e9, volt[lj]);
fprintf(fp1, "%13.5E %13.5E ", field[lj]/1.0e5, cjn/1.0e4);
fprintf(fp1, "%13.5E %13.5E¥n",cjp/1.0e4, (cjn+cjp)/1.0e4);
dvl_el = 0.0;
if (elec[lj] != 0.0)
dvl_el = cjn/elec[lj]/ECHARGE;
dvl_hl = 0.0;
if (hole[lj] != 0.0)
dvl_hl = cjp/hole[lj]/ECHARGE;
fprintf(fp2, "%8.1f %10.5f ", lj*dltex*1e9, volt[lj]);
fprintf(fp2, "%13.5E %13.5E ¥n", elec[lj]/1e6,hole[lj]/1e6);
//fprintf(fp2, "%13.5E %13.5E¥n", dvl_el*1e2,dvl_hl*1e2);
}
/* at each step */
sarea = 1.0;
cr1 = cr2 = ecr = hcr = 0.0;
for (ib = 0; ib < ibmax; ib++) {
if (ib >= 0) {
cr1 = ecurrt[ib];
cr2 = hcurrt[ib];
ecr = sarea*cr1/1e4;
hcr = sarea*cr2/1e4;
}
fprintf(fp3,"%5d %9.3f %15.5E %15.5E %15.5E¥n",ib+1, ib*vdlt, ecr, hcr, ecr+hcr);
}
ページ(6)
simdd.c
fclose(fp0); fclose(fp1); fclose(fp2); fclose(fp3);
}
/*****************************************************************************/
/*
配列の全要素を0にする
*/
/*****************************************************************************/
void clear(double *a, int n)
{
int i;
for (i = 0; i < n; i++) {
a[i] = 0.0;
}
}
/*****************************************************************************/
/*
三重対角行列の決定(定常状態)
*/
/*
ndm == 0: ポアソン方程式
*/
/*
ndm == 1: 電子電流連続の式
*/
/*
ndm == 2: 正孔電流連続の式
*/
/*****************************************************************************/
void tdmtrx(int ndm, double *at, double *bt, double *ct, double *ft, int ll2)
{
int mt, nt;
extern double vt;
double beta[NM], doud[NM], rcmb[NM];
double dnmnt;
/* for Poisson's equation */
if (ndm != 0){
smpmob(ndm, ll2+2, doud); /* mobility
*/
recomb(ndm, ll2+2, rcmb); /* recombination */
clear(beta,NM);
for (mt = 0; mt < ll2+1; mt++) {
beta[mt] = vt*(volt[mt] - volt[mt+1]);
}
}
/* for electron continuity equation */
if (ndm == 1) {
clear(dl1,NM);
clear(dl2,NM);
clear(dudn,NM);
for (mt = 0; mt < ll2+1; mt++) {
dnmnt = taun*(hole[mt] + cdin) + taup*(elec[mt] + cdin);
dudn[mt] = (elec[mt]-taup*rcmb[mt])/dnmnt;
dl1[mt] = doud[mt]*(-brn(beta[mt]))/vt;
dl2[mt] = doud[mt]*brn(-beta[mt])/vt;
}
}
/* for hole continuity equation */
if (ndm == 2) {
clear(dp1,NM);
clear(dp2,NM);
clear(dudp,NM);
for (mt = 0; mt < ll2+1; mt++) {
dnmnt = taun*(hole[mt] + cdin) + taup*(elec[mt] + cdin);
dudp[mt] = (elec[mt] - taun*rcmb[mt])/dnmnt;
dp1[mt] = doud[mt]*brn(-beta[mt])/vt;
dp2[mt] = doud[mt]*(-brn(beta[mt]))/vt;
}
}
for (nt = 1; nt < ll2+1; nt++) {
/* make matrix for Poisson's equation */
if (ndm == 0) {
at[nt-1] = gam1[nt];
bt[nt-1] = gam2[nt] - ECHARGE*vt/eps*(elec[nt]+hole[nt]);
ct[nt-1] = gam3[nt];
ページ(7)
simdd.c
ft[nt-1] = -(cgm[nt]-elec[nt]+hole[nt])*ECHARGE/eps
-gam1[nt]*volt[nt-1]-gam2[nt]*volt[nt] - gam3[nt]*volt[nt+1];
}
else if (ndm == 1) {
at[nt-1] = -dl1[nt-1]*gam1[nt];
bt[nt-1] = gam3[nt]*dl1[nt]-gam1[nt]*dl2[nt-1]-dudn[nt];
ct[nt-1] = dl2[nt]*gam3[nt];
ft[nt-1] = -(dl1[nt]*elec[nt]+dl2[nt]*elec[nt+1])*gam3[nt]
+ (dl1[nt-1]*elec[nt-1]+dl2[nt-1]*elec[nt])*gam1[nt]+rcmb[nt];
}
else if (ndm == 2) {
at[nt-1] = -dp1[nt-1]*gam1[nt];
bt[nt-1] = gam3[nt]*dp1[nt]-gam1[nt]*dp2[nt-1]+dudp[nt];
ct[nt-1] = dp2[nt]*gam3[nt];
ft[nt-1] = -(dp1[nt]*hole[nt]+dp2[nt]*hole[nt+1])*gam3[nt]
+(dp1[nt-1]*hole[nt-1]+dp2[nt-1]*hole[nt])*gam1[nt]-rcmb[nt];
}
}
}
/*****************************************************************************/
/*
Tridiagonal Matrix Solver
*/
/* a[i] : 三重対角行列の対角要素
*/
/* b[i] : 三重対角行列の非対角要素(上側)
*/
/* c[i] : 三重対角行列の非対角要素(下側)
*/
/* 注意点: 行列がN次元の場合,非対角要素はそれぞれN-1個になる。
*/
/* この場合、配列は */
/* b[0]-b[i]-b[N-2]に要素を格納する。
*/
/*****************************************************************************/
void tridag(double *a, double *b, double *c, double *f, double *u, int n)
{
double gam[NM], bet;
int j;
if (b[0] == 0.0) {
fprintf(stderr,"determinant in tridag() 1¥n");
exit(1);
}
bet = b[0];
u[0] = f[0]/bet;
for (j = 1; j < n; j++) {
gam[j] = c[j-1]/bet;
bet = b[j] - a[j]*gam[j];
if (bet == 0.0) {
fprintf(stderr,"determinant in tridag() 2 ...¥n");
exit(1);
}
u[j] = (f[j] - a[j]*u[j-1])/bet;
}
for (j = n-2; j >= 0; j--) {
u[j] = u[j] - gam[j+1]*u[j+1];
}
}
/*****************************************************************************/
/*
Bernoulli function
*/
/*****************************************************************************/
#define X1 -38.0
#define X2 -2.0e-6
#define X3 2.0e-6
#define X4 38.0
#define X5 40.0
double brn(double x)
{
if (x < X1)
return -x;
else if (X1< x && x <= X2)
return x/(exp(x)-1.0);
ページ(8)
simdd.c
simdd.c
else if (X2 < x && x <= X3)
return 1.0-x/2.0;
else if (X3 < x && x <= X4)
return x*exp(-x)/(1.0 - exp(-x));
else if (X4 < x && x <= X5)
return x*exp(-x);
else
return 0.0;
}
/*****************************************************************************/
/*
Recombination model (Shockley-Read-Hall model)
*/
/*****************************************************************************/
void recomb(int ndm, int ll, double *rcmb)
{
extern double cdin, taup, taun;
int mt;
double a1, b1;
if (ndm == 0) {
fprintf(stderr, "???: ndm == 0 in recomb()¥n");
exit(1);
}
if (ll > NM) {
fprintf(stderr, "???: ll > NM in recomb()¥n");
exit(1);
}
if (ll < 3) {
fprintf(stderr, "???: ll < 3 in recomb()¥n");
exit(1);
}
clear(rcmb,ll);
for (mt = 1; mt <= ll-1; mt++) {
a1 = elec[mt]*hole[mt] - cdin*cdin;
b1 = taup*(elec[mt] + cdin) + taun*(hole[mt]+cdin);
rcmb[mt] = a1/b1;
}
}
/*****************************************************************************/
/*
Mobility model (Simplified model)
*/
/*****************************************************************************/
void smpmob(int ndm, int ll, double *doud)
{
extern double temp;
int mt;
if (ndm == 0) {
printf("smpmob ndm = 0¥n");
exit(1);
}
if (ll > NM) {
printf("ll > NM in smpmob()¥n");
exit(1);
}
if (ll < 1) {
printf("ll < 1 in smpmob()¥n");
exit(1);
}
clear(doud,ll);
for (mt = 0; mt <= ll-1; mt++) {
if (ndm == 1)
doud[mt] = 4.0e5*pow(temp, -2.6);
else if (ndm == 2)
doud[mt] = 2.5e4*pow(temp, -2.3);
}
doud[ll] = doud[ll-1];
}
ページ(9)
ページ(10)