5. Function Approximation, Interpolation, and Data Fitting

Introduction and Context

Engineering analysis rarely provides a function in an exact, convenient form. A temperature field may be known only at sensor locations, a constitutive law may be available as scattered experimental data, or a simulation may produce values at selected design points rather than an explicit formula. Computational methods provide ways to construct a useful mathematical representation from such incomplete information.

Function approximation is the broader task of representing a complicated or unknown function by a simpler one. Interpolation constructs a representation that passes through specified data points, whereas data fitting seeks a representation that captures an overall trend and may deliberately avoid passing through every measurement. These distinctions are important. Exact interpolation is appropriate when data are regarded as accurate, while regression is generally preferable when measurements contain noise.

This chapter develops polynomial and spline interpolation, least-squares regression, orthogonal-function approximations, response surfaces, and regularization. It also examines approximation error, conditioning, validation, overfitting, and the professional judgment required to select an appropriate model.

Learning Objectives

After studying this chapter, the reader should be able to:

1. Distinguish interpolation, approximation, regression, and extrapolation.

2. Construct and evaluate polynomial interpolants using Newton and Lagrange forms.

3. Explain numerical instability in high-degree polynomial interpolation and the role of node selection.

4. Formulate and interpret piecewise linear and cubic spline interpolants.

5. Derive and solve least-squares regression problems, including multivariable models.

6. Use orthogonal bases to improve numerical conditioning and interpret model coefficients.

7. Build response-surface models for engineering design and optimization.

8. Diagnose approximation error, bias, variance, and overfitting.

9. Apply regularization and validation techniques to obtain reliable predictive models.

10. Select methods responsibly in experimental, computational, and industrial settings.

5.1 Approximation as a Modelling Decision

Let f(x) denote an unknown or expensive function and let g(x;θ) be a computationally manageable approximation. The parameters θ may be determined from known function values, derivatives, measurements, or simulation results. The quality of an approximation depends not only on the algebraic form of g but also on the intended use. A model used to estimate a value between measured points has different requirements from one used to optimize a design or predict behavior outside the measured range.

Interpolation uses data pairs (x_i,y_i) and imposes

g(x_i)=y_i, i=0,1,...,n.

When the y_i values are exact, this constraint is desirable. When they contain random error, forcing the curve through each point can reproduce measurement noise. Regression instead minimizes a specified discrepancy between observations and model predictions. Approximation therefore involves a balance between fidelity, smoothness, computational cost, interpretability, and extrapolation risk.

The distinction between interpolation and extrapolation must be maintained. Interpolation estimates values within the range of the data; extrapolation estimates values beyond it. Even a highly accurate interpolant can produce unreliable extrapolations because the data impose little constraint outside the observed domain.

5.2 Polynomial Interpolation

Lagrange and Newton Forms

For n+1 distinct nodes, there is a unique polynomial p_n(x) of degree at most n satisfying p_n(x_i)=y_i. In Lagrange form,

p_n(x)=Σ(i=0 to n) y_i L_i(x),

where

L_i(x)=Π(j≠i) (x-x_j)/(x_i-x_j).

Each L_i is one at x_i and zero at all other nodes. This form is conceptually transparent and useful for deriving error formulas, but direct recomputation can be inefficient when data are added.

The Newton divided-difference form is

p_n(x)=f[x_0]+f[x_0,x_1](x-x_0)+...+f[x_0,...,x_n]Π(k=0 to n-1)(x-x_k).

Divided differences are computed recursively:

f[x_i,x_{i+1}]=(f[x_{i+1}]-f[x_i])/(x_{i+1}-x_i),

and

f[x_i,...,x_{i+k}]=[f[x_{i+1},...,x_{i+k}]-f[x_i,...,x_{i+k-1}]]/(x_{i+k}-x_i).

Newton form is particularly useful when nodes are introduced sequentially. Evaluation is normally performed using nested multiplication, reducing arithmetic operations and round-off effects.

Interpolation Error

If f has n+1 continuous derivatives on an interval containing x and the nodes, then

