SPLINE_FITTING

FNLMath : Interpolation and Approximation : SPLINE_FITTING

SPLINE_FITTING

Weighted least-squares fitting by B-splines to discrete One-Dimensional data is performed. Constraints on the spline or its derivatives are optional. The spline function

its derivatives, or the square root of its variance function are evaluated after the fitting.

Required Arguments

DATA = DATA(1:3,:) (Input/Output)
An assumed-shape array with size(data,1) = 3. The data are placed in the array: data(1,i) = xi , data(2,i) = yi, and data(3,i) = σi, i = 1,…,ndata. If the variances are not known but are proportional to an unknown value, users may set data(3,i) = 1, i = 1,…,ndata.

KNOTS = KNOTS (Input)
A derived type, ?_spline_knots, that defines the degree of the spline and the breakpoints for the data fitting interval.

Optional Arguments

CONSTRAINTS = SPLINE_CONSTRAINTS (Input)
A rank-1 array of derived type ?_spline_constraints that give constraints the output spline is to satisfy.

COVARIANCE = G (Output)
An assumed-shape rank-2 array of the same precision as the data. This output is the covariance matrix of the coefficients. It is optionally used to evaluate the square root of the variance function.

IOPT = IOPT(:) (Input/Output)
Derived type array with the same precision as the input array; used for passing optional data to SPLINE_FITTING. The options are as follows:

Packaged Options for SPLINE_FITTING
Prefix = None	Option Name	Option Value
	SPLINE_FITTING_TOL_EQUAL	1
	SPLINE_FITTING_TOL_LEAST	2

IOPT(IO) = ?_OPTIONS(SPLINE_FITTING_TOL_EQUAL, ?_VALUE)
This resets the value for determining that equality constraint equations are rank-deficient. The default is ?_value = 10-4.

IOPT(IO) = ?_OPTIONS(SPLINE_FITTING_TOL_LEAST, ?_VALUE)
This resets the value for determining that least-squares equations are rank-deficient. The default is ?_value = 10-4.

FORTRAN 90 Interface

Generic: CALL SPLINE_FITTING (DATA, KNOTS [, …])

Specific: The specific interface names are S_SPLINE_FITTING and D_SPLINE_FITTING.

Description

This routine has similar scope to CONFT found in IMSL (2003, pp 734-743). We provide the square root of the variance function, but we do not provide for constraints on the integral of the spline. The least-squares matrix problem for the coefficients is banded, with band-width equal to the spline order. This fact is used to obtain an efficient solution algorithm when there are no constraints. When constraints are present the routine solves a linear-least squares problem with equality and inequality constraints. The processed least-squares equations result in a banded and upper triangular matrix, following accumulation of the spline fitting equations. The algorithm used for solving the constrained least-squares system will handle rank-deficient problems. A set of reference are available in Hanson (1995) and Lawson and Hanson (1995). The CONFT routine uses QPROG (loc cit., p. 959), which requires that the least-squares equations be of full rank.

Fatal and Terminal Error Messages

See the messages.gls file for error messages for SPLINE_FITTING. These error messages are numbered 1340–1367.

Examples

Example 1: Natural Cubic Spline Interpolation to Data

The function

is interpolated by cubic splines on the grid of points

Those natural conditions are

Our program checks the term const. appearing in the maximum truncation error term

error ≈ const.×Δx4

at a finer grid.

USE spline_fitting_int

USE show_int

USE norm_int

implicit none

! This is Example 1 for SPLINE_FITTING, Natural Spline

! Interpolation using cubic splines. Use the function

! exp(-x**2/2) to generate samples.

integer :: i

integer, parameter :: ndata=24, nord=4, ndegree=nord-1, &

nbkpt=ndata+2*ndegree, ncoeff=nbkpt-nord, nvalues=2*ndata

real(kind(1e0)), parameter :: zero=0e0, one=1e0, half=5e-1

real(kind(1e0)), parameter :: delta_x=0.15, delta_xv=0.4*delta_x

real(kind(1e0)), target :: xdata(ndata), ydata(ndata), &

spline_data (3, ndata), bkpt(nbkpt), &

ycheck(nvalues), coeff(ncoeff), &

xvalues(nvalues), yvalues(nvalues), diffs

real(kind(1e0)), pointer :: pointer_bkpt(:)

type (s_spline_knots) break_points

type (s_spline_constraints) constraints(2)

xdata = (/((i-1)*delta_x, i=1,ndata)/)

ydata = exp(-half*xdata**2)

xvalues =(/(0.03+(i-1)*delta_xv,i=1,nvalues)/)

ycheck= exp(-half*xvalues**2)

spline_data(1,:)=xdata

spline_data(2,:)=ydata

spline_data(3,:)=one

! Define the knots for the interpolation problem.

bkpt(1:ndegree) = (/(i*delta_x, i=-ndegree,-1)/)

bkpt(nord:nbkpt-ndegree) = xdata

bkpt(nbkpt-ndegree+1:nbkpt) = &

(/(xdata(ndata)+i*delta_x, i=1,ndegree)/)

! Assign the degree of the polynomial and the knots.

pointer_bkpt => bkpt

break_points=s_spline_knots(ndegree, pointer_bkpt)

! These are the natural conditions for interpolating cubic

! splines. The derivatives match those of the interpolating

! function at the ends.

constraints(1)=spline_constraints &

(derivative=2, point=bkpt(nord), type='==', value=-one)

constraints(2)=spline_constraints &

(derivative=2,point=bkpt(nbkpt-ndegree), type= '==', &

value=(-one+xdata(ndata)**2)*ydata(ndata))

coeff = spline_fitting(data=spline_data, knots=break_points,&

constraints=constraints)

yvalues=spline_values(0, xvalues, break_points, coeff)

diffs=norm(yvalues-ycheck,huge(1))/delta_x**nord

if (diffs <= one) then

write(*,*) 'Example 1 for SPLINE_FITTING is correct.'

end if

end

Output

Example 1 for SPLINE_FITTING is correct.

Example 2: Shaping a Curve and its Derivatives

The function

is fit by cubic splines on the grid of equally spaced points

The term noise is uniform random numbers from the normalized interval [-

] where

= 0.01. The spline curve is constrained to be convex down for 0 ≤ x ≤ 1 convex upward for 1 < x ≤ 4, and have the second derivative exactly equal to the value zero at x = 1. The first derivative is constrained with the value zero at x = 0 and is non-negative at the right and of the interval, x = 4. A sample table of independent variables, second derivatives and square root of variance function values is printed.

use spline_fitting_int

use show_int

use rand_int

use norm_int

implicit none

! This is Example 2 for SPLINE_FITTING. Use 1st and 2nd derivative

! constraints to shape the splines.

integer :: i, icurv

integer, parameter :: nbkptin=13, nord=4, ndegree=nord-1, &

nbkpt=nbkptin+2*ndegree, ndata=21, ncoeff=nbkpt-nord

real(kind(1e0)), parameter :: zero=0e0, one=1e0, half=5e-1

real(kind(1e0)), parameter :: range=4.0, ratio=0.02, tol=ratio*half

real(kind(1e0)), parameter :: delta_x=range/(ndata-1), &

delta_b=range/(nbkptin-1)

real(kind(1e0)), target :: xdata(ndata), ydata(ndata), ynoise(ndata),&

sddata(ndata), spline_data (3, ndata), bkpt(nbkpt), &

values(ndata), derivat1(ndata), derivat2(ndata), &

coeff(ncoeff), root_variance(ndata), diffs

real(kind(1e0)), dimension(ncoeff,ncoeff) :: sigma_squared

real(kind(1e0)), pointer :: pointer_bkpt(:)

type (s_spline_knots) break_points

type (s_spline_constraints) constraints(nbkptin+2)

xdata = (/((i-1)*delta_x, i=1,ndata)/)

ydata = exp(-half*xdata**2)

ynoise = ratio*ydata*(rand(ynoise)-half)

ydata = ydata+ynoise

sddata = ynoise

spline_data(1,:)=xdata

spline_data(2,:)=ydata

spline_data(3,:)=sddata

