### 9 One-dimensional Interpolation and Fitting, Splines

The routines for this sort of application (and their origins) are:

• PDA_BSPDOC (SLATEC/CAMSUN)
Documentation for BSPLINE, a package of subprograms for working with piecewise polynomial functions in B-representation.
• PDA_DBINTK (SLATEC/CAMSUN)
Compute the B-representation of a spline which interpolates given data. The knots must be given.
• PDA_DEFC (SLATEC/CAMSUN)
Fit a piecewise polynomial curve to discrete data. The piecewise polynomials are represented as B-splines. The fitting is done in a weighted least squares sense.
• PDA_DBVALU (SLATEC/CAMSUN)
Evaluate the B-representation of a B-spline at X for the function value or any of its derivatives.
Compute the integral of a K-th order B-spline using the B-representation.
• PDA_CURFIT (DIERCKX/NETLIB)
Determine a smooth spline approximation of degree k to the given set of data points. The knots can be given, or can be determined by the routine.
• PDA_SPLEV (DIERCKX/NETLIB)
Evaluates in a number of points x(i) a spline s(x) of degree k, given in its B-spline representation.
• PDA_SPLDER (DIERCKX/NETLIB)
Evaluates in a number of points x(i) the derivative of order NU of a spline s(x) of degree k, given in its B-spline representation.
• PDA_SPLINT (DIERCKX/NETLIB)
Calculates the integral of a spline function s(x) of degree k, which is given in its normalised B-spline representation.
• PDA_DPLINT (SLATEC/CAMSUN)
Produce the polynomial which interpolates a set of discrete data points.
• PDA_DPOLVL (SLATEC/CAMSUN)
Calculate the value of a polynomial and its first NDER derivatives where the polynomial was produced by a previous call to PDA_DPLINT.
• PDA_DPOLCF (SLATEC/CAMSUN)
Compute the coefficients of the polynomial fit (including Hermite polynomial fits) produced by a previous call to PDA_DPLINT.
• PDA_DPOLFT (SLATEC/CAMSUN)
Fit discrete data in a least squares sense by polynomials in one variable. Uses weights.
• PDA_DP1VLU (SLATEC/CAMSUN)
Use the coefficients generated by PDA_DPOLFT to evaluate the polynomial fit of degree L, along with the first NDER of its derivatives, at a specified point.
• PDA_DPCOEF (SLATEC/CAMSUN)
Convert the PDA_DPOLFT coefficients to Taylor series form.

#### 9.1 B-splines

E01BAF finds the interpolating cubic spline interpolant f(x) for a set of points (x,y). The interpolant is evaluated with E02BBF, evaluated with derivatives by E02BCF, and integrated by E02BDF. PDA_DBINTK does this, the order of the splines can be changed as well. This routine needs to be given the knots, while E01BAF set them itself. Evaluation of the interpolant and its derivatives is done by PDA_DBVALU, integration by PDA_DBSQAD.

E02BAF finds the fitting cubic spline f(x) for a set of points and weights (x,y,w). The interior knots 5 ... n+3 must be given and are fixed. The function is evaluated with E02BBF, with derivatives by E02BCF, and integrated with E02BDF. PDA_DEFC does this. All knots must be given, not just the interior ones. Instead of weights PDA_DEFC takes standard deviations (x,y,sigma) and uses 1/sigma as weight. Evaluation of the function and its derivatives is done by PDA_DBVALU, integration by PDA_DBSQAD.

E02BBF and E02BCF evaluate an interpolating or fitting cubic spline and its derivatives. They follow a call to E01BAF or E02BAF. This function is taken over by PDA_DBVALU. E02BDF integrates an interpolating or fitting cubic spline. It follows a call to E01BAF or E02BAF. This function is taken over by PDA_DBSQAD.

E02BEF finds the fitting cubic spline f(x) for a set of points and weights (x,y,w). The knots are located automatically. The function is evaluated with E02BBF, with derivatives by E02BCF, and integrated with E02BDF. PDA_CURFIT (from the DIERCKX package) solves this problem. While the other routines are from SLATEC and for double precision, PDA_CURFIT is for single precision. Hence, PDA_SPLEV, PDA_SPLDER and PDA_SPLINT should be used to evaluate the spline, its n-th derivative, and its integral.

#### 9.2 Ordinary polynomials