f(x)-p_n(x)=f^(n+1)(ξ)/(n+1)! Π(i=0 to n)(x-x_i)

for some ξ in the interval. The error depends on two factors: the smoothness of f and the distribution of nodes. Increasing polynomial degree does not automatically reduce error because the product term may become large and the derivatives may grow rapidly.

For equally spaced nodes, high-degree interpolation may exhibit strong oscillations near interval boundaries, known as the Runge phenomenon. Chebyshev nodes, clustered near the endpoints, reduce the maximum interpolation error and are often preferred for global polynomial interpolation. In practice, barycentric interpolation provides a numerically stable way to evaluate a polynomial in Lagrange form without explicitly forming its coefficients.

Limitations of Global Polynomials

A single high-degree polynomial is globally coupled: changing one data point can alter the approximation throughout the interval. This is undesirable for large data sets, local irregularities, or data arriving incrementally. Polynomial coefficients may also be poorly conditioned in the monomial basis 1,x,x²,..., especially when x has a large magnitude or the degree is high. Scaling x to a standard interval such as [-1,1] is a simple but important improvement.

5.3 Piecewise Interpolation and Splines

Piecewise methods avoid the excessive global flexibility of a high-degree polynomial. The interval is divided into subintervals, and a low-degree polynomial is fitted locally.

Piecewise linear interpolation joins neighboring data points by straight lines. On [x_i,x_(i+1)],

p(x)=y_i+(y_(i+1)-y_i)(x-x_i)/(x_(i+1)-x_i).

It is easy to implement and preserves monotonicity between points, but its first derivative is discontinuous at interior nodes. This can be unacceptable in mechanics, where slopes represent velocity, heat flux, stress gradients, or sensitivity.

A cubic spline uses a cubic polynomial on each interval,

s_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)²+d_i(x-x_i)³,

with continuity of the function, first derivative, and second derivative at interior nodes. A standard natural cubic spline imposes s''(x_0)=s''(x_n)=0. Other boundary conditions include clamped splines, in which endpoint derivatives are specified, and not-a-knot conditions, which effectively merge the first two and last two intervals.

The cubic spline coefficients are obtained from a tridiagonal linear system involving the second derivatives at the nodes. The tridiagonal structure permits efficient solution and makes splines attractive for large data sets. Cubic splines are smooth, local, and generally much less oscillatory than high-degree global polynomials.

A smoothing spline extends interpolation to noisy data by minimizing

Σ_i [y_i-s(x_i)]² + λ∫[s''(x)]²dx.

The first term measures data misfit and the second penalizes curvature. The parameter λ controls the trade-off: λ=0 gives an interpolating spline, while larger λ produces a smoother curve. This is a direct example of regularization.

5.4 Least-Squares Regression

Suppose observations y_i are approximated by a model linear in its parameters:

y_i≈Σ(j=1 to p) θ_j φ_j(x_i).

The basis functions φ_j may be 1,x,x², trigonometric functions, or physically motivated functions. In matrix form,

Aθ≈y,

where A_ij=φ_j(x_i). The ordinary least-squares estimate minimizes

J(θ)=||Aθ-y||²_2.

Setting the gradient to zero produces the normal equations,

A^T A θ=A^T y.

If A has full column rank, the solution is θ=(A^TA)^−1A^Ty in exact arithmetic. However, explicitly forming and inverting A^TA is numerically undesirable because it squares the condition number. QR factorization is generally preferred, while singular-value decomposition is especially valuable for rank-deficient or ill-conditioned problems.

For a straight-line model y≈β_0+β_1x, β_1 measures the fitted change in y per unit x and β_0 is the fitted intercept. A coefficient should not be interpreted causally merely because it is statistically or numerically significant; the model may omit important physical variables.

Weighted and Generalized Least Squares

Measurements often have unequal reliability. Weighted least squares minimizes

J(θ)=Σ_i w_i[y_i-g(x_i;θ)]²,

with larger weights assigned to more precise observations. In matrix notation, the equations are A^TWAθ=A^TWy. If errors are correlated, a covariance-based formulation is more appropriate, and generalized least squares accounts for that covariance structure.

