◆ Bsplev()

Sub Bsplev	(	T() As	Double,
		Ad() As	Double,
		N As	Long,
		K As	Long,
		Nderiv As	Long,
		X As	Double,
		Inev As	Long,
		Svalue() As	Double,
		Info As	Long
	)

B-スプライン関数値および微分値

目的: 本ルーチンはB-形式(T(), A(), N, K)のB-スプライン関数のXにおける関数値および微分値を求め, Svalue(i) (i = 0〜Nderiv-1)に格納する. ただし, T(K-1) <= X <= T(N) である.

Ad(i)は, Nderiv = 1の場合, B-スプライン係数A(i) (i = 0〜N-1)であってよい. その他の場合, 次の呼び出しにより前もって計算しておかなけらばならない.
Call Bspldr(T(), A(), N, K, Nderiv, Ad())

Bspldr
Sub Bspldr(T() As Double, A() As Double, N As Long, K As Long, Nderiv As Long, Ad() As Double, Info As Long)
B-スプライン係数より差分商を構成(Bsplevによる微分値の計算用)

X = T(i) (K-1 <= i <= N-1)であれば, 右極限値が求められる. T(i)において左極限値(左微分値)を求めるためには, Nをi-1で置き替えて X = T(i) (K <= i <= N) とせよ.

引数

[in]	T()	配列 T(LT - 1) (LT >= N + K) ノットベクトル.
[in]	Ad()	配列 Ad(LAd - 1) (LAd >= (2N - Nderiv + 1)Nderiv/2) 差分商の表.
[in]	N	B-スプライン係数の数. (N = ノット多重度の合計 - K)
[in]	K	B-スプラインの次数. (K >= 1)
[in]	Nderiv	微分係数の数. (1 <= Nderiv <= K) Nderiv = 1 であれば 0次微分係数 = 関数値を返す.
[in]	X	引数. (T(K-1) <= X <= T(N))
[in,out]	Inev	初期化パラメータ. [in] 最初の呼び出し時, 1に設定しておかなければならない. [out] 最初の呼び出し後の処理のための情報. ユーザーが変更してはならない. 異なるスプラインは異なったInevパラメータを持つ.
[out]	Svalue()	配列 Svalue(LSvalue - 1) (LSvalue >= Nderiv) Svalue(0)にspline関数値, 残りにNderiv-1個の微分係数を返す.
[out]	Info	= 0: 正常終了. = -1: パラメータ T() の誤り. = -2: パラメータ Ad() の誤り. = -3: パラメータ N の誤り. (N < K) = -4: パラメータ K の誤り. (K < 1) = -5: パラメータ Nderiv の誤り. (Nderiv < 1 or Nderiv > K) = -6: パラメータ X の誤り. (X < T(K) or X > T(N+1)) = -8: パラメータ Svalue() の誤り. = 1: T(K-1)において左極限値は求められない.

出典: SLATEC

使用例

次の自然対数表を3次B-スプライン補間(自然スプライン)して, ln(0.115) およびその微分値 in'(0.115) を求める.

  x       ln(x)
------ ---------
 0.10   -2.3026
 0.11   -2.2073
 0.12   -2.1203
 0.13   -2.0402
------ ---------

Sub Ex_Bsplev()
    Const Ndata = 4
    Dim X(Ndata - 1) As Double, Y(Ndata - 1) As Double, Xe 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 Nderiv As Long, Ad() As Double
    Dim Inev As Long, Svalue() As Double, Info As Long
    '-- Data
    X(0) = 0.1: Y(0) = -2.3026
    X(1) = 0.11: Y(1) = -2.2073
    X(2) = 0.12: Y(2) = -2.1203
    X(3) = 0.13: Y(3) = -2.0402
    '-- B-representation of cubic 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
    '-- Prepare divided difference table for Bsplev
    Nderiv = K
    ReDim Ad((2 * N - Nderiv + 1) * Nderiv / 2 - 1)
    Call Bspldr(T(), Bcoef(), N, K, Nderiv, Ad(), Info)
    If Info <> 0 Then
        Debug.Print "Error in Bspldr: Info =", Info
        Exit Sub
    End If
    '-- Evaluate by Bsplev
    Xe = 0.115
    Inev = 1
    ReDim Svalue(Nderiv - 1)
    Call Bsplev(T(), Ad(), N, K, Nderiv, Xe, Inev, Svalue(), Info)
    If Info <> 0 Then
        Debug.Print "Error in Bsplev: Info =", Info
        Exit Sub
    End If
    Debug.Print "ln(" + CStr(Xe) + ") =", Svalue(0), "ln'(" + CStr(Xe) + ") =", Svalue(1)
    Debug.Print "Info =", Info
End Sub

実行結果: ln(0.115) = -2.16266 ln'(0.115) = 8.68833333333335

Info = 0