11.1.1 Polynomial interpolation

To begin with, consider the problem of determining a polynomial of degree one that passes through the distinct points (x₀, y₀) and (x₁, y₁). The linear polynomial

has the required property: it is easy to verify that P₁(x₀) = y₀ and P₁ (x₁) = y₁. To generalize the concept of linear interpolation to higher-degree polynomials, consider the construction of a polynomial of degree at most N that has the property

(11.1)

where points x_n and y_n are given. This polynomial is unique and may be represented in the form

(11.2)

with the factors

Polynomial (11.2) is called the Nth Lagrange interpolating polynomial. For an evaluation of the interpolating polynomial at a single point x without explicitly computing the coefficients of the polynomial, the following Neville scheme is very practical. Firstly we define

then interpolating polynomials are determined according to the recursive relation

(11.3)

and P₀^(N)(x) = P_N(x)

Representation (11.2) is very convenient for theoretical investigations because of its simple structure. However, for practical computations it is suitable only for small N. For large N the Lagrange factors l_n(x) become very large and highly oscillatory. There is another representation of the interpolating polynomial that is more practical for computational purposes. We first need to introduce divided differences. Assume that N + 1 distinct points x₀, …, x_N and N + 1 values y₀, …, y_N are given. The zeroth divided differences are simply the values of y_n:

The divided differences of order k at the point x_n are recursively defined by

With this notation, the uniquely determined interpolating polynomial with property (11.1) is given by

(11.4)

Polynomial (11.4) is called the Newton interpolating polynomial. This polynomial can also be written in a nested form:

Thus, the value of the Newton interpolating polynomial at a point x can be calculated using the following algorithm (the Horner scheme):

and P_N(x) = a₀.

As noted earlier, the Lagrange interpolating polynomial is a linear combination of Lagrange coefficients. Before implementing the function to calculate the polynomial, we implement the function to evaluate the Lagrange coefficient in the following. This function takes four arguments. The first argument is the value of x where we want to calculate the Lagrange coefficient. The second argument is a two column matrix with the first column containing the values of x and the second column containing y. The third argument, n, is the order of the coefficient and the last argument is the order of the polynomial.

In the following, we implement the function to evaluate Lagrange polynomial at a given value of x. The function takes similar parameter except we do not need the order of coefficients.

We apply the above function on the following set of data.

We apply the Lagrange polynomial method in the following to determine the curve representing this dataset.

11.1.2 Trigonometric interpolation

Periodic functions occur quite frequently in applications. These functions have the property f(x + T) = f(x) for some T > 0. For example, functions defined on closed planar or spatial curves may always be viewed as periodic functions. Polynomial interpolation is not appropriate for such functions, because algebraic polynomials are not periodic. Therefore, we turn to a consideration of interpolation by trigonometric polynomials. Without loss of generality we assume that the period T is equal to 2π. Assume that 2 N + l distinct points x₀, …, x_2N ∈ [0, 2π) and 2 N + 1 values y₀, …, y_2N are given. Then a uniquely determined trigonometric polynomial Q_n(x) exists and it has the property

In the Lagrange representation, this polynomial is given by

(11.5)

With the factors

When the interpolation points are equally spaced, representation (11.5) can be simplified. For an equidistant subdivision with an odd number 2 N + l of interpolation points

there exists a unique trigonometric polynomial

(11.6)

satisfying the interpolation property

Its coefficients are given by

For an equidistant subdivision with an even number 2 N of interpolation points

there exists a unique trigonometric polynomial

(11.7)

satisfying the interpolation property

Its coefficients are given by

Trigonometric interpolation polynomials (11.6) and (11.7) may be viewed as the discretized versions of the Fourier series, where the integrals giving the coefficients of the Fourier series are approximated by the rectangular rule at an equidistant partition.

An efficient numerical evaluation of trigonometric polynomials can be done analogously to the Horner scheme for algebraic polynomials. The recursion scheme has the form

(11.8)

starting with a_N = a_N and _N = 0. It delivers

If we use b_n instead of a_n and start with a_N = b_N and b_N = 0, then recursion scheme (11.8) delivers

Hence the evaluation of a trigonometric polynomial at a point x can be reduced to the evaluation of sin(x) and cos(x), and 6 N additions and 8 N multiplications.

The trigonometric polynomial is a linear combination of Lagrange like trigonometric coefficients. Before implementing the function to calculate the polynomial, we implement the function to evaluate the trigonometric coefficients in the following. The implementation is very similar to the implementation in the previous section. This function takes four arguments. The first argument is the value of x where we want to calculate the coefficient. The second argument is a two column matrix with the first column containing the values of x and the second column containing y. The third argument, n, is the order of the coefficient and the last argument is the order of the polynomial.

In the following, we implement the function to evaluate the trigonometric polynomial at a given value of x. The function takes a similar parameter except that we do not need the order of coefficients.

We apply the above function on the following set of data.

We apply the trigonometric polynomial method in the following and determine the curve representing this dataset.

11.1.3 Interpolation by splines

When a number of data points N becomes relatively large, polynomials have serious disadvantages. They all have an oscillatory nature, and a fluctuation over a small portion of the interval can induce large fluctuations over the entire range.