Model Assessment

The residual at observation i is r_i=y_i-g(x_i;θ). Useful diagnostics include residual plots, root-mean-square error, mean absolute error, maximum error, and coefficient of determination R². No single metric is sufficient. A low training error can coexist with poor predictive performance, especially for flexible models.

Residuals should be examined for systematic patterns. Curvature suggests an inadequate functional form; changing variance suggests heteroscedasticity; clusters or trends against time may indicate dependence; and isolated extreme residuals may signal outliers, sensor faults, or an unmodelled regime.

5.5 Orthogonal Functions and Stable Approximation

An inner product defines orthogonality. For functions on [a,b], one common definition is

⟨u,v⟩=∫_a^b w(x)u(x)v(x)dx,

where w(x)>0 is a weight function. Functions φ_i and φ_j are orthogonal if ⟨φ_i,φ_j⟩=0 for i≠j. Expanding a function as

f(x)≈Σ(k=0 to n)c_k φ_k(x)

is advantageous because, for an orthogonal basis, the coefficients are obtained independently:

c_k=⟨f,φ_k⟩/⟨φ_k,φ_k⟩.

In numerical work, Chebyshev and Legendre polynomials are widely used. Chebyshev polynomials are particularly useful for minimizing approximation error on a finite interval, whereas Legendre polynomials are orthogonal with unit weight on [-1,1]. Fourier bases are natural for periodic phenomena.

Orthogonalization reduces parameter correlation and improves conditioning compared with the monomial basis. Nevertheless, numerical orthogonality may be lost if algorithms are implemented carelessly; modified Gram–Schmidt or QR-based methods are safer than naive Gram–Schmidt.

5.6 Response Surfaces for Engineering Design

A response surface is an approximate relationship between controllable design variables and a performance response. For two variables, a second-order model may be written as

y=β_0+β_1x_1+β_2x_2+β_{12}x_1x_2+β_{11}x_1²+β_{22}x_2²+ε.

The linear terms describe local trends, interaction terms indicate that the effect of one variable depends on another, and quadratic terms represent curvature. The coefficients are usually estimated by least squares using a designed set of experiments or simulations.

Response-surface methodology combines design of experiments, regression, visualization, and optimization. Factor ranges should be selected using engineering knowledge, and variables are often coded to [-1,1] so that coefficients have comparable numerical scales. Central composite and Box–Behnken designs are examples of structured designs for estimating second-order models.

A fitted surface can identify a region of promising designs, but it is a local approximation. Optimization should remain within the experimental or simulation domain unless additional validation supports extrapolation. Confirmation runs are essential: the predicted optimum should be evaluated independently, preferably with a higher-fidelity model or physical test.

5.7 Regularization, Bias, Variance, and Overfitting

A model with too few parameters may be unable to represent important behavior and therefore has high bias. A model with too many parameters may respond strongly to noise and therefore has high variance. The objective is not to minimize training error alone, but to minimize expected error on new data.

Ridge regression adds an L2 penalty:

min_θ ||Aθ-y||²_2+λ||θ||²_2.

Its solution is

θ=(A^TA+λI)^−1A^Ty.

The penalty shrinks coefficients and stabilizes nearly dependent systems. It is useful when predictors are correlated or when the number of parameters is large.

Lasso regression uses an L1 penalty:

min_θ ||Aθ-y||²_2+λ||θ||_1.

The L1 penalty can produce exactly zero coefficients and hence perform a form of variable selection. Elastic net combines L1 and L2 penalties and is useful when both sparsity and stability are desired.

Regularization does not remove the need for physical reasoning. Variables should be scaled before penalization, and the intercept is commonly left unpenalized. The parameter λ should be selected using validation, cross-validation, information criteria, or a physically justified error tolerance. A final test set must not be repeatedly used to tune the model.

5.8 A Practical Modelling Workflow

A reliable approximation study generally follows these stages:

1. Define the purpose, domain, accuracy requirement, and acceptable computational cost.

2. Inspect the data, units, ranges, missing values, outliers, and measurement uncertainty.

