◆ Pfqad_r()

Sub Pfqad_r	(	C() As	Double,
		Xi() As	Double,
		Lxi As	Long,
		K As	Long,
		Id As	Long,
		X1 As	Double,
		X2 As	Double,
		Tol As	Double,
		Quad As	Double,
		Info As	Long,
		XX As	Double,
		YY As	Double,
		IRev As	Long
	)

任意関数×B-スプラインの積分値 (PP(区分多項式)形式) (リバースコミュニケーション版)

目的: 本ルーチンは任意の関数f(x)とPP(区分多項式)形式(C(), Xi(), Lxi, K)のK次B-スプラインのId次微分値の積の区間[X1, X2]における積分を求める. [X1, X2]は通常は[Xi(0), Xi(Lxi)]の部分区間である.

8点ガウス公式を使用した適応求積ルーチンにより, 区間内に含まれる区切り点により構成される部分区間[X1, X2]の積分を求める. fが定義されていれば区間[Xi(0), Xi(Lxi)]外の積分も可能である.

引数

[in]	C()	配列 C(LC1 - 1, LC2 - 1) (LC1 >= K, LC2 >= Lxi) 区分多項式の区切り点における右微分係数.
[in]	Xi()	配列 Xi(LXi - 1) (LXi >= Lxi + 1) 区分多項式の区切り点.
[in]	Lxi	区分多項式の小区間数.
[in]	K	B-スプラインの次数. (K >= 1)
[in]	Id	スプラインの微分係数の次数. (0 <= Id <= K - 1) Id = 0の場合, 関数値となる.
[in]	X1	積分区間の下限値. (通常は Xi(0) <= X1 <= Xi(Lxi))
[in]	X2	積分区間の上限値. (通常は Xi(0) <= X2 <= Xi(Lxi))
[in]	Tol	積分の要求精度. (Eps < Tol <= 0.1) (ただし, Epsは単位丸め誤差(= D1mach(4))である)
[out]	Quad	f(X)*(K次B-スプラインのId次微分係数) の区間[X1, X2]における積分値.
[out]	Info	= 0: 正常終了. = -1: パラメータ C() の誤り. = -2: パラメータ Xi() の誤り. = -3: パラメータ Lxi の誤り. (Lxi < 1) = -4: パラメータ K の誤り. (K < 1) = -5: パラメータ Id の誤り. (Id < 0 or Id >= K) = -8: パラメータ Tol の誤り. (Tol < DTol or Tol > 0.1) = 1: 区間(X1, X2)におけるいくつかの積分値が要求精度を満たさなかった.
[out]	XX	IRev = 1の場合, 関数値を求めるべき点を返す.
[in]	YY	IRev = 1の場合, 再呼び出し時に関数値f(XX)を与えること.
[in,out]	IRev	リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の時には下記処理を行ってから再び本ルーチンを呼び出すこと. = 0: 処理終了. 正常終了かどうかはInfoをチェックすること. = 1: XXにおける関数値f(XX)を求め YYに設定すること. YY以外の変数を変更してはならない.

出典: SLATEC

使用例

次の数表を用いて S = ∫ 1/(1 + x^2) dx [0, 4] (= atan(4)) を求める. (f(X) = 1)

  x    1/(1 + x^2)
----- -------------
 -1        0.5
  0        1 
  1        0.5 
  2        0.2 
  3        0.1 
  4        0.05882
  5        0.03846
----- -------------

Sub Ex_Pfqad_r()
    Const Ndata = 7, A = 0, B = 4
    Dim X(Ndata - 1) As Double, Y(Ndata - 1) As Double, D(Ndata - 1) As Double
    Dim Ibcl As Long, Ibcr As Long, Fbcl As Double, Fbcr As Double, Kntopt As Long
    Dim T(Ndata + 5) As Double, Bcoef(Ndata + 1) As Double, N As Long, K As Long
    Dim Lxi As Long, C(3, Ndata - 2) As Double, Xi(Ndata - 1) As Double
    Dim Id As Long, Tol As Double, Info As Long, S As Double
    Dim XX As Double, YY As Double, IRev As Long
    '-- Data
    X(0) = -1: Y(0) = 0.5
    X(1) = 0: Y(1) = 1
    X(2) = 1: Y(2) = 0.5
    X(3) = 2: Y(3) = 0.2
    X(4) = 3: Y(4) = 0.1
    X(5) = 4: Y(5) = 0.05882
    X(6) = 5: Y(6) = 0.03846
    '-- B-spline interpolation
    Ibcl = 2: Fbcl = 0: Ibcr = 2: Fbcr = 0 '-- Natural spline
    Kntopt = 1
    Call Bint4(X(), Y(), Ndata, Ibcl, Ibcr, Fbcl, Fbcr, Kntopt, T(), Bcoef(), N, K, Info)
    If Info <> 0 Then
        Debug.Print "Error in Bint4: Info =", Info
        Exit Sub
    End If
    '-- Convert to PP form
    Call Bsplpp(T(), Bcoef(), N, K, C(), Xi(), Lxi, Info)
    If Info <> 0 Then
        Debug.Print "Error in Bsplpp: Info =", Info
        Exit Sub
    End If
    '-- Compute integral 1/(1 + x^2) dx [0, 4] (= atan(4))
    Id = 0: Tol = 0.0000000001 '1.0e-10
    IRev = 0
    Do
        Call Pfqad_r(C(), Xi(), Lxi, K, Id, A, B, Tol, S, Info, XX, YY, IRev)
        If IRev = 1 Then YY = 1
    Loop While IRev <> 0
    Debug.Print "S =", S, "S(true) =", Atn(4)
    Debug.Print "Info =", Info
End Sub

実行結果: S = 1.32679961538462 S(true) = 1.32581766366803

Info = 0