◆ Ppqad()

Sub Ppqad	(	C() As	Double,
		Xi() As	Double,
		Lxi As	Long,
		K As	Long,
		X1 As	Double,
		X2 As	Double,
		Pquad As	Double,
		Info As	Long
	)

Integral of PP (piecewise polynomial) form of B-spline

Purpose: This routine computes the integral on [X1, X2] of a K-th order B-spline using the piecewise polynomial representation (C(), Xi(), Lxi, K). Here the Taylor expansion about the left end point Xi(j-1) of the j-th interval is integrated and evaluated on subintervals of [X1, X2] which are formed by included break points. Integration outside [Xi(0), Xi(Lxi)] is permitted.

Parameters

[in]	C()	Array C(LC1 - 1, LC2 - 1) (LC1 >= K, LC2 >= Lxi) Right derivatives at break points.
[in]	Xi()	Array Xi(LXi - 1) (LXi >= Lxi + 1) Break points.
[in]	Lxi	Number of polynomial pieces.
[in]	K	Order of the B-spline. (K >= 1)
[in]	X1	Lower end point of quadrature interval. (Normally Xi(0) <= X1 <= Xi(Lxi))
[in]	X2	Upper end point of quadrature interval. (Normally Xi(0) <= X2 <= Xi(Lxi))
[out]	Pquad	Integral of the B-spline over [X1, X2].
[out]	Info	= 0: Successful exit. = -1: The argument C() is invalid. = -2: The argument Xi() is invalid. = -3: The argument Lxi had an illegal value. (Lxi < 1) = -4: The argument K had an illegal value. (K < 1)

Reference: SLATEC

Example Program

Using the following table, compute S = integral of 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_Ppqad()
    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 Lxi As Long, C(3, Ndata - 2) As Double, Xi(Ndata - 1) As Double
    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
    '-- Convert to PP form
    Call Bsplpp(T(), Bcoef(), N, K, C(), Xi(), Lxi, Info)
    If Info <> 0 Then
        Debug.Print "Error in Bsplpp: Info =", Info
        Exit Sub
    End If
    '-- Compute integral 1/(1 + x^2) dx [0, 4] (= atan(4))
    Call Ppqad(C(), Xi(), Lxi, K, A, B, S, Info)
    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