3. Select a basis or model guided by smoothness, periodicity, locality, and physical laws.

4. Scale variables and choose nodes or an experimental design that covers the domain.

5. Compute coefficients using stable methods such as QR or singular-value decomposition.

6. Evaluate residuals and errors using data not used for fitting whenever possible.

7. Test sensitivity to nodes, model degree, regularization, and influential observations.

8. Validate predictions with new experiments, simulations, or conservation and dimensional checks.

9. Document the domain of validity, assumptions, uncertainty, and extrapolation risks.

This workflow prevents a common error: treating an approximation algorithm as a substitute for model judgement.

5.9 Applied Case Study: Calibrating a Thermal Process

Background

Consider an illustrative manufacturing process in which the maximum temperature of a heat-treated component depends on heating power P and conveyor speed v. Sixteen design experiments are conducted over a specified operating region. The thermocouple measurements contain noise, and the engineering team needs a fast model for process adjustment.

Situation and Problem

A fifth-degree polynomial interpolant through selected measurements gives nearly zero error at those points but oscillates sharply near the limits of the operating region. When new measurements are added, the predicted optimum changes substantially. The team must obtain a stable approximation and determine whether an interior operating condition can meet the temperature target.

Analysis

The data are noisy, so exact interpolation is inappropriate. A second-order response surface is selected because the process is expected to be smooth over the restricted operating region. P and v are coded to [-1,1]. The design matrix includes 1, P, v, P v, P², and v². QR factorization is used to estimate the coefficients, avoiding direct formation of the normal equations.

Residual plots show a modest pattern with speed, suggesting that the quadratic model may be incomplete at one boundary. A third-order term is tested, but cross-validation indicates only a small predictive improvement and increased coefficient sensitivity. Ridge regression is also evaluated; a small positive penalty reduces sensitivity without materially increasing validation error.

Approach and Solution

The chosen model is the regularized quadratic response surface. The predicted target contour is used to identify a feasible operating region rather than a single mathematically optimal point. A confirmation run is conducted at the proposed setting and at two nearby settings. The results are compared with the response-surface predictions and with an independent higher-fidelity thermal simulation.

Outcome

The response surface provides a smooth and computationally inexpensive representation within the tested region. The confirmation results support its use for local process adjustment, but not for predicting temperatures at substantially higher powers or lower speeds. The team therefore records the calibration domain and requires new experiments before extending it.

Learning Points

The case illustrates why noisy engineering data should generally be fitted rather than interpolated, why coding variables improves conditioning, why validation must be separate from fitting, and why a response surface is a local decision tool rather than a universal physical law. It also demonstrates that a slightly more complex model is not automatically better if its predictions are less stable.

5.10 Challenges, Risks, and Professional Practice

Common mistakes include using a high-degree polynomial simply because many points are available, fitting a model before plotting the data, ignoring units and scaling, interpreting R² as proof of correctness, and extrapolating beyond the calibration domain. Duplicate or nearly duplicate nodes can make interpolation numerically unstable. Outliers should not be deleted automatically; their cause must be investigated.

Data leakage is a major validation risk. If the same observations are repeatedly used to select the model and report its performance, the reported error becomes optimistic. Time-dependent data require chronological rather than random splitting. Spatially correlated data may require spatial validation so that nearby points do not appear artificially independent.

Ethically and professionally, engineers must report uncertainty, limitations, and the conditions under which a model was developed. A fitted model should not conceal sensor failure, selectively exclude inconvenient observations, or be used outside its validated domain without disclosure. In safety-critical applications, approximation error should be incorporated into design margins and reviewed under applicable standards.

5.11 Digital Tools and Emerging Developments

Modern scientific-computing libraries provide stable interpolation, spline, optimization, and linear-algebra routines. Automatic differentiation can supply derivatives for sensitivity analysis and gradient-based calibration. Gaussian-process surrogate models provide predictions with uncertainty estimates and are increasingly used when high-fidelity simulations are expensive. Neural networks and other machine-learning models can approximate highly nonlinear responses, but they require careful validation, regularization, scaling, and uncertainty assessment.