E02ADF finds the fitting Chebyshev series minimising r.m.s. The Chebyshev series is equivalent to an ordinary polynomial, but cannot be extrapolated. The polynomial is evaluated by E02AEF or E02AKF. In the latter routine the Chebyshev coefficients can be one column of an array. It will also take the real-world x argument instead of the normalised x-bar argument within the range $-$1 ... +1.

In this library, PDA_DPOLFT fits an ordinary polynomial as a sum of orthogonal polynomials. The representation returned is somewhat special. It can be converted to coefficients of a Taylor series with PDA_DPCOEF or directly evaluated with PDA_DP1VLU. PDA_DP1VLU will return in one call as many derivatives as requested.

#### 9.3 Replacing calls to E01BAF

The SLATEC equivalent of this routine is PDA_DBINTK with order K = 4. The NAG code would look like

INTEGER M, IFAIL
DOUBLE PRECISION X(M), Y(M), T(M+4), C(M+4), WRK(6*M+16)
IFAIL = 1
CALL E01BAF( M, X, Y, T, C, M+4, WRK, 6*M+16, IFAIL )
IF ( IFAIL .NE. 0 ) THEN
An error has occurred
END IF

While the size of the coefficient vector can be reduced to M for PDA_DBINTK, you now need two work spaces. The major difference is that PDA_DBINTK needs to be given the knots. So you have to calculate them in the same way as E01BAF would have done. The handling of the status is different.

INTEGER I, K, M, IFAIL
PARAMETER ( K = 4 )
DOUBLE PRECISION X(M), Y(M), T(M+K), C(M)
DOUBLE PRECISION WRK1( (2*K-1)*M ), WRK2( 2*K )
DO 1 I = 1, K
T(I)   = X(1)
T(M+I) = X(M)
1 CONTINUE
DO 2 I = K+1, M
T(I) = X(I-K/2)      ! Note: K is even
2 CONTINUE
IFAIL = 0
CALL PDA_DBINTK( X, Y, T, M, K, C, WRK1, WRK2, IFAIL )
IF ( IFAIL .NE. 0 ) THEN
An error has occurred
END IF

#### 9.4 Replacing calls to E02BAF

PDA_DEFC has not yet been used anywhere to replace E02BAF. Thus the migration hints given here may contain errors or may be based on misunderstandings.

The SLATEC equivalent of this routine is PDA_DEFC with order K = 4. The NAG code would look like

INTEGER M, N, IFAIL
DOUBLE PRECISION X(M), Y(M), W(M), T(N+7), C(N+7), SS
DOUBLE PRECISION WORK1(M), WORK2( 4*(N+7) )
IFAIL = 1
CALL E02BAF( M, N+7, X, Y, W, T, WORK1, WORK2, C, SS, IFAIL )
IF ( IFAIL .NE. 0 ) THEN
An error has occurred
END IF

Here N+7 is the number of knots. Since the order is 4 (cubic), the number of interior knots is then N+7$-$8 = N$-$1. N is the number of intervals. The interior knots are T(5) ... T(N+3). The weights W are reciprocal errors of Y. Although C has a length of N+7 only N+3 coefficients are returned by E02BAF.

The interior knots T(5) ... T(N+3) are given arguments, but the remaining knots are set by E02BAF and thus returned arguments.

In PDA_DEFC the order is K = 4. There are N+K+3 knots T(). T(1) ... T(K$-$1) and T(N+5) ... T(N+K+3) are end knots. The next inner knots T(K) and T(N+4) are presumably the first and last x value. Then T(K+1) ... T(N+3) would be truly interior knots just as in the NAG code.

PDA_DEFC does not generate knots by itself. Contrary to the NAG code above, the first K and last K knots must be calculated before the call.

The size of the work space is more complex to calculate. PDA_DEFC needs the standard deviation in Y instead of the weights SD = 1/W. PDA_DEFC returns a diagnostic J, which should have value J = 1 if no error occurred.

INTEGER I, J, K, L, M, N, IFAIL
PARAMETER ( K = 4 )
PARAMETER ( L = (N+6) * (K+1)  + (N+K+4) * (K+1)
:              + 2*MAX(M,N+K+3) + N+K+3 + K**2 )
DOUBLE PRECISION X(M), Y(M), W(M), SD(M)
DOUBLE PRECISION T(N+K+3), C(N+3)
DOUBLE PRECISION WORK(L)
DO 1 I = 1, K
T(I)     = MIN( X() )
T(N+3+I) = MAX( X() )
1 CONTINUE
DO 2 I = 1, M
SD(I) = 1D0 / W(I)
2 CONTINUE
IFAIL = 0
CALL PDA_DEFC( M, X, Y, SD, K, N+K+3, T, 1, J, C, L, WORK, IFAIL )
IF ( J .NE. 1 .OR. IFAIL .NE. 0 ) THEN
An error has occurred
END IF

