◆ Bsqad()

Sub Bsqad	(	T() As	Double,
		Bcoef() As	Double,
		N As	Long,
		K As	Long,
		X1 As	Double,
		X2 As	Double,
		Bquad As	Double,
		Info As	Long
	)

B-スプラインの積分値

目的: 本ルーチンはB-形式(T(), Bcoef(), N, K)のK次B-スプラインの, 区間[x1, x2]における積分を求める. 次数Kは, 2, 6または10点ガウス公式を区間[x1, x2]に(区間内に含まれるノットを使って)適用することにより, 最大20まで許される. 20を超える次数Kが必要な場合, f(x) = 1 としてBfqadを使用せよ.

引数

[in]	T()	配列 T(LT - 1) (LT >= N + K) ノットベクトル.
[in]	Bcoef()	配列 Bcoef(LBcoef - 1) (LBcoef >= N) B-スプライン係数.
[in]	N	B-スプライン係数の数. (N = ノット多重度の合計 - K)
[in]	K	B-スプラインの次数. (1 <= K <= 20)
[in]	X1	積分区間の下限値. (T(K) <= X1 <= T(N+1))
[in]	X2	積分区間の上限値. (T(K) <= X2 <= T(N+1))
[out]	Bquad	B-スプラインの区間[X1, X2]における積分値.
[out]	Info	= 0: 正常終了. = -1: パラメータ T() の誤り. = -2: パラメータ Bcoef() の誤り. = -3: パラメータ N の誤り. (N < K) = -4: パラメータ K の誤り. (K < 1 or K > 20) = -5: パラメータ X1 の誤り. (X1 < T(K) or X1 > T(N+1)) = -6: パラメータ X2 の誤り. (X2 < T(K) or X2 > T(N+1))

出典: SLATEC

使用例

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

  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_Bsqad()
    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 Info As Long, S As Double
    '-- 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
    '-- Compute integral 1/(1 + x^2) dx [0, 4] (= atan(4))
    Call Bsqad(T(), Bcoef(), N, K, A, B, S, Info)
    Debug.Print "S =", S, "S(true) =", Atn(4)
    Debug.Print "Info =", Info
End Sub

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

Info = 0