Physics-informed and hybrid models combine data fitting with governing equations, conservation laws, or known symmetries. Such approaches can reduce data requirements and improve extrapolation, although they do not eliminate the need to check identifiability and numerical stability. The growing use of automated model selection makes transparency especially important: engineers should retain information about data provenance, preprocessing, model versions, validation results, and the domain of applicability.

Best Practices

• Match interpolation to trustworthy data and regression to noisy observations.
• Prefer local splines or carefully chosen nodes over unnecessarily high-degree global polynomials.
• Scale variables and use stable factorizations rather than explicit matrix inversion.
• Inspect residuals and validate with independent data or confirmation experiments.
• Use orthogonal bases when polynomial terms are strongly correlated.
• Select regularization strength through validation and report its effect.
• Treat response surfaces as local approximations with explicit domains of validity.
• Include uncertainty, sensitivity, and extrapolation risk in engineering decisions.
• Preserve physical constraints, dimensional consistency, and professional accountability.

Chapter Summary

Function approximation converts discrete, noisy, or computationally expensive information into a usable mathematical representation. Polynomial interpolation is exact at prescribed nodes, but its error and stability depend strongly on node placement and degree. Piecewise linear interpolation is simple, while cubic splines provide smooth local approximations and smoothing splines introduce a curvature penalty for noisy data.

Least-squares regression determines model parameters by minimizing residual error and supports weighted observations, multivariable models, and response surfaces. Orthogonal functions improve numerical conditioning and provide interpretable expansions. Response-surface models are valuable for local design exploration, provided that variables are scaled and the fitted surface is confirmed experimentally.

Approximation quality must be judged on more than training error. Bias, variance, conditioning, residual structure, validation error, and extrapolation risk all matter. Ridge, lasso, elastic-net, and smoothing penalties control overfitting by restricting model complexity. The central professional principle is to choose the simplest stable model that satisfies the engineering purpose and to document its uncertainty and domain of validity.

Key Terms

Interpolation: Construction of a function that passes exactly through specified data points.

Extrapolation: Prediction outside the range covered by the available data.

Approximation error: Difference between the true function and its computational representation.

Lagrange polynomial: Interpolating polynomial expressed as a weighted sum of cardinal basis polynomials.

Newton polynomial: Interpolating polynomial expressed using divided differences.

Divided difference: Recursive quantity used to construct Newton interpolation coefficients.

Runge phenomenon: Oscillatory behavior of high-degree interpolation, especially with equally spaced nodes.

Chebyshev nodes: Nonuniform nodes clustered near interval endpoints to improve polynomial stability.

Spline: Piecewise polynomial joined with specified continuity conditions.

Natural cubic spline: Cubic spline with zero second derivative at both endpoints.

Residual: Difference between an observed value and its fitted or predicted value.

Least squares: Parameter estimation by minimizing the sum of squared residuals.

Orthogonal basis: Set of functions whose pairwise inner products are zero.

Response surface: Approximate relationship between design variables and an engineering response.

Regularization: Addition of a penalty or constraint to control model complexity and instability.

Ridge regression: Least-squares regression with an L2 coefficient penalty.

Lasso regression: Regression with an L1 penalty that can produce sparse coefficients.

Overfitting: Excessive adaptation of a model to training data and noise.

Conditioning: Sensitivity of a numerical problem or solution to perturbations in data.

Cross-validation: Repeated partitioning of data to estimate predictive performance and tune models.

Review Questions

Short Answer Questions

1. Distinguish interpolation from regression in terms of how data points are treated.

2. What is the purpose of divided differences in Newton interpolation?

3. Why can equally spaced nodes cause difficulty in high-degree interpolation?

4. State the continuity conditions normally imposed on a cubic spline.

5. What effect does the regularization parameter have in ridge regression?

Long Answer Questions

1. Derive the Lagrange interpolation formula and explain its error term.

2. Explain the construction, boundary conditions, and computational advantages of cubic splines.

3. Develop the least-squares normal equations for a multivariable linear model and discuss stable solution methods.