#### 9.5 Replacing calls to E02BBF

The equivalent of this routine in SLATEC is PDA_DBVALU with the requested derivative being zero and the order being K = 4 for a cubic spline. The NAG code would look like

INTEGER N, IFAIL
DOUBLE PRECISION T(N+7), C(N+7), X, S
IFAIL = 1
CALL E02BBF( N+7, T, C, X, S, IFAIL )
IF ( IFAIL .NE. 0 ) THEN
An error has occurred
END IF

Here N, T and C are the same as in E02BAF. If T and C originate from a call to E01BAF then for N+7 read M+4 with M the number of data points given to the interpolation.

PDA_DBVALU is a function rather than a subroutine. The dimension passed to PDA_DBVALU is not that of T, but that of C, i.e. N+3 (or M after interpolation). PDA_DBVALU returns the value of any derivative, the fifth argument is zero so that it returns the function value itself. INVB must be given 1 in the first call. For several evaluations of the same spline it should not be changed between calls. It is changed by PDA_DBVALU. So if PDA_DBVALU is called in a DO loop, the statement INVB = 1 is typically before and outside the loop.

INTEGER INVB, K, N
PARAMETER ( K = 4 )
DOUBLE PRECISION T(N+K+3), C(N+3), X, S
DOUBLE PRECISION WORK(3*K)
DOUBLE PRECISION PDA_DBVALU
INVB = 1
IFAIL = 0
S = PDA_DBVALU( T, C, N+3, K, 0, X, INVB, WORK, IFAIL )
IF ( IFAIL .NE. 0 ) THEN
An error has occurred
END IF

#### 9.6 Replacing calls to E02BCF

E02BCF has not yet been replaced anywhere. Thus the migration hints given here may contain errors or may be based on misunderstandings.

The equivalent of this routine in SLATEC is PDA_DBVALU. Several calls are necessary, one for each derivative. For the function value itself set the number of derivative to zero. The order is K = 4 for a cubic spline. The NAG code would look like

INTEGER LEFT, N, IFAIL
DOUBLE PRECISION T(N+7), C(N+7), X, S(O:3)
IFAIL = 1
CALL E02BCF( N+7, T, C, X, LEFT, S, IFAIL )
IF ( IFAIL .NE. 0 ) THEN
An error has occurred
END IF

Here N, T and C are the same as in E02BAF. If T and C originate from a call to E01BAF then for N+7 read M+4 with M the number of data points given to the interpolation. S(I) returns the I-th derivative. In the case that X coincides with a knot and the derivatives are not continuous at that knot, LEFT is used to decide which side of the knot to use.

PDA_DBVALU is a function rather than a subroutine. The dimension passed to PDA_DBVALU is not that of T, but that of C, i.e. N+3 (or M after interpolation). PDA_DBVALU returns the value of any derivative, as specified in the fifth argument.

There is no equivalent to the LEFT parameter in NAG. PDA_DBVALU returns right limiting values, except at the right end point.

INVB must be given 1 in the first call. For several evaluations of the same spline it should not be changed between calls. It is changed by PDA_DBVALU. So if PDA_DBVALU is called in a DO loop, the statement INVB = 1 is typically before and outside the loop. In the code below, IFAIL is reset inside the DO loop. Assuming that an error will quit the loop, the IFAIL = 0 statement could be before and outside the DO loop as well.

INTEGER INVB, I, K, N
PARAMETER ( K = 4 )
DOUBLE PRECISION T(N+K+3), C(N+3), X, S(0:K-1)
DOUBLE PRECISION WORK(3*K)
DOUBLE PRECISION PDA_DBVALU
INVB = 1
DO 1 I = 0, K-1
IFAIL = 0
S(I) = PDA_DBVALU( T, C, N+3, K, I, X, INVB, WORK, IFAIL )
IF ( IFAIL .NE. 0 ) THEN
An error has occurred
END IF
1 CONTINUE

#### 9.7 Replacing calls to E02BDF

