今回のコラムでは

　　f(x)=∫(0,x)[erfc(t)]^ndt

　　erfc(x)=2/(π)^(1/2)∫(x,∞)exp(-t^2)dt=1-erf(x)

を求める近似式について考えてみることにします．近似の範囲としては

　　ｘ：０〜５，ｎ：０．５〜１

の近似値が得られれば実用上十分です．

　もちろん最近のMathematicaなどのソフトを使えばこのような積分の数値処理

　　F[n_,x_]:=Nintegrate[Erfc[t]^n,{t,0,x}]

は無論，数式処理も相当できるのですが，ここでは近似式で代用できるかどうかについて考えてみます．

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

【１】問題の背景（熱拡散による半導体基板中の不純物濃度分布）

　モノリシック半導体集積回路では，ドナーやアクセプタなどの不純物元素をシリコン基板中に添加することによりｐ型やｎ型の領域を作っている．熱拡散はこの不純物ドーピングのための基本プロセスのひとつで，基板を不純物の雰囲気中で1000〜1300℃に加熱し，不純物をしみこませる方法である．

　基板表面から任意の距離と時間における不純物濃度の分布Ｎ(x,t)を不純物の供給が十分であり，表面濃度が常に一定値Ｎsに保たれているとして，熱拡散方程式

　　∂Ｎ／∂ｔ＝Ｄ∂^2Ｎ／∂ｘ^2

を解くと，不純物濃度の分布を表す式は

　　Ｎ(x,t)＝Ｎserfc(x/2(Dt)^(1/2))

と求められ，誤差補関数を用いて表されることがわかる．

−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

　また，電気伝導度Ｃと不純物濃度の関係は，実験的に

　　Ｃ＝αＮ^β　　　（Irvinの関係）

で表されることが知られている（α，βはＮの値によって異なるのだが，簡単のため一定として扱う）．

　所定温度，所定時間（Ｄ，ｔは一定）の熱処理後の拡散層の平均電気伝導度は

　　∫(0,x)Ｃ(u)du／∫(0,x)du

　＝αＮs^β／ｘ∫(0,x)erfc(u/2(Dt)^(1/2))du

となって，

　　f(x)=∫(0,x)[erfc(t)]^ndt

を求める問題に帰着されるというわけである．

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

【２】試行錯誤（その１）

　BASICには，組み込み関数として誤差関数

　　erf(x)=2/(π)^(1/2)∫(0,x)exp(-t^2)dt

や誤差補関数erfc(x)が含まれていないので，

　　f(x)=∫(0,x)[erfc(t)]^ndt

の数値積分は非常に扱いにくく，erfc(x)の値をコラム「超幾何関数を用いた確率分布の計算」に掲げたような級数展開で求めてから，引き続き台形則やシンプソン則で数値積分するという手順になります．すなわち，２重の級数の形ΣΣになってしまいますので，大層時間のかかるうえに計算誤差が大きいことが予想され大変始末が悪いのです．

　まず思いつくのは，ｘ＝０の周りのテイラー展開（マクローリン展開）を求める方法です．

　　f'(x)=[erfc(x)]^n

　　f''(x)=n[erfc(x)]^(n-1)[erfc(x)]'

　　f(3)(x)=n(n-1)[erfc(x)]^(n-2){[erfc(x)]'}^2+n[erfc(x)]^(n-1)[erfc(x)]''

　　f(4)(x)=n(n-1)(n-2)[erfc(x)]^(n-3){[erfc(x)]'}^3+3n(n-1)[erfc(x)]^(n-2)[erfc(x)]'[erfc(x)]''+n[erfc(x)]^(n-1)[erfc(x)]'''

　　[erfc(x)]'=-2/π^(1/2)exp(-x^2)

　　[erfc(x)]''=4/π^(1/2)exp(-x^2)x

　　[erfc(x)]'''=4/π^(1/2)exp(-x^2)(1-2x^2)

　　[erfc(x)](4)=-8/π^(1/2)exp(-x^2)(-3x+4x^3)

より

　　f(x)=f(0)+f'(0)x+f''(0)x^2/2+f(3)(0)x^3/6+f(4)(0)x^4/24+･･･

=x-2n/π^(1/2)x^2/2+4n(n-1)/πx^3/6-8n(n-1)(n-2)/π^(3/2)x^4/24+･･･

となりますが，積分上限ｘが大きくなるにつれ収束速度が問題になりそうです．そのような場合，漸近展開式

　　erfc(x)=exp(-x^2)/xπ^(1/2)[1+Σ(-1)^k(2k-1)!!/(2x^2)^k]

もありますが・・・．

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

【３】試行錯誤（その２）

　誤差関数の近似式として，中心確率で定義した近似式

　　Φc(x)=1/(2π)^(1/2)∫(-x,x)exp(-t^2/2)dt

　　　　　≒{1-exp(-2x^2/π)}^(1/2)

が，簡単なわりには誤差関数のかなりよい近似になるので知っておくと便利です．

　これを用いると

　　erfc(x)=2{1-Φ(2^(1/2)x)}

　　　　　≒1-{1-exp(-4x^2/π)}^(1/2)

となるのですが，

　　f(x)=∫(0,x)[erfc(t)]^ndt

を計算しやすい形に変形することはできませんでした．

　次に，正規分布の代用としてロジスティック分布を用いてみるというアイディアが浮かびました．ロジスティック分布は正規分布によく似た釣り鐘型曲線で，しかも積分値を求めやすい関数となっています．

f(x)=exp(-(x-α)/β)/{1+exp(-(x-α)/β)}^2

=1/4βsech^2{(x-α)/2β}

F(x)=1/{1+exp(-(x-α)/β)}

=1/2[1+tanh{(x-α)/2β}]

-∞<x<∞

mean=α (mode,medianとも)

variance=β^2π^2/3

　累積分布関数F(x)がロジスティック曲線になる分布がロジスティック分布です．α=0,β=1の場合が標準ロジスティック分布

f(x)=exp(-x)/{1+exp(-x)}^2

でその平均値０，分散はπ^2／３となります．

　ロジスティック分布も裾の重い分布の１つで（−２，２）の外の確率が０．２３８４，（−３，３）の外の確率が０．０９４９です．この分布の裾の重さはほぼ自由度８〜９のｔ分布に相当します．また，尖度は正規分布が３に対してロジスティック分布が４．２です．

　正規分布もロジスティック分布もベル型曲線を描きますが，正規分布

　　f(x)=1/√(2π)exp(-x^2/2)

の粗い近似式として，平均値０，分散１に規格化したロジスティック分布

　　f(x)=exp(-πx/√3)/{1+exp(-πx/√3)}^2

　　F(x)=1/{1+exp(-πx/√3)}

はしばしば用いられます．

　しかし，ロジスティック近似

　　erfc(x)=2{1-Φ(2^(1/2)x)}

　　　　　≒2{1-F(2^(1/2)x)}=2{1-1/{1+exp(-√2πx/√3)}}

をもってしても計算しやすい形に変形したり，収束を速くしたりすることには役立ちませんでした．

　この問題は２年前からときどき思い出しては考えてみるのですが，いまだ解決には至っておりません．何か「よい近似式」はないものでしょうか？

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