This tutorial covers shape-preserving spline fitting with fitpack_convex_curve, which enforces local convexity or concavity constraints on a cubic spline. Use it when the data represents a physical quantity that must be convex or concave — e.g. moisture content vs depth, dose-response curves, or thermodynamic isotherms — and an unconstrained spline would introduce spurious wiggles.

Mathematical Background

A cubic B-spline \( s(x) = \sum_j c_j B_j(x) \) is locally convex ( \( s''(x) \geq 0 \)) when the coefficient second differences satisfy

\[ c_{j-1} - 2\,c_j + c_{j+1} \geq 0 \]

and locally concave ( \( s''(x) \leq 0 \)) when

\[ c_{j-1} - 2\,c_j + c_{j+1} \leq 0 \]

FITPACK expresses these requirements through per-data-point flags \( v_i \in \{-1, 0, +1\} \) with the constraint \( s''(x_i)\,v_i \leq 0 \):

Flag	Constraint	Shape
\( v_i = +1 \)	\( s''(x_i) \leq 0 \)	Concave (curves downward)
\( v_i = -1 \)	\( s''(x_i) \geq 0 \)	Convex (curves upward)
\( v_i = 0 \)	none	Unconstrained

The fitting always uses cubic splines ( \( k = 3 \)). The core solver (concon) places knots automatically to satisfy the smoothing constraint while respecting the requested convexity at every data point.

See also: Dierckx, Ch. 8, §8.3–8.4 (pp. 173–196)

Basic Example

Consider moisture-content measurements as a function of depth. The data rises quickly then levels off — a concave profile.

use fitpack, only: fitpack_convex_curve, fp_real, fp_flag, &
                   fitpack_success, fitpack_message, one
implicit none
 
integer, parameter :: m = 16
type(fitpack_convex_curve) :: curve
integer(FP_FLAG) :: ierr
real(FP_REAL)    :: x(m), y(m), v(m)
 
! Moisture content vs depth (16 points, concave trend)
x = [0.1_fp_real,  0.3_fp_real,  0.5_fp_real,  0.7_fp_real,  0.9_fp_real,  &
     1.25_fp_real, 1.75_fp_real, 2.25_fp_real, 2.75_fp_real, 3.5_fp_real,  &
     4.5_fp_real,  5.5_fp_real,  6.5_fp_real,  7.5_fp_real,  8.5_fp_real,  9.5_fp_real]
y = [0.124_fp_real, 0.234_fp_real, 0.256_fp_real, 0.277_fp_real, 0.278_fp_real, &
     0.291_fp_real, 0.308_fp_real, 0.311_fp_real, 0.315_fp_real, 0.322_fp_real, &
     0.317_fp_real, 0.326_fp_real, 0.323_fp_real, 0.321_fp_real, 0.322_fp_real, &
     0.328_fp_real]

Step 1 — Load Data and Set Convexity Flags

! Load data points
call curve%new_points(x, y)
 
! Concave everywhere: v = +1 means s''(x_i) <= 0
v = one
ierr = curve%set_convexity(v)
if (.not. fitpack_success(ierr)) error stop fitpack_message(ierr)

Step 2 — Fit with Constraints

Call curve\fit with a smoothing parameter. The solver selects knots automatically while enforcing \( s''(x_i) \leq 0 \) at every data point:

ierr = curve%fit(smoothing=0.04_fp_real)

if (.not. fitpack_success(ierr)) error stop fitpack_message(ierr)

Step 3 — Evaluate and Inspect

real(FP_REAL) :: y_eval, d2y
integer :: i
 
do i = 1, m
    y_eval = curve%eval(x(i), ierr)
    d2y    = curve%dfdx(x(i), order=2, ierr=ierr)
    print '(f8.3, f10.4, es12.3)', x(i), y_eval, d2y   ! d2y <= 0
end do
 
call curve%destroy()

Comparing Constrained vs Unconstrained

Without constraints, a standard fitpack_curve may oscillate through noisy data, creating regions where the second derivative changes sign. The constrained fit eliminates these wiggles:

use fitpack, only: fitpack_curve, fitpack_convex_curve, fp_real, fp_flag, one
 
type(fitpack_curve)        :: free_curve
type(fitpack_convex_curve) :: constrained
integer(FP_FLAG) :: ierr
real(FP_REAL)    :: v(m)
 
free_curve = fitpack_curve(x, y, ierr=ierr)          ! unconstrained
 
call constrained%new_points(x, y)                     ! constrained concave
v = one
ierr = constrained%set_convexity(v)
ierr = constrained%fit(smoothing=0.04_fp_real)
 
! Second derivative at an interior point
print *, free_curve%dfdx(3.0_fp_real, order=2, ierr=ierr)    ! may be > 0
print *, constrained%dfdx(3.0_fp_real, order=2, ierr=ierr)   ! guaranteed <= 0
 
call free_curve%destroy()
call constrained%destroy()

The constrained spline guarantees \( s''(x) \leq 0 \) throughout.

Mixed Constraints

The constraint vector need not be uniform. Assign different flags per point for data that changes curvature sign:

v(1:6)  = -one   ! convex region: s''(x_i) >= 0
v(7:m)  =  one   ! concave region: s''(x_i) <= 0
ierr = curve%set_convexity(v)

Setting \( v_i = 0 \) leaves individual points unconstrained, which is useful at transition regions where the curvature sign is unknown.

Constructor Shortcut

A one-step constructor combines new_points, set_convexity, and fit:

v = one

curve = fitpack_convex_curve(x, y, v, smoothing=0.04_fp_real, ierr=ierr)

API Summary

fitpack_convex_curve extends fitpack_curve with convexity constraints. The spline degree is always cubic ( \( k = 3 \)).

Method	Description
`curve\new_points(x, y, w)`	Load data points and optional weights
`curve\set_convexity(v)`	Set per-point convexity flags ( \( +1 \), \( -1 \), or \( 0 \))
`curve\fit(smoothing)`	Fit with automatic knots (calls `concon`)
`curve\least_squares()`	Least-squares fit with current knots (calls `cocosp`)
`curve\eval(x)`	Evaluate the spline
`curve\dfdx(x, order)`	Evaluate derivatives
`curve\integral(a, b)`	Definite integral over \( [a, b] \)
`curve\destroy()`	Release allocated memory

Complete Example

See examples/example_convex_curve.f90 (examples/example_convex_curve.f90) for a full working program that loads moisture-content data with endpoint-emphasized weights, fits at multiple smoothing levels, and prints a comparison table of knot count, residual, and maximum second derivative.

See also: fitpack_convex_curves, Dierckx Ch. 8 (pp. 173–196)