BSLS2
Computes a two-dimensional tensor-product spline approximant using least squares, returning the tensor-product B-spline coefficients.
Required Arguments
XDATA — Array of length NXDATA containing the data points in the X-direction. (Input)
XDATA must be nondecreasing.
YDATA — Array of length NYDATA containing the data points in the Y-direction. (Input)
YDATA must be nondecreasing.
FDATA — Array of size NXDATA by NYDATA containing the values on the X − Y grid to be interpolated. (Input)
FDATA(I, J) contains the value at (XDATA(I), YDATA(I)).
KXORD — Order of the spline in the X-direction. (Input)
KYORD — Order of the spline in the Y-direction. (Input)
XKNOT — Array of length KXORD + NXCOEF containing the knots in the X-direction. (Input)
XKNOT must be nondecreasing.
YKNOT — Array of length KYORD + NYCOEF containing the knots in the Y-direction. (Input)
YKNOT must be nondecreasing.
BSCOEF — Array of length NXCOEF * NYCOEF that contains the tensor product B-spline coefficients. (Output)
BSCOEF is treated internally as an array of size NXCOEF by NYCOEF.
Optional Arguments
NXDATA — Number of data points in the X-direction. (Input)
Default: NXDATA = size (XDATA,1).
NYDATA — Number of data points in the Y-direction. (Input)
Default: NYDATA = size (YDATA,1).
LDF — Leading dimension of FDATA exactly as specified in the dimension statement of calling program. (Input)
Default: LDF = size (FDATA,1).
NXCOEF — Number of B-spline coefficients in the X-direction. (Input)
Default: NXCOEF = size (XKNOT,1) - KXORD.
NYCOEF — Number of B-spline coefficients in the Y-direction. (Input)
Default: NYCOEF = size (YKNOT,1) - KYORD.
XWEIGH — Array of length NXDATA containing the positive weights of XDATA. (Input)
Default: XWEIGH = 1.0.
YWEIGH — Array of length NYDATA containing the positive weights of YDATA. (Input)
Default: YWEIGH = 1.0.
FORTRAN 90 Interface
Generic: CALL BSLS2 (XDATA, YDATA, FDATA, KXORD, KYORD, XKNOT, YKNOT, BSCOEF [, …])
Specific: The specific interface names are S_BSLS2 and D_BSLS2.
FORTRAN 77 Interface
Single: CALL BSLS2 (NXDATA, XDATA, NYDATA, YDATA, FDATA, LDF, KXORD, KYORD, XKNOT, YKNOT, NXCOEF, NYCOEF, XWEIGH, YWEIGH, BSCOEF)
Double: The double precision name is DBSLS2.
Description
The routine BSLS2 computes the coefficients of a tensor-product spline least-squares approximation to weighted tensor-product data. The input for this subroutine consists of data vectors to specify the tensor-product grid for the data, two vectors with the weights, the values of the surface on the grid, and the specification for the tensor-product spline. The grid is specified by the two vectors x = XDATA and y = YDATA of length n = NXDATA and m = NYDATA, respectively. A two-dimensional array f = FDATA contains the data values that are to be fit. The two vectors wx = XWEIGH and wy = YWEIGH contain the weights for the weighted least-squares problem. The information for the approximating tensor-product spline must also be provided. This information is contained in kx = KXORD, tx = XKNOT, and N = NXCOEF for the spline in the first variable, and in ky = KYORD, ty = YKNOT and M = NYCOEF for the spline in the second variable. The coefficients of the resulting tensor-product spline are returned in c = BSCOEF, which is an N * M array. The procedure computes coefficients by solving the normal equations in tensor-product form as discussed in de Boor (1978, Chapter 17). The interested reader might also want to study the paper by E. Grosse (1980).
The final result produces coefficients c minimizing
where the function Bkl is the tensor-product of two B-splines of order kx and ky. Specifically, we have
The spline
can be evaluated using
BS2VL and its partial derivatives can be evaluated using
BS2DR.
Comments
1. Workspace may be explicitly provided, if desired, by use of B2LS2/DB2LS2. The reference is:
CALL B2LS2 (NXDATA, XDATA, NYDATA, YDATA, FDATA, LDF, KXORD, KYORD, XKNOT, YKNOT, NXCOEF, NYCOEF, XWEIGH, YWEIGH, BSCOEF, WK)
The additional argument is:
WK — Work array of length (NXCOEF + 1) * NYDATA + KXORD * NXCOEF + KYORD * NYCOEF + 3 * MAX(KXORD, KYORD).
2. Informational errors
Type | Code | Description |
|---|
3 | 14 | There may be less than one digit of accuracy in the least squares fit. Try using higher precision if possible. |
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 | 9 | The data point abscissae must be nondecreasing. |
4 | 10 | The smallest element of the data point array must be greater than or equal to the K_ORDth knot. |
4 | 11 | The largest element of the data point array must be less than or equal to the (N_COEF + 1)st knot. |
Example
The data for this example arise from the function ex sin(x + y) + ɛ on the rectangle [0, 3] × [0, 5]. Here, ɛ is a uniform random variable with range [−1, 1]. We sample this function on a 100 × 50 grid and then try to recover it by using cubic splines in the x variable and quadratic splines in the y variable. We print out the values of the function ex sin(x + y) on a 3 × 5 grid and compare these values with the values of the tensor-product spline that was computed using the IMSL routine BSLS2.
USE IMSL_LIBRARIES
IMPLICIT NONE
INTEGER KXORD, KYORD, LDF, NXCOEF, NXDATA, NXVEC, NYCOEF,&
NYDATA, NYVEC
PARAMETER (KXORD=4, KYORD=3, NXCOEF=15, NXDATA=100, NXVEC=4,&
NYCOEF=7, NYDATA=50, NYVEC=6, LDF=NXDATA)
!
INTEGER I, J, NOUT
REAL BSCOEF(NXCOEF,NYCOEF), EXP, F, FDATA(NXDATA,NYDATA),&
FLOAT, RNOISE, SIN, VALUE(NXVEC,NYVEC), X,&
XDATA(NXDATA), XKNOT(NXCOEF+KXORD), XVEC(NXVEC),&
XWEIGH(NXDATA), Y, YDATA(NYDATA),&
YKNOT(NYCOEF+KYORD), YVEC(NYVEC), YWEIGH(NYDATA)
INTRINSIC EXP, FLOAT, SIN
! Define function
F(X,Y) = EXP(X)*SIN(X+Y)
! Set random number seed
CALL RNSET (1234579)
! Set up X knot sequence.
DO 10 I=1, NXCOEF - KXORD + 2
XKNOT(I+KXORD-1) = 3.0*(FLOAT(I-1)/FLOAT(NXCOEF-KXORD+1))
10 CONTINUE
XKNOT(NXCOEF+1) = XKNOT(NXCOEF+1) + 0.001
! Stack knots.
DO 20 I=1, KXORD - 1
XKNOT(I) = XKNOT(KXORD)
XKNOT(I+NXCOEF+1) = XKNOT(NXCOEF+1)
20 CONTINUE
! Set up Y knot sequence.
DO 30 I=1, NYCOEF - KYORD + 2
YKNOT(I+KYORD-1) = 5.0*(FLOAT(I-1)/FLOAT(NYCOEF-KYORD+1))
30 CONTINUE
YKNOT(NYCOEF+1) = YKNOT(NYCOEF+1) + 0.001
! Stack knots.
DO 40 I=1, KYORD - 1
YKNOT(I) = YKNOT(KYORD)
YKNOT(I+NYCOEF+1) = YKNOT(NYCOEF+1)
40 CONTINUE
! Set up X-grid.
DO 50 I=1, NXDATA
XDATA(I) = 3.0*(FLOAT(I-1)/FLOAT(NXDATA-1))
50 CONTINUE
! Set up Y-grid.
DO 60 I=1, NYDATA
YDATA(I) = 5.0*(FLOAT(I-1)/FLOAT(NYDATA-1))
60 CONTINUE
! Evaluate function on grid and
! introduce random noise in [1,-1].
DO 70 I=1, NYDATA
DO 70 J=1, NXDATA
RNOISE = RNUNF()
RNOISE = 2.0*RNOISE - 1.0
FDATA(J,I) = F(XDATA(J),YDATA(I)) + RNOISE
70 CONTINUE
! Use default weights equal to 1.
!
! Compute least squares approximation.
CALL BSLS2 (XDATA, YDATA, FDATA, KXORD, KYORD, &
XKNOT, YKNOT, BSCOEF)
! Get output unit number
CALL UMACH (2, NOUT)
! Write heading
WRITE (NOUT,99999)
! Print interpolated values
! on [0,3] x [0,5].
DO 80 I=1, NXVEC
XVEC(I) = FLOAT(I-1)
80 CONTINUE
DO 90 I=1, NYVEC
YVEC(I) = FLOAT(I-1)
90 CONTINUE
! Evaluate spline
CALL BS2GD (0, 0, XVEC, YVEC, KXORD, KYORD, XKNOT,&
YKNOT, BSCOEF, VALUE)
DO 110 I=1, NXVEC
DO 100 J=1, NYVEC
WRITE (NOUT,'(5F15.4)') XVEC(I), YVEC(J),&
F(XVEC(I),YVEC(J)), VALUE(I,J),&
(F(XVEC(I),YVEC(J))-VALUE(I,J))
100 CONTINUE
110 CONTINUE
99999 FORMAT (13X, 'X', 14X, 'Y', 10X, 'F(X,Y)', 9X, 'S(X,Y)', 10X,&
'Error')
END
Output
X Y F(X,Y) S(X,Y) Error
0.0000 0.0000 0.0000 0.2782 -0.2782
0.0000 1.0000 0.8415 0.7762 0.0653
0.0000 2.0000 0.9093 0.8203 0.0890
0.0000 3.0000 0.1411 0.1391 0.0020
0.0000 4.0000 -0.7568 -0.5705 -0.1863
0.0000 5.0000 -0.9589 -1.0290 0.0701
1.0000 0.0000 2.2874 2.2678 0.0196
1.0000 1.0000 2.4717 2.4490 0.0227
1.0000 2.0000 0.3836 0.4947 -0.1111
1.0000 3.0000 -2.0572 -2.0378 -0.0195
1.0000 4.0000 -2.6066 -2.6218 0.0151
1.0000 5.0000 -0.7595 -0.7274 -0.0321
2.0000 0.0000 6.7188 6.6923 0.0265
2.0000 1.0000 1.0427 0.8492 0.1935
2.0000 2.0000 -5.5921 -5.5885 -0.0035
2.0000 3.0000 -7.0855 -7.0955 0.0099
2.0000 4.0000 -2.0646 -2.1588 0.0942
2.0000 5.0000 4.8545 4.7339 0.1206
3.0000 0.0000 2.8345 2.5971 0.2373
3.0000 1.0000 -15.2008 -15.1079 -0.0929
3.0000 2.0000 -19.2605 -19.1698 -0.0907
3.0000 3.0000 -5.6122 -5.5820 -0.0302
3.0000 4.0000 13.1959 12.6659 0.5300
3.0000 5.0000 19.8718 20.5170 -0.6452
