Download Lecture Notes on Interpolation - Computational Methods | MCEN 3030 and more Study notes Mechanical Engineering in PDF only on Docsity!
Chapter 3
Interpolation
Interpolation is defined as a continuous representation of a discrete set of data points. Applications might be for differentiation or integration or simply estimating the value of the function between two data points. There are two distinct interpolation methods to obtain such a representation:
- standard - where the curve passes smoothly through the data points, and
- least squares - where a smooth curve passes sufficiently close to data points that have some uncer- tainty.
Note that with standard interpolation we actually pass a curve through the data. If data is from an experiment with some uncertainty, it is best to use the method of least squares which does not require the approximating function to pass through all the data points. Our use of interpolation will be primarily in deriving differentiation and integration formulas; thus we will use standard interpolation exclusively. When the interpolating function is evaluated outside the defining range, then it is called an extrapolation.
3.1 General Interpolation Problem
The general standard interpolation problem is defined as follows: Given:
- a set of data points xi, i = 0, 1 ,... , N ;
- a set of values fi ≡ f (xi), i = 0, 1 ,... , N ; and
- a set of basis functions bj (x), j = 0, 1 ,... , N ;
Find: coefficients aj , j = 0, 1 ,... , N , so that
f (xi) =
∑^ N
j=
aj bj (xi), i = 0, 1 ,... , N.
Note that for simplicity we have assumed linearity in the coefficients aj , and subsequently the above equation is just a system of linear algebraic equations which can be written as
fi = Bij aj ,
where Bij ≡ bj (xi). The set of basis functions bj (x) is still arbitrary at this point (as long as they are complete in the space of f ), but for convenience we impose two restrictions on them:
- they should be easy to evaluate, and
- they should be easy to differentiate and integrate.
The most common basis functions used that satisfy these restrictions are polynomials of various kinds.
14 CHAPTER 3. INTERPOLATION
3.2 Lagrange Interpolation
Suppose we have a set of N + 1 (non-equally spaced) data (xi, fi). By choosing the basis functions bj (x) = xj^ , j = 0, 1 ,... , N , we can construct a polynomial of degree N that passes through the data:
fi =
∑^ N
j=
aj xji = a 0 + a 1 xi + a 2 x^2 i + · · · + aN xNi , i = 0, 1 , · · · , N.
Here, we have N + 1 equations for the N + 1 unknowns a 0 , a 1 , · · · , aN. This procedure for finding the coefficients of the polynomial is not very attractive. It involves solving a system of algebraic equations which is generally ill-conditioned for large N (the resulting matrix is called a Vandermonde matrix). In practice, we will define the polynomial in an explicit way (as opposed to solving a system of equations). Define a polynomial of degree N associated with each point xj :
L( jN )(x) = αj (x − x 0 )(x − x 1 ) · · · (x − xj− 1 )(x − xj+1 ) · · · (x − xN ),
where αj is a constant to be determined. Let’s use the product notation for L( jN^ )
L( jN )(x) = αj
∏N
i= i 6 =j
(x − xi).
If x is equal to any of the data points except xj , L( jN )(x) = 0. Note
L( jN )(xj ) = αj
∏N
i= i 6 =j
(xj − xi).
Let
αj =
∏^ N
i= i 6 =j
(xj − xi)
− 1
,
so that
L( jN )(x) = (x − x 0 ) · · · (x − xj− 1 )(x − xj+1) · · · (x − xN ) (xj − x 0 ) · · · (xj − xj− 1 )(xj − xj+1) · · · (xj − xN )
Then
L (N ) j (xi) =^ δij^ =
0 if i 6 = j, 1 if i = j. ,
where δij is the Kronecker delta function. Next form the linear combination:
f (x) =
∑^ N
j=
L
(N ) j (x)fj^.
This is a polynomial of degree N because it is a linear combination of polynomials of degree N. It is the classic Lagrange polynomial. By construction, it passes through all the data points. For example, at x = xi f (xi) = L( 0 N )(xi)f 0 + L( 1 N )(xi)f 1 + · · · + L( jN )(xi)fj + · · · + L( NN )(xi)fN.
Since L (N ) j (xi) is equal to zero except for^ j^ =^ i, and^ L
(N ) i (xi) = 1, f (xi) = fi.
Thus, f (x) is the desired polynomial. Note that the polynomial interpolation is unique. That is, there is only one polynomial of degree N passing through a set of N + 1 points. The Lagrange polynomial is
16 CHAPTER 3. INTERPOLATION
- The curve is piecewise cubic, that is, the coefficients of the polynomial are different on each interval (xi, xi+1).
- The curve passes through the given data.
- With the exception of the end points, the first and second derivatives are continuous at the points xi.
From continuity of f ′′(x), using the piecewise linear Lagrange formula, we have
f (^) i′′ (x) ≃
xi+1 − x xi+1 − xi
f ′′(xi) +
x − xi xi+1 − xi
f ′′(xi+1). (3.1)
Integrating the equation twice and fitting to the data f (xi) and f (xi+1) results in
fi(x) =
(xi+1 − x)^3 6 hi
f ′′(xi) +
(x − xi)^3 6 hi
f ′′(xi+1 ) + (xi+1 − x)
[
f (xi) hi
hi 6
f ′′(xi)
]
(x − xi)
[
f (xi+1) hi
hi 6
f ′′(xi+1)
]
where hi = xi+1 −xi. If we can find f ′′(xi), the formula will be complete. This is obtained from continuity of f ′^ using two adjacent formulas. The result is
hi− 1 6
f ′′(xi− 1 ) + (hi− 1 + hi) 3
f ′′(xi) + hi 6
f ′′(xi+1 ) = f (xi+1) − f (xi) hi
f (xi) − f (xi− 1 ) hi− 1
This is a tridiagonal system for f ′′(xi) that is numerically easily solved. This equation cannot be applied at i = 0 or N , so that we only have N − 1 equations for N + 1 unknowns. We get two additional equations by making some assumptions and/or approximations in the end intervals. Some choices of boundary conditions are:
- Periodic: f ′′(x 0 ) = f ′′(xN − 1 ), f ′′(xN ) = f ′′(x 1 ).
- Parabolic: f ′′(x 0 ) = f ′′(x 1 ), f ′′(xN ) = f ′′(xN − 1 ).
- Natural (or free): f ′′(x 0 ) = f ′′(xN ) = 0.
- Cantilever: f ′′(x 0 ) = λf ′′(x 1 ), f ′′(xN ) = λf ′′(xN − 1 ), 0 ≤ λ ≤ 1.
For equally spaced intervals h, we can show that if fi are the exact values a smooth function y(x) takes at xi, then
y(x) − f (x) =
h^4 96 ymax(iv)(ξ), x 0 < ξ < xN.
Example: The simplest case requires 4 points (N = 3): two in the middle and two end points. Also we assume equally spaced intervals and cantilever end conditions. Then for i = 1
f ′′(x 0 ) + 4f ′′(x 1 ) + f ′′(x 2 ) =
h^2
[f (x 0 ) − 2 f (x 1 ) + f (x 2 )]
or
f 0 ′′ + 4f 1 ′′ + f 2 ′′ =
h^2 (f 0 − 2 f 1 + f 2 ).
Applying the end condition at i = 0,
(4 + λ)f 1 ′′ + f 2 ′′ =
h^2 (f 0 − 2 f 1 + f 2 ).
Similarly, at i = 2
f 1 ′′ + (4 + λ)f 2 ′′ =
h^2
(f 1 − 2 f 2 + f 3 ).
3.5. PARAMETRIC INTERPOLATION 17
These equations can be readily solved for f 1 ′′ and f 2 ′′. To make things simple, we will evaluate the function in the middle of the interval at x = (x 1 + x 2 )/2, so
f
x 1 + x 2 2
(f 1 + f 2 ) −
h^2 6
(f 1 ′′ + f 2 ′′ ),
[
8(5 + λ)
]
(f 1 + f 2 ) −
8 (5 + λ)
(f 0 + f 3 ).
For f (xi) = y(xi) = exi^ , 0 ≤ x ≤ 1:
λ = 0 =⇒ |f (0.5) − e^0.^5 | ≈ 5 × 10 −^3 ,
λ = 1 =⇒ |f (0.5) − e^0.^5 | ≈ 5 × 10 −^4. In most cases the spline gives a smooth fit to the data. In some cases, however, it wiggles too much. A cure for this is to put the spline in tension. The effect of tension is to replace equation (3.1) by
f (^) i′′ (x) − σ^2 fi(x) = xi+1 − x xi+1 − xi
[
f ′′(xi) − σ^2 f (xi)
]
x − xi xi+1 − xi
[
f ′′(xi+1 ) − σ^2 f (xi+1)
]
On integrating twice and requiring that the curve pass through the points, we obtain
fi(x) =
sinh [σ(xi+1 − x)] sinh [σ(xi+1 − xi)]
f ′′(xi) σ^2
(xi+1 − x) (xi+1 − xi)
[
f (xi) −
f ′′(xi) σ^2
]
sinh [σ(x − xi)] sinh [σ(xi+1 − xi)]
f ′′(xi+1) σ^2
(x − xi) (xi+1 − xi)
[
f (xi+1) −
f ′′(xi+1 ) σ^2
]
Then on applying continuity of the first derivative, we find the equation for f ′′(xi) to be [ 1 hi− 1
σ sinh σhi− 1
]
f ′′(xi− 1 ) σ^2
[
σ cosh σhi− 1 sinh σhi− 1
hi− 1
σ cosh σhi sinh σhi
hi
]
f ′′(xi) σ^2
[
hi
σ sinh σhi
]
f ′′(xi+1 ) σ^2
f (xi+1) − f (xi) hi
f (xi) − f (xi− 1 ) hi− 1
Again we have a tridiagonal system of equations which is easily solved. Generally, it is recommended that if the natural spline is not satisfactory, one should try the tension spline with the tension parameter set equal to unity as a first guess. If the results still wiggle too much, the tension σ can be increased. Endpoint conditions for the tension spline may be modeled after those for the natural spline by replacing f ′′(xi) by f ′′(xi) − σ^2 f (xi) in them.
3.5 Parametric Interpolation
None of the methods discussed is capable of interpolating curves having infinite slopes, such as in com- puting contours f (x, y) = constant. One way out of this difficulty is to use parametric interpolation. We regard both x and y as dependent on a third variable s. The choice of the parameter s is free as long as it increases monotonically along the curve. A simple choice is the arclength. Then we take s = s 0 = 0, s = s 1 = s 0 + [(x 1 − x 0 )^2 + (y 1 − y 0 )^2 ]^1 /^2 , s = s 2 = s 1 + [(x 2 − x 1 )^2 + (y 2 − y 1 )^2 ]^1 /^2 , and so on. The parametric curves x(s) and y(s) are treated by the usual spline method. For a closed curve, periodic end conditions are appropriate; for other cases, other end conditions may be applied.
3.6 Multidimensional Interpolation
The simplest methods in multidimensional interpolation are extensions of one-dimensional methods. For example, if a function of x and y is given on a regular array of points, then one can do the following. For each line x = xi on which the data are given, compute the second derivatives (which we call f (^) i′′ (yj )) of