An alternative approach is to construct a different interpolation polynomial on each subinterval. The most common piecewise polynomial interpolation uses cubic polynomials between pairs of points. It is called cubic spline interpolation. The term spline originates from the thin wooden or metal strips that were used by draftsmen to fit a smooth curve between specified points. Since the small displacement w of a thin elastic beam is governed by the fourth-order differential equation d⁴w/dx⁴ = 0, cubic splines indeed model the draftsmen splines.

A cubic spline interpolant s(x) is a combination of cubic polynomials (s_n(x)(s(x) = ∪ s_n(x)), which have the following form on each subinterval [x_n, x_n + 1]:

(11.9)

To determine the coefficients of s_n(x), several conditions are imposed. These include

and

Additional conditions follow.

These conditions result in the system of linear equations for coefficients c_n of cubic polynomial (11.9):

where Δx_n = x_n + 1 − x_n The coefficient matrix of this system is a tridiagonal one, and it also has the property of diagonal dominance. Hence, the solution of this system can be obtained by the sweep method. Once the values of c_n are known, other coefficients of cubic polynomial s_n(x) are determined as follows:

11.2.1 Least-squares approximation

Consider the problem of estimating the values of a function at nontabulated points, given the experimental data in Table 11.1.

Interpolation requires a function that assumes the value of y_n at x_n for each n = 1, …, 11. Figure 11.1 shows a graph of the values in Table 11.1. From this graph, it appears that the actual relationship between x and y is linear. It is evident that no line precisely fits the data because of errors in the data collection procedure. In this case, it is unreasonable to require that the approximating function agrees exactly with the given data. A better approach for a problem of this type would be to find the “best” (in some sense) approximating line, even though it might not agree precisely with the data at any point.

Let fn = ax_n + b denote the nth value on the approximating line and y_n the nth given y-value. The general problem of fitting the best line to a collection of data involves minimizing the total error

with respect to the parameters a and b. For a minimum to occur, we need

and

These equations are simplified to the normal equations

(11.10)

The solution to this system is as follows:

Where

This approach is called the method of least-squares.

The problem of approximating a set of data (x_n,y_n), n = 1, …, N with an algebraic polynomial

of degree M < N-1 using least-squares is resolved in a similar manner. It requires choosing the coefficients a₀, …, a_M to minimize the total error:

Hence, the conditions ∂e/∂a_m = 0 for each m = 0, …, M result in a system of M + 1 linear equations for M + 1 unknowns a_m,

(11.11)

This system will have a unique solution, provided that the x_n are distinct. We apply the above method on the following set of data.

There are 7 data points (N = 7) and we will use a 3^rd order polynomial for least square approximation (M = 3). We build the following linear system of equations based on the least square method.

We solve the above system of equations in the following to find out the unknown values of a’s.

Once we know the polynomial coefficients, we could use the built-in polyval function to evaluate the polynomial for different values of the independent variable, x. We use this function in the following to generate the approximate polynomial values and create a plot showing the given data points and the approximating polynomial.

The approximating polynomial function nicely goes through the data points. The error could be reduced by increasing the polynomial order. But for the high order polynomial the curve shows oscillatory behavior so we could only use the low order polynomial.

11.2.2 Approximation by orthogonal functions

We now turn our attention to the problem of function approximation. The previous section was concerned with the least-squares approximation of discrete data. We now extend this idea to the polynomial approximation of a given continuous function f(x) defined on an interval [a,b]. The summations of previous section are replaced by their continuous counterparts, namely defined integrals, and a system of normal equation is obtained as before. For the least-squares polynomial approximation of degree M, the coefficients in

are chosen to minimize

To determine these coefficients we need to solve the following system of linear equations

(11.12)

These equations always possess a unique solution. The (n,m) entry of the coefficient matrix of system (11.12) is

Matrices of this form are known as Hilbert matrices and they are very ill conditioned. Numerical results obtained using these matrices are sensitive to rounding errors and must be considered inaccurate. In addition, the calculations that were performed in obtaining the Nth degree polynomial do not reduce the amount of work required to obtain polynomials of higher degree.

These difficulties can, however, be avoided by using an alternative representation of an approximating function. We write

(11.13)

where φ_m(x) (m = 0, …, M) is a function or polynomial. We show that, by proper choice of φ_m(x) , the least-squares approximation can be determined without having to solve a system of linear equations. Now the values of a_m(x) (m = 0, …, M) are chosen to minimize

where ρ(x) is called a weight function. It is a nonnegative function, integrable on [a,b], and it is also not identically zero on any subintervals of [a,b]. Proceeding as before gives

(11.14)

This system will be particularly easy to solve if the functions _m(x) are chosen to Satisfy

(11.15)

In this case, all but one of the integrals on the left-hand side of equation (11.14) are zero, and solution can be obtained as

Functions satisfying (11.15) are said to be orthogonal on the interval [a,b] with respect to the weight function ρ(x). Orthogonal function approximations have the advantage that an improvement of the approximation through addition of an extra term a_M + 1φ_M + 1(x) does not affect the previously computed coefficients a₀, …, a_M. Substitution of

(11.16)

where x∈ [a, b] and in equations (11.13) to (11.15) yields a rescaled and/or shifted expansion interval.

Orthogonal functions (polynomials) play an important role in mathematical physics and mechanics, so approximation (11.13) may be very useful in various applications. For example, the Chebyshev polynomials can be applied to the problem of function approximation. The Chebyshev polynomials T_n(x) are orthogonal on [-1, 1] with respect to the weight function ρ(x) = (1-x²)^− 1/2. For x∈ [–1, 1] they are defined as T_n (x) = cos(n arccos(x)) for each n≥ 0 or

and

The first few polynomials T_n(x) are

The orthogonality constant in (11.15) is a_n = π/2. The Chebyshev polynomial T_N(x) of degree N≥ 1 has N simple zeros in [-1, 1] at

(11.17)

for each n = 1, …, N. The error of polynomial approximation (11.13) is essentially equal to the omitted term a_M + 1T_M + 1(x) for x in [-1,1]. Because T_M + 1(x) oscillates with amplitude one in the expansion interval, the maximum excursion of the absolute error in this interval is approximately equal to |a_M + 1|.

Trigonometric functions are used to approximate functions with periodic behavior, that is, functions with the property f(x + T) = f(x) for all x and T > 0. Without loss of generality, we assume that T = 2π and restrict the approximation to the interval [–π,π]. Let us consider the following set of functions

Functions from this set are orthogonal on [-π,π] with respect to the weight function ρ(x) = 1, and α_k are equal to one. Given f(x) is a continuous function on the interval [-π,π], the least square approximation (11.16) is defined by

where

As cos(mx) and sin(mx) oscillate with an amplitude of one in the expansion interval, the maximum absolute error, that is max |f(x) − g(x)|, may be estimated by

The following function finds the coefficients for Chebyshev polynomial for a given input function argument and the degree of the polynomial.

The following function evaluates the Chebyshev polynomial for a given input vector containing polynomial coefficients and the value of x where the polynomial needs to be calculated.

We use the last two function implementations in the following program to determine the Chebyshev approximation of the exponential function e^x for x ∈ [-1, 1].

This program generates the following plot which shows some of the points of the exponential function along with the curve obtained by the approximated Chebyshev polynomial. It is interesting to see how the curve neatly goes through the exponential function points.

11.2.3 Approximation by interpolating polynomials

The interpolating polynomials may also be applied to the approximation of continuous functions f(x) on some interval [a,b]. In this case N + 1 distinct points x₀, …, x_N∈ [a,b] are given, and y_n = f(x_n). First, let us assume that f(x) is defined on the interval [-1, 1], If we construct an interpolating polynomial, which approximates function f(x), then the error for this approximation can be represented in the form

where P_N(x) denotes the Lagrange interpolating polynomial. There is no control over z(x), so minimization of the error is equivalent to minimization of the quantity (x-x)(x-x₁)…(x-x_N). The way to do this is through the special placement of the nodes x₀,…, x_N, since now we are free to choose these argument values. Before we turn to this question, some additional properties of Chebyshev polynomial should be considered.

The monic Chebyshev polynomial (a polynomial with leading coefficient equal to one) T_n^⁎(x) is derived from the Chebyshev polynomial T_n(x) by dividing its leading coefficients by 2ⁿ⁻¹, where n≥ l. Thus

As a result of the linear relationship between T_n^⁎(x) and T_n(x), the zeros of T_n^⁎(x) are defined by expression (11.17). The polynomials T_n^⁎(x), where n≥ 1, have a unique property: they have the least deviation from zero among all polynomials of degree n. Since (x-x₀)(x-x₁)…(x-x_N) is a monic polynomial of degree N + 1, the minimum of this quantity on the interval [-1, 1] is obtained when

Therefore, x_n must be the (n + l)st zero of T_n + 1^⁎(x), for each n = 0. …, N; that is, x_n = x_{N + 1,n + 1}^⁎ This technique for choosing points to minimize the interpolation error can be extended to a general closed interval [a,6] by choosing the change of variable

(11.18)

to transform the numbers x_n = x_{N + 1,n + 1}^⁎ in the interval [-1, 1] into corresponding numbers in the interval [a,b].

(a) Calculate coefficients of Lagrange interpolating polynomial that passes through the given points x_n and y_n (tables 11.2-11.7). Plot the polynomial together with initial data.

(b) Calculate coefficients of trigonometric interpolating polynomial that passes through the given points x_n and y_n (tables 11.8-11.13). Plot the polynomial together with initial data.

(c) Calculate coefficients of the Lagrange interpolating polynomial that approximates given function on some interval. Plot the function and polynomial. Suggested functions and intervals are

(1). Bessel function J₀(x), x ∈ [0,4],

(2). Bessel function J1(x), x ∈ [0,3],

(3). Errorfunction erf

(4). Function , x ∈ [0,1,2],

(5). Inverse hyperbolic sinus arcsh(x), x ∈ [− 5, 5],

(6). nverse hyperbolic cosine arcth(x), x ∈ [− 0.9, 0.9].