4. Explain how orthogonal functions improve polynomial approximation and numerical conditioning.

5. Discuss the role of response-surface methodology in engineering design, including design, fitting, validation, and limitations.

Analytical / Application Questions

1. A temperature data set contains substantial sensor noise and is sampled at irregular operating points. Would you use a high-degree interpolating polynomial, a cubic spline, or a regularized regression model? Justify your choice and identify validation steps.

2. A fitted response surface has very low training error but gives inconsistent optima when one experiment is removed. Diagnose the likely causes and propose a revised modelling strategy.

3. A design matrix for a polynomial model is nearly rank deficient because two input variables are strongly correlated. Compare the usefulness of variable scaling, an orthogonal basis, QR factorization, ridge regression, and lasso regression.

22. Multiple Choice Questions with Answers

Question 1

Which statement best distinguishes interpolation from least-squares regression?

A. Interpolation always uses fewer parameters than regression.

B. Interpolation passes through specified data values, while regression minimizes a discrepancy that may not be zero at every point.

C. Regression is applicable only to linear functions.

D. Interpolation is used only for extrapolation.

Correct Answer: B. Interpolation passes through specified data values, while regression minimizes a discrepancy that may not be zero at every point.

Question 2

What is a principal cause of the Runge phenomenon?

A. Using cubic splines with too few intervals

B. Using equally spaced nodes in a high-degree global polynomial

C. Applying a small ridge penalty

D. Using weighted least squares

Correct Answer: B. Using equally spaced nodes in a high-degree global polynomial

Question 3

Which boundary condition defines a natural cubic spline?

A. The function is zero at both endpoints

B. The first derivative is zero at all interior nodes

C. The second derivative is zero at both endpoints

D. The third derivative is continuous at the endpoints

Correct Answer: C. The second derivative is zero at both endpoints

Question 4

Why is QR factorization generally preferred to explicitly solving the normal equations?

A. It always produces an exact interpolation

B. It avoids the conditioning problems associated with forming A-transpose-A

C. It eliminates the need for validation

D. It forces all coefficients to zero

Correct Answer: B. It avoids the conditioning problems associated with forming A-transpose-A

Question 5

In an orthogonal expansion, how is a coefficient commonly obtained?

A. By dividing the function by its maximum value

B. By projecting the function onto the corresponding basis function

C. By setting all other coefficients equal to one

D. By minimizing only the maximum residual

Correct Answer: B. By projecting the function onto the corresponding basis function

Question 6

In a quadratic response surface, the term beta12 x1 x2 represents:

A. An intercept

B. A pure curvature effect of x1

C. An interaction between x1 and x2

D. Measurement variance

Correct Answer: C. An interaction between x1 and x2

Question 7

What is the primary effect of ridge regularization?

A. It forces every coefficient to be exactly zero.

B. It shrinks coefficients and improves stability in correlated or ill-conditioned models.

C. It guarantees zero training error.

D. It permits reliable extrapolation under all conditions.

Correct Answer: B. It shrinks coefficients and improves stability in correlated or ill-conditioned models.

Question 8

A model has an excellent training error but poor validation error. Which interpretation is most appropriate?

A. The model may be overfitting the training data.

B. The data must be exact.

C. The model is necessarily physically correct.

D. Increasing model degree will always solve the problem.

Correct Answer: A. The model may be overfitting the training data.

Question 9

Why should engineering variables often be scaled before regularized regression?

A. Scaling makes all measurements identical.

B. Scaling ensures that the penalty treats coefficients more comparably across variables.

C. Scaling removes all measurement noise.

D. Scaling converts regression into interpolation.

Correct Answer: B. Scaling ensures that the penalty treats coefficients more comparably across variables.

Question 10

Which practice is most appropriate when using a response surface to locate an engineering optimum?

A. Use the predicted optimum without confirmation.

B. Extrapolate far beyond the design region to find a stronger optimum.

C. Confirm the predicted optimum with independent simulation or experiment.

D. Select the model solely by its training R-squared value.

Correct Answer: C. Confirm the predicted optimum with independent simulation or experiment.