BSLSQ
Computes the least-squares spline approximation, and returns the B-spline coefficients.
Required Arguments
XDATA — Array of length NDATA containing the data point abscissas. (Input)
FDATA — Array of length NDATA containing the data point ordinates. (Input)
KORDER — Order of the spline. (Input)
KORDER must be less than or equal to NDATA.
XKNOT — Array of length NCOEF + KORDER containing the knot sequence. (Input)
XKNOT must be nondecreasing.
NCOEF — Number of B-spline coefficients. (Input)
NCOEF cannot be greater than NDATA.
BSCOEF — Array of length NCOEF containing the B-spline coefficients. (Output)
Optional Arguments
NDATA — Number of data points. (Input)
Default: NDATA = size(XDATA, 1)
WEIGHT — Array of length NDATA containing the weights. (Input)
Default: WEIGHT = 1.0.
FORTRAN 90 Interface
Generic: CALL BSLSQ (XDATA, FDATA, KORDER, XKNOT, NCOEF, BSCOEF [, …])
Specific: The specific interface names are S_BSLSQ and D_BSLSQ.
FORTRAN 77 Interface
Single: CALL BSLSQ (NDATA, XDATA, FDATA, WEIGHT, KORDER, XKNOT, NCOEF, BSCOEF)
Double: The double precision name is DBSLSQ.
Description
The routine BSLSQ is based on the routine L2APPR by de Boor (1978, page 255). The IMSL routine BSLSQ computes a weighted discrete L2 approximation from a spline subspace to a given data set (xi, fi) for i = 1, …, N (where N = NDATA). In other words, it finds B-spline coefficients, a = BSCOEF, such that
is a minimum, where m = NCOEF and Bj denotes the j-th B-spline for the given order, KORDER, and knot sequence, XKNOT. This linear least squares problem is solved by computing and solving the normal equations. While the normal equations can sometimes cause numerical difficulties, their use here should not cause a problem because the B-spline basis generally leads to well-conditioned banded matrices.
The choice of weights depends on the problem. In some cases, there is a natural choice for the weights based on the relative importance of the data points. To approximate a continuous function (if the location of the data points can be chosen), then the use of Gauss quadrature weights and points is reasonable. This follows because BSLSQ is minimizing an approximation to the integral
The Gauss quadrature weights and points can be obtained using the IMSL routine GQRUL (see Chapter 4, “Integration and Differentiation”).
Comments
1. Workspace may be explicitly provided, if desired, by use of B2LSQ/DB2LSQ. The reference is:
CALL B2LSQ (NDATA, XDATA, FDATA, WEIGHT, KORDER, XKNOT, NCOEF, BSCOEF, WK1, WK2, WK3, WK4, IWK)
The additional arguments are as follows:
WK1 — Work array of length (3 + NCOEF) * KORDER.
WK2 — Work array of length NDATA.
WK3 — Work array of length NDATA.
WK4 — Work array of length NDATA.
IWK — Work array of length NDATA.
2. Informational errors
|
Type |
Code |
Description |
|
4 |
5 |
Multiplicity of the knots cannot exceed the order of the spline. |
|
4 |
6 |
The knots must be nondecreasing. |
|
4 |
7 |
All weights must be greater than zero. |
|
4 |
8 |
The smallest element of the data point array must be greater than or equal to the KORDth knot. |
|
4 |
9 |
The largest element of the data point array must be less than or equal to the (NCOEF + 1)st knot. |
3. The B-spline representation can be evaluated using BSVAL, and its derivative can be evaluated using BSDER.
Example
In this example, we try to recover a quadratic polynomial using a quadratic spline with one interior knot from two different data sets. The first data set is generated by evaluating the quadratic at 50 equally spaced points in the interval (0, 1) and then adding uniformly distributed noise to the data. The second data set includes the first data set, and, additionally, the values at 0 and at 1 with no noise added. Since the first and last data points are uncontaminated by noise, we have chosen weights equal to 105 for these two points in this second problem. The quadratic, the first approximation, and the second approximation are then evaluated at 11 equally spaced points. This example illustrates the use of the weights to enforce interpolation at certain of the data points.
USE IMSL_LIBRARIES
IMPLICIT NONE
INTEGER KORDER, NCOEF
PARAMETER (KORDER=3, NCOEF=4)
!
INTEGER I, NDATA, NOUT
REAL ABS, BSCOF1(NCOEF), BSCOF2(NCOEF), F,&
FDATA1(50), FDATA2(52), FLOAT, RNOISE, S1,&
S2, WEIGHT(52), X, XDATA1(50), XDATA2(52),&
XKNOT(KORDER+NCOEF), XT, YT
INTRINSIC ABS, FLOAT
!
DATA WEIGHT/52*1.0/
! Define function
F(X) = 8.0*X*(1.0-X)
! Set random number seed
CALL RNSET (12345679)
NDATA = 50
! Set up interior knots
DO 10 I=1, NCOEF - KORDER + 2
XKNOT(I+KORDER-1) = FLOAT(I-1)/FLOAT(NCOEF-KORDER+1)
10 CONTINUE
! Stack knots
DO 20 I=1, KORDER - 1
XKNOT(I) = XKNOT(KORDER)
XKNOT(I+NCOEF+1) = XKNOT(NCOEF+1)
20 CONTINUE
! Set up data points excluding
! the endpoints 0 and 1.
! The function values have noise
! introduced.
DO 30 I=1, NDATA
XDATA1(I) = FLOAT(I)/51.0
RNOISE = RNUNF()
RNOISE = RNOISE - 0.5
FDATA1(I) = F(XDATA1(I)) + RNOISE
30 CONTINUE
! Compute least squares B-spline
! representation.
CALL BSLSQ (XDATA1, FDATA1, KORDER, XKNOT, NCOEF, BSCOF1)
! Now use same XDATA values but with
! the endpoints included. These
! points will have large weights.
NDATA = 52
CALL SCOPY (50, XDATA1, 1, XDATA2(2:), 1)
CALL SCOPY (50, FDATA1, 1, FDATA2(2:), 1)
!
WEIGHT(1) = 1.0E5
XDATA2(1) = 0.0
FDATA2(1) = F(XDATA2(1))
WEIGHT(NDATA) = 1.0E5
XDATA2(NDATA) = 1.0
FDATA2(NDATA) = F(XDATA2(NDATA))
! Compute least squares B-spline
! representation.
CALL BSLSQ (XDATA2, FDATA2, KORDER, XKNOT, NCOEF, BSCOF2, &
WEIGHT=WEIGHT)
! Get output unit number
CALL UMACH (2, NOUT)
! Write heading
WRITE (NOUT,99998)
! Print the two interpolants
! at 11 points.
DO 40 I=1, 11
XT = FLOAT(I-1)/10.0
YT = F(XT)
! Evaluate splines
S1 = BSVAL(XT,KORDER,XKNOT,NCOEF,BSCOF1)
S2 = BSVAL(XT,KORDER,XKNOT,NCOEF,BSCOF2)
WRITE (NOUT,99999) XT, YT, S1, S2, (S1-YT), (S2-YT)
40 CONTINUE
!
99998 FORMAT (7X, 'X', 9X, 'F(X)', 6X, 'S1(X)', 5X, 'S2(X)', 7X,&
'F(X)-S1(X)', 7X, 'F(X)-S2(X)')
99999 FORMAT (' ', 4F10.4, 4X, F10.4, 7X, F10.4)
END
X F(X) S1(X) S2(X) F(X)-S1(X) F(X)-S2(X)
0.0000 0.0000 0.0515 0.0000 0.0515 0.0000
0.1000 0.7200 0.7594 0.7490 0.0394 0.0290
0.2000 1.2800 1.3142 1.3277 0.0342 0.0477
0.3000 1.6800 1.7158 1.7362 0.0358 0.0562
0.4000 1.9200 1.9641 1.9744 0.0441 0.0544
0.5000 2.0000 2.0593 2.0423 0.0593 0.0423
0.6000 1.9200 1.9842 1.9468 0.0642 0.0268
0.7000 1.6800 1.7220 1.6948 0.0420 0.0148
0.8000 1.2800 1.2726 1.2863 -0.0074 0.0063
0.9000 0.7200 0.6360 0.7214 -0.0840 0.0014
1.0000 0.0000 -0.1878 0.0000 -0.1878 0.0000