bkpt=(/((i-nord)*delta_b, i=1,nbkpt)/)

! Assign the degree of the polynomial and the knots.

pointer_bkpt => bkpt

break_points=s_spline_knots(ndegree, pointer_bkpt)

icurv=int(one/delta_b)+1

! At first shape the curve to be convex down.

do i=1,icurv-1

constraints(i)=spline_constraints &

(derivative=2, point=bkpt(i+ndegree), type='<=', value=zero)

end do

! Force a curvature change.

constraints(icurv)=spline_constraints &

(derivative=2, point=bkpt(icurv+ndegree), type='==', value=zero)

! Finally, shape the curve to be convex up.

do i=icurv+1,nbkptin

constraints(i)=spline_constraints &

(derivative=2, point=bkpt(i+ndegree), type='>=', value=zero)

end do

! Make the slope zero and value non-negative at right.

constraints(nbkptin+1)=spline_constraints &

(derivative=1, point=bkpt(nord), type='==', value=zero)

constraints(nbkptin+2)=spline_constraints &

(derivative=0, point=bkpt(nbkptin+ndegree), type='>=', value=zero)

coeff = spline_fitting(data=spline_data, knots=break_points, &

constraints=constraints, covariance=sigma_squared)

! Compute value, first two derivatives and the variance.

values=spline_values(0, xdata, break_points, coeff)

root_variance=spline_values(0, xdata, break_points, coeff, &

covariance=sigma_squared)

derivat1=spline_values(1, xdata, break_points, coeff)

derivat2=spline_values(2, xdata, break_points, coeff)

call show(reshape((/xdata, derivat2, root_variance/),(/ndata,3/)),&

"The x values, 2-nd derivatives, and square root of variance.")

! See that differences are relatively small and the curve has

! the right shape and signs.

diffs=norm(values-ydata)/norm(ydata)

if (all(values > zero) .and. all(derivat1 < epsilon(zero))&

.and. diffs <= tol) then

write(*,*) 'Example 2 for SPLINE_FITTING is correct.'

end if

end

Output

Example 2 for SPLINE_FITTING is correct.

Example 3: Splines Model a Random Number Generator

The function

is an unnormalized probability distribution. This function is similar to the standard Normal distribution, with specific choices for the mean and variance, except that it is truncated. Our algorithm interpolates g(x) with a natural cubic spline, f(x). The cumulative distribution is approximated by precise evaluation of the function

Gauss-Legendre quadrature formulas, IMSL (1994, pp. 621-626), of order two are used on each polynomial piece of f(t) to evaluate q(x) cheaply. After normalizing the cubic spline so that q(1) = 1, we may then generate random numbers according to the distribution f(x ≅ g(x. The values of x are evaluated by solving q(x) = u, ‑1 < x < 1. Here u is a uniform random sample. Newton’s method, for a vector of unknowns, is used for the solution algorithm. Recalling the relation

we believe this illustrates a method for generating a vector of random numbers according to a continuous distribution function having finite support.

use spline_fitting_int

use linear_operators

use Numerical_Libraries

implicit none

! This is Example 3 for SPLINE_FITTING. Use splines to

! generate random (almost normal) numbers. The normal distribution

! function has support (-1,+1), and is zero outside this interval.

! The variance is 0.5.

integer i, niterat

integer, parameter :: iweight=1, nfix=0, nord=4, ndata=50

integer, parameter :: nquad=(nord+1)/2, ndegree=nord-1

integer, parameter :: nbkpt=ndata+2*ndegree, ncoeff=nbkpt-nord

integer, parameter :: last=nbkpt-ndegree, n_samples=1000

integer, parameter :: limit=10

real(kind(1e0)), dimension(n_samples) :: fn, rn, x, alpha_x, beta_x

INTEGER LEFT_OF(n_samples)

real(kind(1e0)), parameter :: one=1e0, half=5e-1, zero=0e0, two=2e0

real(kind(1e0)), parameter :: delta_x=two/(ndata-1)

real(kind(1e0)), parameter :: qalpha=zero, qbeta=zero, domain=two

real(kind(1e0)) qx(nquad), qxi(nquad), qw(nquad), qxfix(nquad)

real(kind(1e0)) alpha_, beta_, quad(0:ndata-1)

real(kind(1e0)), target :: xdata(ndata), ydata(ndata),&

coeff(ncoeff), spline_data(3, ndata), bkpt(nbkpt)

real(kind(1e0)), pointer :: pointer_bkpt(:)

type (s_spline_knots) break_points

type (s_spline_constraints) constraints(2)

! Approximate the probability density function by splines.

xdata = (/(-one+(i-1)*delta_x, i=1,ndata)/)

ydata = exp(-half*xdata**2)

spline_data(1,:)=xdata

spline_data(2,:)=ydata

spline_data(3,:)=one

bkpt=(/(-one+(i-nord)*delta_x, i=1,nbkpt)/)

! Assign the degree of the polynomial and the knots.

pointer_bkpt => bkpt

break_points=s_spline_knots(ndegree, pointer_bkpt)

! Define the natural derivatives constraints:

constraints(1)=spline_constraints &

(derivative=2, point=bkpt(nord), type='==', &

value=(-one+xdata(1)**2)*ydata(1))

constraints(2)=spline_constraints &

(derivative=2, point=bkpt(last), type='==', &

value=(-one+xdata(ndata)**2)*ydata(ndata))

! Obtain the spline coefficients.

coeff=spline_fitting(data=spline_data, knots=break_points,&

constraints=constraints)

! Compute the evaluation points 'qx(*)' and weights 'qw(*)' for

! the Gauss-Legendre quadrature. This will give a precise

! quadrature for polynomials of degree <= nquad*2.

call gqrul(nquad, iweight, qalpha, qbeta, nfix, qxfix, qx, qw)

! Compute pieces of the accumulated distribution function:

quad(0)=zero

do i=1, ndata-1

alpha_= (bkpt(nord+i)-bkpt(ndegree+i))*half

beta_ = (bkpt(nord+i)+bkpt(ndegree+i))*half

! Normalized abscissas are stretched to each spline interval.

! Each polynomial piece is integrated and accumulated.

qxi = alpha_*qx+beta_

quad(i) = sum(qw*spline_values(0, qxi, break_points,&

coeff))*alpha_&

+ quad(i-1)

end do

! Normalize the coefficients and partial integrals so that the

! total integral has the value one.

coeff=coeff/quad(ndata-1); quad=quad/quad(ndata-1)

rn=rand(rn)

x=zero; niterat=0

solve_equation: do

! Find the intervals where the x values are located.

LEFT_OF=NDEGREE; I=NDEGREE

I=I+1; if(I >= LAST) EXIT

WHERE(x >= BKPT(I))LEFT_OF = LEFT_OF+1

end do

! Use Newton's method to solve the nonlinear equation:

! accumulated_distribution_function - random_number = 0.

alpha_x = (x-bkpt(LEFT_OF))*half

beta_x = (x+bkpt(LEFT_OF))*half

FN=QUAD(LEFT_OF-NORD)-RN

DO I=1,NQUAD

FN=FN+QW(I)*spline_values(0, alpha_x*QX(I)+beta_x,&

break_points, coeff)*alpha_x

END DO

! This is the Newton method update step:

x=x-fn/spline_values(0, x, break_points, coeff)

niterat=niterat+1

! Constrain the values so they fall back into the interval.

! Newton's method may give approximates outside the interval.

where(x <= -one .or. x >= one) x=zero

if(norm(fn,1) <= sqrt(epsilon(one))*norm(x,1))&

exit solve_equation

end do solve_equation

! Check that Newton's method converges.

if (niterat <= limit) then

write (*,*) 'Example 3 for SPLINE_FITTING is correct.'

end if

end

Output

Example 3 for SPLINE_FITTING is correct.

Example 4: Represent a Periodic Curve

The curve tracing the edge of a rectangular box, traversed in a counter-clockwise direction, is parameterized with a spline representation for each coordinate function, (x(t), y(t)). The functions are constrained to be periodic at the ends of the parameter interval. Since the perimeter arcs are piece-wise linear functions, the degree of the splines is the value one. Some breakpoints are chosen so they correspond to corners of the box, where the derivatives of the coordinate functions are discontinuous. The value of this representation is that for each t the splines representing (x(t), y(t)) are points on the perimeter of the box. This “eases” the complexity of evaluating the edge of the box. This example illustrates a method for representing the edge of a domain in two dimensions, bounded by a periodic curve.

use spline_fitting_int

use norm_int

implicit none

! This is Example 4 for SPLINE_FITTING. Use piecewise-linear

! splines to represent the perimeter of a rectangular box.

integer i, j

integer, parameter :: nbkpt=9, nord=2, ndegree=nord-1, &

ncoeff=nbkpt-nord, ndata=7, ngrid=100, &

nvalues=(ndata-1)*ngrid

real(kind(1e0)), parameter :: zero=0e0, one=1e0

real(kind(1e0)), parameter :: delta_t=one, delta_b=one, delta_v=0.01

real(kind(1e0)) delta_x, delta_y

real(kind(1e0)), dimension(ndata) :: sddata=one, &

! These are redundant coordinates on the edge of the box.

xdata=(/0.0, 1.0, 2.0, 2.0, 1.0, 0.0, 0.0/), &

ydata=(/0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0/)

real(kind(1e0)) tdata(ndata), xspline_data(3, ndata), &

yspline_data(3, ndata), tvalues(nvalues), &

xvalues(nvalues), yvalues(nvalues), xcoeff(ncoeff), &

ycoeff(ncoeff), xcheck(nvalues), ycheck(nvalues), diffs

real(kind(1e0)), target :: bkpt(nbkpt)

real(kind(1e0)), pointer :: pointer_bkpt(:)

type (s_spline_knots) break_points

type (s_spline_constraints) constraints(1)

tdata = (/((i-1)*delta_t, i=1,ndata)/)

xspline_data(1,:)=tdata; yspline_data(1,:)=tdata

xspline_data(2,:)=xdata; yspline_data(2,:)=ydata

xspline_data(3,:)=sddata; yspline_data(3,:)=sddata

bkpt(nord:nbkpt-ndegree)=(/((i-nord)*delta_b, &

i=nord, nbkpt-ndegree)/)

! Collapse the outside knots.

bkpt(1:ndegree)=bkpt(nord)

bkpt(nbkpt-ndegree+1:nbkpt)=bkpt(nbkpt-ndegree)

! Assign the degree of the polynomial and the knots.

pointer_bkpt => bkpt

break_points=s_spline_knots(ndegree, pointer_bkpt)

! Make the two parametric curves also periodic.

constraints(1)=spline_constraints &

(derivative=0, point=bkpt(nord), type='.=.', &

value=bkpt(nbkpt-ndegree))

xcoeff = spline_fitting(data=xspline_data, knots=break_points, &

constraints=constraints)

ycoeff = spline_fitting(data=yspline_data, knots=break_points, &

constraints=constraints)

! Use the splines to compute the coordinates of points along the perimeter.

! Compare them with the coordinates of the edge points.

tvalues= (/((i-1)*delta_v, i=1,nvalues)/)

xvalues=spline_values(0, tvalues, break_points, xcoeff)

yvalues=spline_values(0, tvalues, break_points, ycoeff)

do i=1, nvalues

j=(i-1)/ngrid+1

delta_x=(xdata(j+1)-xdata(j))/ngrid

delta_y=(ydata(j+1)-ydata(j))/ngrid

xcheck(i)=xdata(j)+mod(i+ngrid-1,ngrid)*delta_x

ycheck(i)=ydata(j)+mod(i+ngrid-1,ngrid)*delta_y

end do

diffs=norm(xvalues-xcheck,1)/norm(xcheck,1)+&

norm(yvalues-ycheck,1)/norm(ycheck,1)

if (diffs <= sqrt(epsilon(one))) then

write(*,*) 'Example 4 for SPLINE_FITTING is correct.'

end if

end

Output

Example 4 for SPLINE_FITTING is correct.