PDA_DBSQAD has not yet been used anywhere to replace E02BDF. Thus the migration hints given here may contain errors or may be based on misunderstandings.

The equivalent of this routine in SLATEC is PDA_DBSQAD. The NAG code would look like

INTEGER N, IFAIL
DOUBLE PRECISION T(N+7), C(N+7), DEFINT
IFAIL = 1
CALL E02BDF( N+7, T, C, DEFINT, IFAIL )
IF ( IFAIL .NE. 0 ) THEN
An error has occurred
END IF

Here N, T and C are the same as in E02BAF. If T and C originate from a call to E01BAF then for N+7 read M+4 with M the number of data points given to the interpolation. DEFINT returns the integral over the whole x range where the spline is defined. This is from T(4) to T(N+4), which are most probably the smallest and largest X used in the fit or interpolation.

The dimension passed to PDA_DBSQAD is not that of T, but that of C, i.e. N+3 (or M after interpolation). PDA_DBSQAD calculates the integral for any interval on which the spline is defined. For the same interval as in the NAG code, the two limiting knots are given to PDA_DBSQAD.

INTEGER K, N
PARAMETER ( K = 4 )
DOUBLE PRECISION T(N+K+3), C(N+3), DEFINT
DOUBLE PRECISION WORK(3*K)
CALL PDA_DBSQAD( T, C, N+3, K, T(4), T(N+4), DEFINT, WORK, IFAIL )
IF ( IFAIL .NE. 0 ) THEN
An error has occurred
END IF

#### 9.8 Replacing calls to E02BEF

E02BEF is a more advanced routine than the other NAG spline routines. It places knots automatically while fitting a cubic spline. This library includes for this case a number of routines from DIERCKX. PDA_CURFIT performs the spline approximation with given or automatic knots. The fitted function is evaluated with PDA_SPLEV, its derivatives with PDA_SPLDER, its integral with PDA_SPLINT. The routines exist only for single precision arguments.

These routines are so far unused, so there are no migration hints.

#### 9.9 Replacing calls to E02ADF

Superficially, the equivalent SLATEC routine is PDA_DPOLFT. Since NAG has no routine to do a polynomial extrapolation, Figaro usurped this routine with a peculiar set of weights and very few degrees of freedom to extrapolate a polynomial. Although PDA_DPOLFT can be used even in that case, it is PDA_DPLINT that is intended for polynomial interpolation.

The NAG code would look like

INTEGER I, M, K, NROWS, IFAIL
DOUBLE PRECISION X(M), Y(M), W(M)
DOUBLE PRECISION WORK1(3*M), WORK2( 2*(K+1) )
DOUBLE PRECISION A1(NROWS,K+1), S(K+1)
DOUBLE PRECISION A2(K+1)
IFAIL = 1
CALL E02ADF( M, K+1, NROWS, X, Y, W, WORK1, WORK2, A1, S, IFAIL )
IF ( IFAIL .NE. 0 ) THEN
An error has occurred
END IF
DO 1 I = 1, K+1
A2(I) = A1(K+1,I)
1 CONTINUE

Here W are the weights proportional to the reciprocal of the standard deviation of Y. A1 returns a matrix of Chebyshev coefficients, one set of coefficients for each polynomial degree from 0 to K. The DO loop extracts the coefficients for degree K into A2. Note that these coefficients form a column rather than a row in A1. S returns the r.m.s. for the fit of each degree from 0 to K.

The behaviour of PDA_DPOLFT is controlled by the given value of EPS, passing zero (0D0) makes it perform fits for all degrees from 0 to K. EPS is also a returned argument, it returns the r.m.s. for the highest degree fitted. What degree that was is returned in NDEG. An indication of the success is returned in IFAIL1.

The weights should be proportional to the reciprocal of the variance, i.e. the square of the weights used for NAG.

The returned description of the polynomials A3 is rather different from the Chebyshev coefficients returned by E02ADF. A3 would be passed on to PDA_DP1VLU to evaluate the polynomial or to PDA_DPCOEF to convert A3 to coefficients of a Taylor series. A3 contains sufficient information to evaluate the polynomial of any degree from 0 to K. The desired degree is specified to PDA_DP1VLU or PDA_DPCOEF.

R is a returned vector in which the fit of highest degree is evaluated at all given X. Often this may render any further evaluation calls obsolete.

