[VBA] 指数積分関数 Ei(x) の数値計算とグラフ

 この記事のマクロを使用する場合は 数値計算用コード をコピーしてモジュールに貼りつけてください(定期的に更新されるので常に最新版をコピーして上書きするようにしてください)。この記事に掲載されるマクロも全てそこに入っています。

指数積分関数

 指数積分関数 とは

\[\mathrm{Ei}\:(x)=-\Gamma (0,-x)=-\int_{-x}^{\infty}e^{-t}\:t^{a-1}dx\]

と定義される関数です。これは
 
\[-\mathrm{Ei}\:(-x)=\Gamma (0,x)=\int_{z}^{\infty}e^{-t}\:t^{a-1}dx\]
と書き直すこともできます。ここに $\Gamma (a,x)$ は第 2 種不完全ガンマ関数です。つまり右辺の積分値を求めるために不完全ガンマ関数を利用できるということです。ただし精度を確保するために $x$ が大きいところでは以下の漸近展開を用いることにします。
 
\[\mathrm{Ei}\:(x)\simeq\frac{e^x}{x}\left(1+\frac{1!}{x}+\frac{2!}{x^2}+\frac{3!}{x^3}\:+\cdots\:\right)\]
 

指数積分を計算する EI関数

 指数積分関数を計算するマクロです。

 '指数積分関数 (C)BlogCat

 Function EI(x As Double) As Double

 Dim d As Double, k As Double, s As Double
 Dim z1 As Complex, mz As Complex
 Dim gmz As Complex, mgmz As Complex

 'x が小さいところでは不完全ガンマ関数を用います

 If x <= 50 Then

 'x を複素数に変換
 z1.re = x
 z1.im = 0

 '-z
 mz = CPXCLC(CPX(-1, 0), z1, 3)

 'Γ(0,-z)
 gmz = INCGMA2(0, mz)

 '-Γ(0,-z)
 mgmz = CPXCLC(CPX(-1, 0), gmz, 3)

 EI = mgmz.re

 'x が大きいところでは漸近展開を用います

 Else

 d = 1
 s = 1

 For k = 1 To 10

 d = d * k / x

 s = s + d

 Next k

 EI = s * exp(x) / x

 End If

 End Function

 EI関数をワークシートで使用するときは

=EI(x)

のように入力してください。x は数値です。たとえば

=EI(10)

と入力すると

2492.22897624188

という値が返りますが、これは 12 桁以上の精度があることを確認済みです。引数 x の上限は 120 あたりまで戻り値の精度 8 桁以上が保たれます(おそらく 160 あたりまでは大丈夫です)。

指数積分関数 EI(x) のグラフ

 EI関数を使って実変数の指数積分関数
 
\[\mathrm{Ei}\:(x)=-\Gamma (0,-x)=-\int_{-x}^{\infty}e^{-t}\:t^{a-1}dx\]
をワークシートで計算してみます。たとえば、セル A1 に x の値が入力されているときには

=EI(A2)

と入力して計算させます。データを作成してグラフを描いてみると次のようになります。

 Excel VBA 指数積分関数Ei(x)のグラフ

 実用では次の積分を求めたいケースが多いと思います。
 
\[-\mathrm{Ei}\:(-x)=\Gamma (0,x)=\int_{x}^{\infty}e^{-t}\:t^{a-1}dx\]
 この場合は

=-EI(-A2)

によって計算することができます。グラフは原点について反転させた形になります。

 Excel VBA 指数積分のグラフ

 原点付近では $1/t$ が発散するので値をもちません。
 また $t$ の増加とともに $e^{-t}/t$ は急減する関数ですから、$x$ の値が十分に大きいところを積分の下限とすると、積分値はほとんど 0 となります。 ≫ VBA 数値計算

スポンサーリンク
末尾広告
末尾広告

コメントをどうぞ

メールアドレスが公開されることはありません。

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください