This tutorial explains how to use the fitpack_curve type for univariate spline smoothing and interpolation of the form \( y = s(x) \). It covers the mathematical background, basic and advanced usage, and practical guidance on choosing a smoothing parameter.

Use fitpack_curve when you have a single-valued relationship between an independent variable \( x \) and a dependent variable \( y \) — that is, each \( x_i \) maps to exactly one \( y_i \). The data must be sorted in strictly increasing order of \( x \).

Mathematical Background

Given \( m \) data points \( (x_i, y_i) \) with positive weights \( w_i \), FITPACK constructs a B-spline \( s(x) \) of degree \( k \) (default \( k = 3 \), cubic) by solving the constrained optimization problem

\[ \min \int_{x_1}^{x_m} \left[ s''(x) \right]^2 dx \quad \text{subject to} \quad \sum_{i=1}^{m} w_i^2 \left( y_i - s(x_i) \right)^2 \leq S \]

The smoothing parameter \( S \geq 0 \) controls the trade-off between roughness and fidelity to the data:

** \( S = 0 \)** — Interpolation. The spline passes through every data point exactly.
**Small \( S \)** — Close to the data but may oscillate between points.
**Large \( S \)** — Very smooth, but may miss genuine features.

Choosing S

When the data satisfy \( y_i = g(x_i) + \varepsilon_i \) with independent errors of variance \( \sigma^2 \) and unit weights, the expected residual sum of squares for the true function is \( m \sigma^2 \). A natural starting point is therefore

\[ S \approx m \, \sigma^2 \]

If no noise estimate is available, \( S = m \) (assuming unit weights) is a reasonable first guess. The algorithm selects knots automatically: fewer knots for larger \( S \), more knots for smaller \( S \). After fitting, inspect curve\mse() and curve\knots to judge the result.

Basic Example

The following program fits a smoothing spline to noisy sine data, evaluates it, and prints the result.

program basic_curve
    use fitpack, only: fitpack_curve, fp_real, fp_flag, &
                       fitpack_success, fitpack_message, &
                       zero, one, half, pi
    implicit none
 
    integer, parameter :: m = 50
    real(FP_REAL), parameter :: two = 2.0_fp_real
 
    real(FP_REAL) :: x(m), y(m), dx
    type(fitpack_curve) :: curve
    integer(FP_FLAG) :: ierr
    real(FP_REAL) :: y_at_pi
    integer :: i
 
    ! Generate data on [0, 2*pi]
    dx = two * pi / (m - 1)
    do i = 1, m
        x(i) = (i - 1) * dx
        y(i) = sin(x(i)) + 0.05_fp_real * sin(17.0_fp_real * x(i))
    end do
 
    ! Fit a smoothing spline (S = m)
    ierr = curve%new_fit(x, y, smoothing=real(m, fp_real))
    if (.not. fitpack_success(ierr)) then
        print *, 'Fit failed: ', fitpack_message(ierr)
        error stop
    end if
 
    ! Evaluate at x = pi
    y_at_pi = curve%eval(pi, ierr)
    print '(a,f10.6)', 's(pi)  = ', y_at_pi
    print '(a,i0)',    'knots  = ', curve%knots
    print '(a,es10.3)', 'mse    = ', curve%mse()
 
    call curve%destroy()
end program basic_curve

Step-by-step

Import § The use fitpack statement exposes the curve type, kind parameters (FP_REAL, FP_FLAG), error-handling helpers (FITPACK_SUCCESS, FITPACK_MESSAGE), and named constants (zero, one, half, pi). Note that two is not a public export — define it locally when needed.
Prepare data § The abscissae x must be strictly increasing. The corresponding values y may contain noise.
Fit § curve\new_fit(x, y, smoothing=S) constructs the spline. It returns an integer(FP_FLAG) error code. Check it with FITPACK_SUCCESS(ierr).
Evaluate § curve\eval(x, ierr) accepts a scalar or an array and returns the spline value(s).
Clean up § curve\destroy() releases all internal allocations.

Advanced Features

Weights

Supply per-point weights via the w argument to emphasize or de-emphasize individual observations. Points with larger weights are fitted more closely.

real(FP_REAL) :: w(m)
 
! Inverse-variance weighting
w = one / sigma
ierr = curve%new_fit(x, y, w=w, smoothing=real(m, fp_real))