INTEGER I, M, K, NDEG, IFAIL1, IFAIL2
DOUBLE PRECISION X(M), Y(M), W(M), W2(M), R(M)
DOUBLE PRECISION A3( 3*M + 3*(K+1) )
DOUBLE PRECISION EPS, S(K+1)
DO 1 I = 1, M
W2(I) = W(I) * W(I)
1 CONTINUE
IFAIL2 = 0
EPS = 0D0
CALL PDA_DPOLFT( M, X, Y, W2, K, NDEG, EPS, R, IFAIL1, A3, IFAIL2 )
IF ( NDEG .NE. K .OR. IFAIL1 .NE. 1 .OR. IFAIL2 .NE. 0 ) THEN
An error has occurred
END IF
DO 2 I = 1, K
S(I) = 0D0
2 CONTINUE
S(K+1) = EPS

#### 9.10 Replacing calls to E02AEF

The equivalent of this routine in SLATEC is PDA_DP1VLU. The NAG code would look like

INTEGER I, K, K2, NROWS, IFAIL
DOUBLE PRECISION A1(NROWS,K+1)
DOUBLE PRECISION A2(K2+1), XX, XCAP, P
DO 1 I = 1, K2+1
A2(I) = A1(K2+1,I)
1 CONTINUE
XCAP = ( ( XX - MIN(X()) ) - ( MAX(X()) - XX ) )
:     / ( MAX(X()) - MIN(X()) )
IFAIL = 1
CALL E02AEF( K2+1, A2, XCAP, P, IFAIL )
IF ( IFAIL .NE. 0 ) THEN
An error has occurred
END IF

Here A1 is the coefficient matrix returned by E02ADF, A2 is the column extracted for the required degree. XX is the x value for which an evaluation is required, it must be scaled into the range $-$1 ... +1, using the original extreme x values passed to E02ADF. P returns the function value.

PDA_DP1VLU returns any number of derivatives in addition to the function value. The second argument specifies how many derivatives are required. No scaling of XX is necessary, and no processing of the coefficients A3 as returned by PDA_DPOLFT.

INTEGER K, K2, IFAIL
DOUBLE PRECISION A3( 3*M + 3*(K+1) )
DOUBLE PRECISION XX, P, DUMMY
IFAIL = 0
CALL PDA_DP1VLU( K2, 0, XX, P, DUMMY, A3, IFAIL )
IF ( IFAIL .NE. 0 ) THEN
An error has occurred
END IF

#### 9.11 Replacing calls to E02AKF

E02AKF is similar to E02AEF but provides a different interface. Namely, A1 (from E02ADF) and the row length can be given instead of the extraction A2. Also XX is given rather than the scaled XCAP. There are no detailed migration hints for this routine yet.

#### 9.12 Replacing calls to GEN_CHB2NO

This is not a NAG routine, but a routine in Figaro’s GEN library. It converts the vector of Chebyshev coefficients into a vector of coefficients of an ordinary polynomial. This is particularly useful when these are to be written to output files. In future the routine will still be needed to read old files that contain Chebyshev coefficients.

This routine is mentioned here really, because it has an equivalent in SLATEC named PDA_DPCOEF. PDA_DPCOEF is more general, in that it returns coefficients for a Taylor series, i.e. a polynomial in (XX-X0) for given expansion point X0.

The old code would look like

INTEGER I, K, K2, NROWS
DOUBLE PRECISION A1(NROWS,K+1)
DOUBLE PRECISION A2(K2+1), C(K2+1)
DO 1 I = 1, K2+1
A2(I) = A1(K2+1,I)
1 CONTINUE
CALL GEN_CHB2NO( K2, MIN(X()), MAX(X()), A2, C )

This would be replaced by

INTEGER K, K2, IFAIL
DOUBLE PRECISION A3( 3*M + 3*(K+1) )
DOUBLE PRECISION C(K2+1)
IFAIL = 0
CALL PDA_DPCOEF( K2, 0D0, C, A3, IFAIL )
IF ( IFAIL .NE. 0 ) THEN
An error has occurred
END IF

Just in case you are tempted to evaluate the polynomial yourself rather than use PDA_DP1VLU, here is a piece of code to get P = f(XX):

INTEGER I, K2
DOUBLE PRECISION XX, P
DOUBLE PRECISION C(K2+1)
P = C(K2+1)
DO 2 I = K2, 1, -1
P = P * XX
P = P + C(I)
2 CONTINUE