"""Numerical first and second derivatives from tabulated data. The program uses local four-point Lagrange interpolation. The data points may be unequally spaced, and the stencil is shifted automatically near either endpoint. Based on deriv_Lg.f90 by Alexander Godunov (October 2009). Revised for the companion website (September 2026). """ from math import exp from typing import Sequence STENCIL_SIZE = 4 def main() -> None: """Generate sample data and calculate first and second derivatives.""" number_of_data = 31 x_test = [0.00, 0.20, 0.70, 1.25, 1.95, 2.75, 3.00] x_data = [0.0] * number_of_data y_data = [0.0] * number_of_data #------------------------------------------------------------------- # Step 1: Generate sample tabulated data. # # The example uses f(x)=exp(x) on a mildly nonuniform grid. Both # exact derivatives equal exp(x), so the results are easy to check. # # For measured or previously calculated data, replace this loop # with assignments to x_data and y_data or with file-input code. #------------------------------------------------------------------- for i in range(number_of_data): fraction = i / (number_of_data - 1) x_data[i] = 3.0 * fraction**1.15 y_data[i] = exp(x_data[i]) #------------------------------------------------------------------- # Step 2: Calculate the derivatives at the requested points. #------------------------------------------------------------------- print("Numerical derivatives from tabulated data") print("Local four-point Lagrange interpolation\n") print( " x exact f' numerical f'" " exact f'' numerical f''" ) for x_eval in x_test: first_derivative, status_first = derivative_from_data( x_eval, x_data, y_data, 1 ) second_derivative, status_second = derivative_from_data( x_eval, x_data, y_data, 2 ) if status_first != 0 or status_second != 0: print( f"Derivative could not be evaluated at x = {x_eval:.3f}; " f"status = {status_first} {status_second}" ) continue exact = exp(x_eval) print( f"{x_eval:9.3f}" f"{exact:18.8e}" f"{first_derivative:18.8e}" f"{exact:18.8e}" f"{second_derivative:18.8e}" ) def derivative_from_data( x_eval: float, x_data: Sequence[float], y_data: Sequence[float], derivative_order: int, ) -> tuple[float, int]: """Estimate a first or second derivative from tabulated data. A cubic Lagrange polynomial is constructed from four consecutive data points near x_eval. Differentiating this polynomial gives the requested derivative. Parameters ---------- x_eval Point where the derivative is required. x_data Strictly increasing data abscissas. y_data Tabulated function values. derivative_order 1 for the first derivative or 2 for the second derivative. Returns ------- derivative, status The numerical derivative and a status code: 0: success 1: x_data and y_data have different sizes 2: fewer than four data points were supplied 3: derivative_order is not 1 or 2 4: x_data is not strictly increasing 5: x_eval is outside the tabulated interval """ derivative = 0.0 number_of_points = len(x_data) # Check the input before selecting the interpolation stencil. if len(y_data) != number_of_points: return derivative, 1 if number_of_points < STENCIL_SIZE: return derivative, 2 if derivative_order not in (1, 2): return derivative, 3 for i in range(number_of_points - 1): if x_data[i + 1] <= x_data[i]: return derivative, 4 if x_eval < x_data[0] or x_eval > x_data[-1]: return derivative, 5 # Select four nearby points. At a boundary, the stencil is shifted # automatically so that all four points remain in the data range. first = select_stencil(x_eval, x_data) x_local = x_data[first : first + STENCIL_SIZE] y_local = y_data[first : first + STENCIL_SIZE] derivative = lagrange_derivative( x_eval, x_local, y_local, derivative_order ) return derivative, 0 def select_stencil(x_eval: float, x_data: Sequence[float]) -> int: """Locate the first point of the local four-point stencil. The binary search finds ``left`` such that x_data[left] <= x_eval < x_data[left+1]. The stencil is centered on this interval as closely as possible. """ number_of_points = len(x_data) left = 0 right = number_of_points - 1 while right > left + 1: middle = (left + right) // 2 if x_eval < x_data[middle]: right = middle else: left = middle first = left - (STENCIL_SIZE - 2) // 2 first = max(0, min(first, number_of_points - STENCIL_SIZE)) return first def lagrange_derivative( x_eval: float, x: Sequence[float], y: Sequence[float], derivative_order: int, ) -> float: """Differentiate the local Lagrange interpolating polynomial. For the local data (x_j, y_j), P(x) = sum_j y_j L_j(x), where L_j(x) = product_{k != j} (x-x_k)/(x_j-x_k). The loops evaluate the first or second derivative of each basis polynomial directly. No constant grid spacing is assumed. """ number_of_points = len(x) value = 0.0 for j in range(number_of_points): # Denominator of the j-th Lagrange basis polynomial. denominator = 1.0 for k in range(number_of_points): if k != j: denominator *= x[j] - x[k] if derivative_order == 1: # In each term of L'_j(x), one factor is differentiated. for k in range(number_of_points): if k == j: continue term = 1.0 for i in range(number_of_points): if i != j and i != k: term *= x_eval - x[i] value += y[j] * term / denominator else: # In each term of L''_j(x), two factors are differentiated. for k in range(number_of_points): if k == j: continue for ell in range(number_of_points): if ell == j or ell == k: continue term = 1.0 for i in range(number_of_points): if i != j and i != k and i != ell: term *= x_eval - x[i] value += y[j] * term / denominator return value if __name__ == "__main__": main()