[VBA] tmexp(-atn) の有限/無限積分

tmexp(-atn) の積分関数

 次のような積分関数 $I(m,n,a;x)$ を定義します。
 
\[I(m,n,a;x)=\int_{0}^{x}t^m\exp (-at^n)dt\quad (a\geq 0)\]
 $a$ は 0 または正の実数、$m,\:n$ はそれぞれ正の整数です。
 非常に広い範囲の積分を含む関数で、たとえば $I(0,2,1;x)$ は誤差関数
 
\[\mathrm{erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}\exp (-t^2)dt\]
において係数 $2/\sqrt{\pi}$ を省いたものとなります。
 

tmexp(-atn) の積分を計算するユーザー定義関数

 $I(m,n,a;x)$ を計算するマクロです。

 'x^m exp(-ax^n)の積分

 Function IEXPXN(x As Double, et1 As Integer, _
 et2 As Integer, a As Double) As Variant

 Dim p As Double, q As Double, r As Double
 Dim h As Double, s As Double
 Dim m As Integer, k As Integer

 'et1, et2, a が負であればエラーを返します
 If et1 < 0 Or et2 < 0 Or a < 0 Then

 IEXPXN = CVErr(xlErrNA)

 Else

 '分割数
 m = 1024

 '区間の幅
 h = x / (2 * m)

 s = 0

 For k = 0 To m - 1

 p = 2 * k * h

 q = (2 * k + 1) * h

 r = (2 * k + 2) * h

 s = s + exp(-a * p ^ et2) * p ^ et1 + _
 4 * exp(-a * q ^ et2) * q ^ et1 + exp(-a * r ^ et2) * r ^ et1

 Next k

 IEXPXN = s * h / 3

 End If

 End Function

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

=IEXPXN(x,m,n,a)

と入力します。第一引数に x の値を指定することに注意してください。
 

I(1,2,1;x) と I(2,2,1;x) のグラフ

 例として
 
\[\begin{align*}I(1,2,1;x)&=\int_{0}^{x}t\exp (-t^2)dt\\[6pt]
I(2,2,1;x)&=\int_{0}^{x}t^2\exp (-t^2)dt\end{align*}\]
を計算してシートにグラフを描いてみると次のようになります。

 指数関数の積分関数
 

xkexp(-ax), xnexp(-ax2)の無限積分

 $x^k\exp (-ax)$ および $x^k\exp (-ax^2)$ の積分は次の式で与えられます。
 (証明はこちらの記事を参照してください

\[I(a,k)=\int_{0}^{\infty} e^{-ax}x^kdx=\frac{k!}{a^{k+1}}\]
\[J(a,k)=\begin{cases}
\displaystyle \int_{0}^{\infty} e^{-ax^2}x^{2n}dx=\frac{(2n-1)!!}{2^{n+1}}\sqrt{\frac{\pi}{a^{2n+1}}} & (k=2n)\\[6pt]
\displaystyle \int_{0}^{\infty} e^{-ax^2}x^{2n+1}dx=\frac{n!}{2a^{n+1}} & (k=2n+1)\end{cases}\]

 ただし $a$ は正の実数、$k,\:n$ は自然数です。この積分は IEXPXN関数で上限を十分大きくとって計算させることもできますが、IEX 関数を使ったほうが精度が高くなります。
 

xkexp(-ax), xnexp(-ax2)の無限積分を計算するマクロ

 $x^k\exp (-ax)$ と $x^k\exp (-ax^2)$ の無限積分を計算するマクロです。

 '(x^k)exp(-ax), (x^k)exp(-ax^2) の積分

 Function IEX(a As Double, n As Integer, _
 Optional t As Boolean = False) As Variant

 Dim k As Integer
 Dim d As Double, ct As Double

 'n が負または a が 0 以下であればエラーを返して計算終了
 If n < 0 Or a <= 0 Then

 IEX = CVErr(xlErrNA)

 Exit Function

 End If

 'オプション引数が True のときは (x^k)exp(-ax) を積分します
 Select Case t

 Case True

 d = 1 / a

 For k = 1 To n

  d = d * k / a

 Next k

 IEX = d

 'オプション引数が False のときは (x^k)exp(-ax^2) を積分します
 Case Else

 d = 1 / (2 * a)
 ct = Sqr(4 * Atn(1) / a) / 2

 If n = 0 Then

 IEX = ct

 ElseIf n Mod 2 = 0 Then

 For k = 1 To n - 1

  d = d * (2 * k + 1) / (2 * a)

 Next k

 IEX = d * ct

 Else

 For k = 1 To n

  d = d * k / a

 Next k

 IEX = d

 End If

 End Select

 End Function

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

=IEX(a,k[,被積分関数の種類])

と入力します。a は正の実数、k は 0 以上の整数を指定してください。3つめの引数はオプションで False または 0 のときは
 
\[J(a,k)=\begin{cases}
\displaystyle \int_{0}^{\infty} e^{-ax^2}x^{2n}dx=\frac{(2n-1)!!}{2^{n+1}}\sqrt{\frac{\pi}{a^{2n+1}}} & (k=2n)\\[6pt]
\displaystyle \int_{0}^{\infty} e^{-ax^2}x^{2n+1}dx=\frac{n!}{2a^{n+1}} & (k=2n+1)\end{cases}\]
を計算します。True もしくは 1 を指定すると
 
\[I(a,k)=\int_{0}^{\infty} e^{-ax}x^kdx=\frac{k!}{a^{k+1}}\]
を計算します。横軸に $a$ をとり、$k=1,\:2,\:3$ として $J(a,k)$ を計算させると次のようなグラフとなります。

 Excel指数関数x累乗の無限積分

 また、a に 1, k に 0 を指定するとガウス積分
 
\[\int_{0}^{\infty} e^{-ax^2}dx=\frac{\sqrt{\pi}}{2}=0.886\]
を得ることができます。 ≫ VBA 数値計算

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

コメントをどうぞ

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

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