◆ Bfqad_r()

Sub Bfqad_r	(	T() As	Double,
		Bcoef() As	Double,
		N 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
	)

Integral of product of arbitrary function and B-representation of B-spline (reverse communication version)

Purpose: This routine computes the integral on [X1, X2] of a product of a function f(x) and the Id-th derivative of a K-th order B-spline, using the B-representation (T(), Bcoef(), N, K). [X1, X2] must be a subinterval of [T(K-1), T(K)].

Then integration routine using adaptive 8-point Gauss-Legendre rule integrates the product on subintervals of [X1, X2] formed by included (distinct) knots.

Parameters

[in]	T()	Array T(LT - 1) (LT >= N + K) Knot vector.
[in]	Bcoef()	Array Bcoef(LBcoef - 1) (LBcoef >= N) B-spline coefficients.
[in]	N	Number of B-spline coefficients. (N = sum of knot multiplicities - K)
[in]	K	Order of the B-spline. (1 <= K)
[in]	Id	Order of the spline derivative. (0 <= Id <= K - 1) Id = 0 gives the spline function.
[in]	X1	Lower end point of quadrature interval. (T(K) <= X1 <= T(N+1))
[in]	X2	Upper end point of quadrature interval. (T(K) <= X2 <= T(N+1))
[in]	Tol	Desired accuracy for the quadrature. (Eps < Tol <= 0.1, where Eps is the double precision unit roundoff (= D1mach(4)))
[out]	Quad	Integral of f(x)*(Id-th derivative of a K-th order B-spline) on [X1, X2].
[out]	Info	= 0: Successful exit. = -1: The argument T() is invalid. = -2: The argument Bcoef() is invalid. = -3: The argument N had an illegal value. (N < K) = -4: The argument K had an illegal value. (K < 1) = -5: The argument Id had an illegal value. (Id < 0 or Id >= K) = -6: The argument X1 had an illegal value. (X1 < T(K) or X1 > T(N+1)) = -7: The argument X2 had an illegal value. (X2 < T(K) or X2 > T(N+1)) = -8: The argument Tol had an illegal value. (Tol < Eps or Tol > 0.1) = 1: Some quadrature on (X1, X2) does not meet the requested tolerance.
[out]	XX	When returned with IRev = 1, XX contains the abscissa where the function value should be evaluated and given in the next call.
[in]	YY	When returned with IRev = 1, the function value f(XX) should be given in YY in the next call.
[in,out]	IRev	Control variable for reverse communication. [in] Before first call, IRev should be initialized to zero. On succeeding calls, IRev should not be altered. [out] If IRev is not zero, set the required values to the specified variable as follows and call the routine again. = 0: Computation finished. See return code in Info. = 1: User should set the function values at XX in YY. Do not alter any variables other than YY.

Reference: SLATEC

Example Program

Using the following table, compute S = integral of 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_Bfqad_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 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
    '-- Compute integral 1/(1 + x^2) dx [0, 4] (= atan(4))
    Id = 0: Tol = 0.0000000001 '1.0e-10
    IRev = 0
    Do
        Call Bfqad_r(T(), Bcoef(), N, 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

Example Results: S = 1.32679961538462 S(true) = 1.32581766366803

Info = 0