◆ Bsplvn()

Sub Bsplvn	(	T() As	Double,
		Jhigh As	Long,
		K As	Long,
		Index As	Long,
		X As	Double,
		Ileft As	Long,
		Vnikx() As	Double,
		Info As	Long
	)

B-スプライン基底関数値

目的: 本ルーチンはXにおけるmax(Jhigh, (j+1)*(Index-1))次のゼロでない基底関数の値を求める. ただし, T(K-1) <= X <= T(N) で, j = iworkはIndex = 1の最初の呼び出し時にルーチン内で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の入力用)

BsplvnはBsplvdと同じアルゴリズムを使って計算を行う. 基底関数だけを求める場合には Jhigh = K および Index = 1 と設定するとBsplvdを呼び出すよりも速く計算できる. 微分係数を求める場合にはBsplvdが適している. 左極限値を求めるやり方はBsplvdを参照せよ.

引数

[in]	T()	配列 T(LT - 1) (LT >= N + K) 長さ(N + K)のノットベクトル. ただし, NはB-スプライン基底関数の数 (= 多重度の合計 - K).
[in]	Jhigh	求めるB-スプラインの次数. (1 <= Jhigh <= K)
[in]	K	可能な最も高い次数. (K >= 1)
[in]	Index	= 1: 最初から計算を行う. 1, 2, ..., Jhigh次のすべてのB-スプラインの値が順に求められる. = 2: j+1, j+2, ..., Jhigh次のB-スプラインの値だけを計算する. 前回呼び出しにより求められたj次までの基底関数は保持されていなければならない.
[in]	X	基底関数の引数. (T(K - 1) <= X <= T(N))
[in]	Ileft	T(Ileft - 1) <= X <= T(Ileft) となるような整数.
[out]	Vnikx()	配列 Vnikx(LVnikx - 1) (LVnikx >= K) 基底関数の値.
[out]	Info	= 0: 正常終了. = -1: パラメータ T() の誤り. = -2: パラメータ Jhigh の誤り. (Jhigh < 1 または Jhigh > K)) = -3: パラメータ K の誤り. (K < 1) = -4: パラメータ Index の誤り. (Index <> 1 かつ Index <> 2) = -7: パラメータ Vnikx() の誤り. = 1: Ileftの値が正しくない. (X < T(Ileft-1) または X > T(Ileft))

出典: SLATEC

使用例

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

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

Sub Ex_Bsplvn()
    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
    Dim 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 Bsplvn
    Xe = 0.115
    Ilo = 1
    Call Interv(T(), N + K, Xe, Ilo, Ileft, Info)
    If Ileft > N Then Ileft = N
    ReDim Vnikx(K - 1)
    Call Bsplvn(T(), K, K, 1, Xe, Ileft, Vnikx(), Info)
    If Info <> 0 Then
        Debug.Print "Error in Bsplvn: Info =", Info
        Exit Sub
    End If
    '-- Compute interpolated value
    F = 0
    For I = 0 To K - 1
        F = F + Bcoef(Ileft - K + I) * Vnikx(I)
    Next
    Debug.Print "ln(" + CStr(Xe) + ") =", F
    Debug.Print "Info =", Info
End Sub

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

Info = 0