◆ Bsplvd()

Sub Bsplvd	(	T() As	Double,
		K As	Long,
		Nderiv As	Long,
		X As	Double,
		Ileft As	Long,
		Vnikx() As	Double,
		Info As	Long
	)

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

目的: 本ルーチンはXにおけるゼロでない基底関数の値および(Nderiv - 1)次以下の微分係数値を求める.

Ileftは T(Ileft-1) <= X <= T(Ileft) となるような入力値である. 次の呼び出しにより正しいIleftを生成することができる.
Call Interv(T(), N + 1, X, Ilo, Ileft, Info)

Interv
Sub Interv(Xt() As Double, Lxt As Long, X As Double, Ilo As Long, Ileft As Long, Info As Long)
Ileftを求める (BsplvdおよびBsplvnの入力用)

Bsplvdの出力は, 少なくてもK×Nderivの大きさの配列Vnikx(i,j)で, その列にXにおけるK個の非ゼロの基底関数およびそのNderiv-1右微分値が入る. ただし, i = 0〜K-1, j = 0〜Nderiv-1 である. これらの基底関数のインデックスは Ileft-K+i, i = 1〜K, K <= Ileft <= N である. i番目の基底関数の非ゼロ部分は区間 (T(i), T(i+K)), i = 0〜N-1) にある. X = T(Ileft) であればVnikxにはT(Ileft)における左極限値(左微分値)が入る. 特に, Ileft = N であれば右端点 X = T(N) において左極限値を求める. T(i) (i = K〜N) において左極限値を求めるためには, Xに次に小さなノットを設定しIntervを呼び出してIleftを求め, XにT(i)を設定しBsplvdを呼び出せばよい.

引数

[in]	T()	配列 T(LT - 1) (LT >= N + K) 長さ(N + K)のノットベクトル. ただし, NはB-スプライン基底関数の数 (= 多重度の合計 - K).
[in]	K	B-スプラインの次数. (K >= 1)
[in]	Nderiv	微分係数の数 + 1. (1 <= Nderiv <= K)
[in]	X	基底関数の引数. (T(K - 1) <= X <= T(N))
[in]	Ileft	T(Ileft - 1) <= X <= T(Ileft) となるような整数.
[out]	Vnikx()	配列 Vnikx(LVnikx1 - 1, LVnikx2 - 1) (LVnikx1 >= K, Lvnikx2 >= Nderiv) Xにおける非ゼロの基底関数およびその微分係数 (横方向に微分係数が並ぶ).
[out]	Info	= 0: 正常終了. = -1: パラメータ T() の誤り. = -2: パラメータ K の誤り. (K < 1) = -3: パラメータ Nderiv の誤り. (Nderiv < 1 または Nderiv > K) = -6: パラメータ Vnikx() の誤り. = 1: Ileftの値が正しくない. (X < T(Ileft-1) または X > T(Ileft))

出典: SLATEC

使用例

次の自然対数表をBint4により3次B-スプライン補間する. その係数を使ってBsplvdにより基底関数を計算し, ln(0.115) および微分値 ln'(0.115) の補間値を求める.

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

Sub Ex_Bsplvd()
    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 Vnikx() As Double, F As Double, Df As Double
    Dim Nderiv As Long, Ilo As Long, Ileft As Long, Info As Long, I 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
    '-- Compute basis function values by Bsplvd
    Xe = 0.115
    Ilo = 1
    Call Interv(T(), N + K, Xe, Ilo, Ileft, Info)
    If Ileft > N Then Ileft = N
    Nderiv = 2
    ReDim Vnikx(K - 1, Nderiv - 1)
    Call Bsplvd(T(), K, Nderiv, Xe, Ileft, Vnikx(), Info)
    If Info <> 0 Then
        Debug.Print "Error in Bsplvd: Info =", Info
        Exit Sub
    End If
    '-- Compute interpolated value
    F = 0: Df = 0
    For I = 0 To K - 1
        F = F + Bcoef(Ileft - K + I) * Vnikx(I, 0)
        Df = Df + Bcoef(Ileft - K + I) * Vnikx(I, 1)
    Next
    Debug.Print "ln(" + CStr(Xe) + ") =", F, "ln'(" + CStr(Xe) + ") =", Df
    Debug.Print "Info =", Info
End Sub

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

Info = 0