XLPack 7.0
XLPack 数値計算ライブラリ (Excel VBA) リファレンスマニュアル
読み取り中…
検索中…
一致する文字列を見つけられません

◆ Defin()

Sub Defin ( F As  LongPtr,
A As  Double,
B As  Double,
Result As  Double,
Info As  Long,
Optional Neval As  Long,
Optional Eps As  Double = -1,
Optional L As  Long = -1 
)

有限区間の積分 (自動積分) (二重指数関数(DE)公式)

目的
本ルーチンは要求精度を満たす[a, b]におけるf(x)の積分値を求める. ここで, f(x)はユーザー定義ルーチンFにより与えられる関数である.
積分値は二重指数関数(DE)公式を使用した自動積分により求められる.

区間[-1, 1]における積分を変数変換により無限区間[-∞, ∞]の積分に変換する.
∫ f(x)dx [-1, 1] = ∫ g(t)dt [-∞, ∞]
ただし, x=φ(t), g(t) = f(φ(t))φ'(t)
ここで, φ(t)は単調で, かつ φ(-∞) = -1, φ(∞) = 1 とする.
φ(t)をうまく選べば、g(t)を指数関数の2乗で速く減衰する関数にすることができ, それに台形則を適用することにより効率よく積分を求めることができる. 本ルーチンでは次の変換関数を使用する.
φ(t) = tanh((π/2)sinh(t))
DE公式は端点に特異性があるときにも適用できる.
引数
[in]F被積分関数f(x)を求めるユーザー定義ルーチンで, 通常(L = 0 の場合)は次のように定義すること.
Function F(X As Double) As Double
F = f(X)の計算値
End Function
f(x)が端点に特異性を持つ場合には, f2(y)が下のように表されるものとして, Fではf(x)ではなくf2(y)を定義しパラメータ L = 1 としてDefinを呼び出すことにより, 通常より精度よく[a, b]におけるf(x)の積分値を計算することができる.
  f2(y) = f(a - y) (-(b - a)/2 < y < 0)
  f2(y) = f(b - y) (0 < y <= (b - a)/2)
Fの中でXを変更してはならない.
[in]A積分区間の下限 a.
[in]B積分区間の上限 b.
[out]Result求められた積分値.
[out]Info= 0: 正常終了.
= 1: 負側の減衰が遅いため収束判定条件を緩めた.
= 2: 正側の減衰が遅いため収束判定条件を緩めた.
= 3: 上の両方.
= 4: 積分計算が収束しなかった. (Epsが小さすぎる, 区間内に特異点があるなどの原因による)
[out]Neval(省略可)
被積分関数の評価回数.
[in]Eps(省略可)
要求絶対誤差. (省略時 = 1.0e-14)
(Eps <= 0 であれば省略時の既定値とみなす)
[in]L(省略可)
Fの定義方法. (省略時 = 0)
= 0: 通常通り.
= 1: 端点での計算を避けるために変数変換して積分を行う. ユーザー定義ルーチンFではf(x)ではなくf2(y)を定義する(Fを参照せよ).
(その他の値であれば省略時の既定値とみなす)
出典
森正武、「FORTRAN77数値計算プログラミング(増補版)」岩波書店 (1987)
使用例
次の定積分を求める.
∫ 1/(1 + x^2) dx [0, 4] (= atan(4))
Function F1(X As Double) As Double
F1 = 1 / (1 + X ^ 2)
End Function
Sub Ex_Defin()
Dim A As Double, B As Double, Result As Double, Info As Long
A = 0: B = 4
Call Defin(AddressOf F1, A, B, Result, Info)
Debug.Print "S =", Result, "S(true) =", Atn(4)
Debug.Print "Info =", Info
End Sub
Sub Defin(F As LongPtr, A As Double, B As Double, Result As Double, Info As Long, Optional Neval As Long, Optional Eps As Double=-1, Optional L As Long=-1)
有限区間の積分 (自動積分) (二重指数関数(DE)公式)
実行結果
S = 1.32581766366803 S(true) = 1.32581766366803
Info = 0