B-Spline and Piecewise Polynomial Examples¶
Example Code¶
Representative script:
"""IMSL BS2/BS3 and piecewise-polynomial family examples with plots.
Covers:
- BS2IN, BS2VL, BS2DR, BS2GD, BS2IG
- BS3IN, BS3VL, BS3DR, BS3GD, BS3IG
- BSCPP, PPVAL, PPDER, PP1GD, PPITG
"""
from __future__ import annotations
from pathlib import Path
from typing import Any, Dict
import sys
import matplotlib.pyplot as plt
import numpy as np
REPO_ROOT = Path(__file__).resolve().parent.parent
if str(REPO_ROOT) not in sys.path:
sys.path.insert(0, str(REPO_ROOT))
from interpolation import (
bspline_2d_derivative,
bspline_2d_evaluate,
bspline_2d_evaluate_grid,
bspline_2d_integrate,
bspline_2d_interpolate,
bspline_3d_derivative,
bspline_3d_evaluate,
bspline_3d_evaluate_grid,
bspline_3d_integrate,
bspline_3d_interpolate,
bspline_evaluate_grid,
bspline_interpolate,
bspline_to_piecewise_poly,
piecewise_poly_derivative,
piecewise_poly_evaluate,
piecewise_poly_evaluate_grid,
piecewise_poly_integrate,
)
def run_demo_imsl_bs23_pp_family() -> Dict[str, Any]:
"""Run website-style demos for BS2/BS3 and PP families."""
out: Dict[str, Any] = {}
out_dir = Path("test_output")
out_dir.mkdir(parents=True, exist_ok=True)
# BS2 family ---------------------------------------------------------------
kx2, ky2 = 5, 2
x2 = np.array([(i - 11) / 10.0 for i in range(1, 22)], dtype=float)
y2 = np.array([(i - 1) / 5.0 for i in range(1, 7)], dtype=float)
xx2, yy2 = np.meshgrid(x2, y2, indexing="ij")
f2 = xx2**3 + xx2 * yy2
rep2 = bspline_2d_interpolate(x2, y2, f2, order_x=kx2, order_y=ky2)
xv = np.linspace(0.0, 1.0, 4)
yv = np.linspace(0.0, 1.0, 4)
s2_grid = bspline_2d_evaluate_grid(rep2, xv, yv)
f2_grid = np.array([[x**3 + x * y for y in yv] for x in xv], dtype=float)
out["BS2IN"] = {"grid_values": s2_grid, "grid_error": f2_grid - s2_grid}
out["BS2VL"] = float(bspline_2d_evaluate(rep2, 1.0, 1.0))
# Derivative example for BS2DR/BS2GD (stable mixed derivative target = 1)
xv_d = np.linspace(0.2, 0.8, 4)
yv_d = np.linspace(0.2, 0.8, 4)
s21 = np.array(
[[bspline_2d_derivative(rep2, float(x), float(y), dx_order=1, dy_order=1) for y in yv_d] for x in xv_d],
dtype=float,
)
f21 = np.ones((xv_d.size, yv_d.size), dtype=float)
out["BS2DR"] = float(bspline_2d_derivative(rep2, 1.0, 1.0, dx_order=1, dy_order=1))
out["BS2GD"] = {"grid_values": s21, "grid_error": f21 - s21}
# BS2IG example
bs2ig_val = float(bspline_2d_integrate(rep2, 0.0, 1.0, 0.5, 1.0))
bs2ig_exact = 0.3125
out["BS2IG"] = {"computed": bs2ig_val, "exact": bs2ig_exact, "error": bs2ig_exact - bs2ig_val}
fig2, ax2 = plt.subplots(figsize=(8, 5))
c2 = ax2.contourf(xv, yv, s2_grid.T, levels=12, cmap="viridis")
fig2.colorbar(c2, ax=ax2, label="S(x,y)")
ax2.set_title("IMSL BS2IN/BS2VL")
ax2.set_xlabel("x")
ax2.set_ylabel("y")
fig2.tight_layout()
path2 = out_dir / "demo_imsl_bs2_family.svg"
fig2.savefig(path2, format="svg")
plt.close(fig2)
out["bs2_plot_path"] = str(path2)
# BS3 family ---------------------------------------------------------------
kx3, ky3, kz3 = 5, 2, 3
x3 = np.array([(i - 11) / 10.0 for i in range(1, 22)], dtype=float)
y3 = np.array([(i - 1) / 5.0 for i in range(1, 7)], dtype=float)
z3 = np.array([(i - 1) / 7.0 for i in range(1, 9)], dtype=float)
xx3, yy3, zz3 = np.meshgrid(x3, y3, z3, indexing="ij")
f3 = xx3**3 + xx3 * yy3 * zz3
rep3 = bspline_3d_interpolate(x3, y3, z3, f3, order_x=kx3, order_y=ky3, order_z=kz3)
xg = np.array([-1.0, -1.0 / 3.0, 1.0 / 3.0, 1.0], dtype=float)
yg = np.array([0.0, 1.0 / 3.0, 2.0 / 3.0, 1.0], dtype=float)
zg = np.array([0.0, 1.0], dtype=float)
s3_grid = bspline_3d_evaluate_grid(rep3, xg, yg, zg)
f3_grid = np.array([[[x**3 + x * y * z for z in zg] for y in yg] for x in xg], dtype=float)
out["BS3IN"] = {"grid_values": s3_grid, "grid_error": f3_grid - s3_grid}
out["BS3VL"] = float(bspline_3d_evaluate(rep3, 1.0, 1.0, 1.0))
# BS3DR/BS3GD example (stable mixed derivative target = 1)
xg_d = np.array([-0.8, -0.4, 0.0, 0.4], dtype=float)
yg_d = np.array([0.2, 0.4, 0.6, 0.8], dtype=float)
zg_d = np.array([0.2, 0.4], dtype=float)
s201_grid = np.array(
[
[
[bspline_3d_derivative(rep3, float(x), float(y), float(z), dx_order=1, dy_order=1, dz_order=1) for z in zg_d]
for y in yg_d
]
for x in xg_d
],
dtype=float,
)
f201_grid = np.ones((xg_d.size, yg_d.size, zg_d.size), dtype=float)
out["BS3DR"] = float(bspline_3d_derivative(rep3, 1.0, 1.0, 1.0, dx_order=1, dy_order=1, dz_order=1))
out["BS3GD"] = {"grid_values": s201_grid, "grid_error": f201_grid - s201_grid}
# Use an asymmetric quadrature grid that preserves accuracy while
# avoiding CI timeout in coverage runs.
bs3ig_val = float(
bspline_3d_integrate(
rep3,
0.0,
1.0,
0.5,
1.0,
0.0,
0.5,
nx=28,
ny=20,
nz=16,
)
)
bs3ig_exact = 11.0 / 128.0
out["BS3IG"] = {"computed": bs3ig_val, "exact": bs3ig_exact, "error": bs3ig_exact - bs3ig_val}
fig3, ax3 = plt.subplots(figsize=(8, 5))
mid_slice = s3_grid[:, :, 1]
c3 = ax3.contourf(xg, yg, mid_slice.T, levels=12, cmap="plasma")
fig3.colorbar(c3, ax=ax3, label="S(x,y,z=1)")
ax3.set_title("IMSL BS3IN/BS3VL")
ax3.set_xlabel("x")
ax3.set_ylabel("y")
fig3.tight_layout()
path3 = out_dir / "demo_imsl_bs3_family.svg"
fig3.savefig(path3, format="svg")
plt.close(fig3)
out["bs3_plot_path"] = str(path3)
# BSCPP / PPVAL / PPDER / PP1GD / PPITG --------------------------------
kpp = 4
npp = 20
xpp = np.linspace(0.0, 1.0, npp)
fpp = xpp * np.exp(xpp)
rep_pp = bspline_interpolate(xpp, fpp, order=kpp)
pp = bspline_to_piecewise_poly(rep_pp)
xvpp = np.linspace(0.0, 1.0, npp)
s_pp = np.asarray([piecewise_poly_evaluate(pp, float(x)) for x in xvpp], dtype=float)
ds_pp = np.asarray([piecewise_poly_derivative(pp, float(x), k=1) for x in xvpp], dtype=float)
exact_pp = xvpp * np.exp(xvpp)
exact_dpp = (xvpp + 1.0) * np.exp(xvpp)
out["BSCPP"] = {"n_intervals": int(pp.coefficients.shape[0]), "order": int(pp.order)}
out["PPVAL"] = {"values": s_pp, "error": exact_pp - s_pp}
out["PPDER"] = {"values": ds_pp, "error": exact_dpp - ds_pp}
out["PP1GD"] = {
"values": piecewise_poly_evaluate_grid(pp, xvpp),
"dvalues": np.asarray([piecewise_poly_derivative(pp, float(x), k=1) for x in xvpp], dtype=float),
}
# PPITG page example uses x^2 reproduction and two intervals.
xq = np.linspace(0.0, 2.0, 10)
yq = xq**2
rep_q = bspline_interpolate(xq, yq, order=3)
pp_q = bspline_to_piecewise_poly(rep_q)
ppitg_1 = piecewise_poly_integrate(pp_q, 0.0, 0.5)
ppitg_2 = piecewise_poly_integrate(pp_q, 0.0, 2.0)
out["PPITG"] = {
"computed": np.array([ppitg_1, ppitg_2], dtype=float),
"exact": np.array([1.0 / 24.0, 8.0 / 3.0], dtype=float),
}
figp, axp = plt.subplots(figsize=(8, 5))
axp.plot(xvpp, exact_pp, "k--", linewidth=1.2, label="x*exp(x)")
axp.plot(xvpp, s_pp, color="#1f77b4", linewidth=2.0, label="PPVAL")
axp.plot(xvpp, ds_pp, color="#ff7f0e", linewidth=1.6, label="PPDER")
axp.set_title("IMSL BSCPP/PPVAL/PPDER")
axp.set_xlabel("x")
axp.grid(True, alpha=0.3)
axp.legend(loc="best")
figp.tight_layout()
pathp = out_dir / "demo_imsl_pp_family.svg"
figp.savefig(pathp, format="svg")
plt.close(figp)
out["pp_plot_path"] = str(pathp)
return out
if __name__ == "__main__":
run_demo_imsl_bs23_pp_family()
Input (Console)¶
Run the B-spline and piecewise scripts from the package root:
python examples/example_imsl_bs23_pp_family.py
python examples/example_imsl_bsint_family.py
python examples/example_imsl_bsitg.py
python examples/example_imsl_bsnak_bsopk.py
Plot Output¶
Generated SVG plots:
Note: The y-axis is logarithmic so BSOPK remains visible when BSNAK spikes.¶
Output Console¶
Summary console output:
script input setup representative numeric outputs
example_imsl_bs23_pp_family.py 2D/3D tensor grids plus PP checks BS2 sample = 2.000000000000001; BS3 sample = 2; PP integrals = [0.04166666666666651, 2.666666666666669]
example_imsl_bsint_family.py 9 eval points for interp/deriv/grid BSINT head = [-1.669845252164844e-15, 0.2960102286930971, 0.4999999999999994]; BSDER = [1.049569367139735, 0.9262096826168076, 0.802849381531523]; max err = 0.05754316190017666
example_imsl_bsitg.py single integral over [-1,1] computed integral = 0.250000957410549; exact = 0.25; abs error = 9.574105490073315e-07
example_imsl_bsnak_bsopk.py spline orders 2 through 7 BSNAK errors = [0.05735177358221177, 5523.630105296993, 3.905851290599733, 28970035.06660343, 494385.4071861197, 12873.49377518255]; BSOPK errors = [0.01143705068921863, 0.003755372172662086, 0.07137867036154814, 1.91726455143818, 2.571673847046684, 49.40658062286245]
Detailed integration result tables:
BSINT family full output
x bsint_s bsint_error bsder_s bsder_ds bsder_err_s bsder_err_ds bs1gd_s bs1gd_ds
0 -1.669845252164844e-15 1.669845252164844e-15 0.4472135954999575 1.049569367139735 3.885780586188048e-16 0.06846462161015943 0.4472135954999575 1.049569367139735
0.125 0.2960102286930971 0.05754316190017666 0.5460025788245315 0.9262096826168076 0.00171997868063456 -0.01333875344153068 0.5460025788245315 0.9262096826168076
0.25 0.4999999999999994 5.551115123125783e-16 0.6324555320336757 0.802849381531523 2.220446049250313e-16 -0.0122799664894282 0.6324555320336757 0.802849381531523
0.375 0.6166686020914135 -0.004296166395619028 0.7069626626983782 0.7029015345882961 0.0001441184881694202 0.004205246598251389 0.7069626626983782 0.7029015345882961
0.5 0.7071067811865479 -3.33066907387547e-16 0.7745966692414832 0.6497785962067404 2.220446049250313e-16 -0.004281371838837589 0.7745966692414832 0.6497785962067404
0.625 0.7901116899682359 0.0004577250738589944 0.8369183819469426 0.5966556578806959 -0.0002583554128670063 0.0009586467865009274 0.8369183819469426 0.5966556578806959
0.75 0.8660254037844388 -2.220446049250313e-16 0.8944271909999161 0.558514425041956 -2.220446049250313e-16 0.0005025693329914205 0.8944271909999161 0.558514425041956
0.875 0.9359596675160309 -0.0005453208225455475 0.9487461145021723 0.527864044919113 -6.281645165850946e-05 -0.0008177682243830642 0.9487461145021723 0.527864044919113
1 1.000000000000003 -2.886579864025407e-15 1.000000000000002 0.4972138181180696 -1.554312234475219e-15 0.002786181881930361 1.000000000000002 0.4972138181180696
BSNAK vs BSOPK full order/error output
order bsnak_max_error bsopk_max_error
3 0.05735177358221177 0.01143705068921863
4 5523.630105296993 0.003755372172662086
5 3.905851290599733 0.07137867036154814
6 28970035.06660343 1.91726455143818
7 494385.4071861197 2.571673847046684
8 12873.49377518255 49.40658062286245