When weights are proportional to \( 1 / \sigma_i \), the smoothing constraint becomes \( \sum (y_i - s(x_i))^2 / \sigma_i^2 \leq S \), and \( S \approx m \) remains a natural choice.

Derivatives

Evaluate derivatives of the fitted spline up to order \( k \) (the spline degree, default 3):

real(FP_REAL) :: dydx, d2ydx2
 
dydx   = curve%dfdx(pi / 4, order=1, ierr=ierr)   ! first derivative
d2ydx2 = curve%dfdx(pi / 4, order=2, ierr=ierr)   ! second derivative

The order parameter selects which derivative: 1 for \( s'(x) \), 2 for \( s''(x) \), 3 for \( s'''(x) \).

Integration

Compute the definite integral of the spline over an arbitrary sub-interval \( [a, b] \):

real(FP_REAL) :: area
 
area = curve%integral(zero, pi)

This evaluates \( \int_a^b s(x) \, dx \) analytically from the B-spline coefficients, so the result is exact for the fitted spline (no numerical quadrature error).

Zero Finding

Find all zeros (roots) of a cubic spline within the data range:

real(FP_REAL), allocatable :: roots(:)
 
roots = curve%zeros(ierr)
print '(a,i0)', 'Number of zeros: ', size(roots)

The curve\zeros(ierr) function returns an allocatable array containing every root of \( s(x) = 0 \). This operation requires a cubic spline ( \( k = 3 \)).

Interpolation

For exact interpolation — a spline that passes through every data point — set smoothing=zero:

ierr = curve%new_fit(x, y, smoothing=zero)

Alternatively, call the dedicated interpolation method:

ierr = curve%interpolate()

Both produce the smoothest spline (minimum roughness integral) that satisfies \( s(x_i) = y_i \) for all \( i \).

Smoothing Sweep

When the appropriate smoothing level is unknown, try a sequence of \( S \) values on a logarithmic scale and inspect the resulting number of knots and residual:

integer :: j
 
print '(a)', '     S       knots   residual'
do j = 1, 5
    ierr = curve%new_fit(x, y, smoothing=10.0_fp_real**(3 - j))
    if (fitpack_success(ierr)) then
        print '(es10.1,i8,es12.3)', &
            curve%smoothing, curve%knots, curve%mse()
    end if
end do

A good fit typically shows a "knee" where adding more knots (decreasing \( S \)) yields diminishing improvement in the residual. Choose the \( S \) value at or just past this transition point.

API Summary

Method	Description
`curve\new_fit(x, y, w, smoothing, order)`	Fit a smoothing spline; returns `integer(FP_FLAG)`
`curve\eval(x, ierr)`	Evaluate \( s(x) \) at a scalar or array of points
`curve\dfdx(x, order, ierr)`	Evaluate derivative of order 1, 2, or 3
`curve\integral(from, to)`	Definite integral \( \int_a^b s(x) \, dx \)
`curve\zeros(ierr)`	All zeros of \( s(x) \) (cubic splines only)
`curve\interpolate()`	Exact interpolation ( \( S = 0 \))
`curve\knots`	Number of knots in the current fit
`curve\mse()`	Weighted residual sum of squares \( f_p \)
`curve\smoothing`	Current smoothing parameter \( S \)
`curve\destroy()`	Release all allocated memory

Complete Example

A fully compilable program demonstrating smoothing, interpolation, derivatives, integration, root finding, and a smoothing sweep is provided in examples/example_curve.f90. Build and run it with:

fpm run --example example_curve

Error Handling

All fitting and evaluation routines return an integer(FP_FLAG) error code. Use FITPACK_SUCCESS(ierr) to test for success and FITPACK_MESSAGE(ierr) to obtain a human-readable description:

if (.not. fitpack_success(ierr)) then
    print *, 'Error: ', fitpack_message(ierr)
end if

Common return values:

Flag	Meaning
`FITPACK_OK` (0)	Success
`FITPACK_INTERPOLATING_OK` ( \(-1\))	Interpolating spline returned
`FITPACK_S_TOO_SMALL` (1)	Smoothing parameter too small for the knot budget
`FITPACK_MAXIT` (2)	Maximum iterations reached
`FITPACK_INPUT_ERROR` (10)	Invalid input data

See also: Curve Fitting Theory; B-Spline Foundations; Periodic Curve Fitting