|
|
◆ Defin_r()
| Sub Defin_r |
( |
A As |
Double, |
|
|
B As |
Double, |
|
|
Result As |
Double, |
|
|
Info As |
Long, |
|
|
XX As |
Double, |
|
|
YY As |
Double, |
|
|
IRev As |
Long, |
|
|
Optional Neval As |
Long, |
|
|
Optional Eps As |
Double = -1, |
|
|
Optional L As |
Long = -1 |
|
) |
| |
有限区間の積分 (自動積分) (二重指数関数(DE)公式) (リバースコミュニケーション版)
- 目的
- 本ルーチンは要求精度を満たす[a, b]におけるf(x)の積分値を求める. ここで, f(x)は任意の関数である. ユーザーは変数IRevに従って必要なf(x)の関数値を計算しなければならない.
積分値は二重指数関数(DE)公式を使用した自動積分により求められる.
区間[-1, 1]における積分を変数変換により無限区間[-∞, ∞]の積分に変換する. ∫ f(x)dx [-1, 1] = ∫ g(t)dt [-∞, ∞]
ただし, x=φ(t), g(t) = f(φ(t))φ'(t)
φ(t)をうまく選べば、g(t)を指数関数の2乗で速く減衰する関数にすることができ, それに台形則を適用することにより効率よく積分を求めることができる. 本ルーチンでは次の変換関数を使用する. φ(t) = tanh((π/2)sinh(t))
- 引数
-
| [in] | A | 積分区間の下限 a. |
| [in] | B | 積分区間の上限 b. |
| [out] | Result | 求められた積分値. |
| [out] | Info | = 0: 正常終了.
= 1: 負側の減衰が遅いため収束判定条件を緩めた.
= 2: 正側の減衰が遅いため収束判定条件を緩めた.
= 3: 上の両方.
= 4: 積分計算が収束しなかった. (Epsが小さすぎる, 区間内に特異点があるなどの原因による) |
| [out] | XX | IRev = 1の場合, 関数値を求めるべき点を返す. |
| [in] | YY | IRev = 1の場合, 再呼び出し時に関数値を与えること. |
| [in,out] | IRev | リバースコミュニケーションの制御変数.
[in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない.
[out] 0 以外の時には下記処理を行ってから再び本ルーチンを呼び出すこと.
= 0: 処理終了. 正常終了かどうかはInfoをチェックせよ.
= 1: XXにおける関数値を求めYYに設定すること. YY以外の変数を変更してはならない. 通常(L = 0 の場合)は, 関数値としてf(XX)を設定する. f(x)が端点に特異性を持つ場合には, L = 1としf(XX)ではなくf2(XX)をYYに設定することにより通常より精度よくf(x)の積分値を計算することができる. ただし, f2(y)は次のように定義される.
f2(y) = f(a - y) (-(b - a)/2 < y < 0)
f2(y) = f(b - y) (0 < y <= (b - a)/2) |
| [out] | Neval | (省略可)
被積分関数の評価回数. |
| [in] | Eps | (省略可)
要求絶対誤差. (省略時 = 1.0e-14)
(Eps <= 0 であれば省略時の既定値とみなす) |
| [in] | L | (省略可)
Fの定義方法. (省略時 = 0)
= 0: 通常通り.
= 1: 端点での計算を避けるために変数変換して積分を行う. ユーザーはf(x)の代わりにf2(y)を計算する(IRevを参照せよ).
(その他の値であれば省略時の既定値とみなす) |
- 出典
- 森正武、「FORTRAN77数値計算プログラミング(増補版)」岩波書店 (1987)
- 使用例
- 次の定積分を求める.
∫ 1/(1 + x^2) dx [0, 4] (= atan(4))
Sub Ex_Defin_r()
Dim A As Double, B As Double, Result As Double, Info As Long
Dim XX As Double, YY As Double, IRev As Long
A = 0: B = 4
IRev = 0
Do
Call Qag_r(A, B, Result, Info, XX, YY, IRev) If IRev = 1 Then YY = 1 / (1 + XX ^ 2)
Loop While IRev <> 0
Debug.Print "S =", Result, "S(true) =", Atn(4)
Debug.Print "Info =", Info
End Sub
Sub Qag_r(A As Double, B As Double, Result As Double, Info As Long, XX As Double, YY As Double, IRev As Long, Optional AbsErr As Double, Optional Neval As Long, Optional EpsAbs As Double=-1, Optional EpsRel As Double=-1, Optional Key As Long=-1, Optional Limit As Long=-1, Optional Last As Long) 有限区間の積分 (適応自動積分) (15/21/31/41/51/61点ガウス・クロンロッド則) (リバースコミュニケーション版)
- 実行結果
S = 1.32581766366803 S(true) = 1.32581766366803
Info = 0
|
1.9.7