Function Interpolation
Marcin Chrząszcz, Danny van Dyk
mchrzasz@cern.ch, danny.van.dyk@gmail.com
Numerical Methods, 21 September, 2016
1
/
21Interpolation - graphical interpretation
⇛ Let’s assume we have a given function in a tabular form:
x i x 1 x 2 x 3 ... x n f i = f (x i ) f 1 f 2 f 3 ... f n
⇛ The x i are data points.
⇛ The problem is how to calculate the function values between the points (interpolation) or outside (extrapolation - see next lecture!) .
(MeV) mχ
0 2000 4000
(MeV)χσ
2 3 4 5 6 7 8 9
2
/
21Interpolation - formal formulation
⇛ On the interval: [a, b] ⊂ R and set of points: D = (x i , y i ), x i ̸ = x j ⇔ i ̸ = j.
⇛ The interpolation is to find a function F : [a, b] → R that satisfies:
∀ (x i , y i ) ∈ D : F (x i ) = y i
⇛ The above equation is called the interpolation equation.
⇛ The interpolation error:
ϵ(x) = f(x) − F (x), where f(x) is the function that is being interpolated.
⇛ Each interpolation method has different errors, each of it can be estimated.
3
/
21Linear Interpolation
⇛ Let’s say we have n + 1 interpolation points: (y i , x i ).
⇛ The linear interpolation is given with with the equation:
F (x) =
∑ n j=0
a j ψ j (x) = a 0 ψ 0 (x) + a 1 ψ 1 (x) + a 2 ψ 2 (x) + ... + a n ψ n (x),
where:
• a j - parameter,
• ψ j (x) - base function.
⇛ Now the only thing one needs to do is to find the a i coefficients. Now there are infinite number such solutions so we need to assume that the number interpolation points is equal to number of base functions.
4
/
21Linear Interpolation
⇛ Now to find all the a i one just needs to solve the linear equation system:
ψ 0 (x 0 ) ψ 1 (x 0 ) ... ψ n (x 0 ) ψ 0 (x 1 ) ψ 1 (x 1 ) ... ψ n (x 1 )
...
ψ 0 (x 0 ) ψ 1 (x 0 ) ... ψ n (x 2 )
a 0
a 1
...
a n
=
y 0
y 1
...
y n
⇛ The only necessary condition to solve the equation is that the:
det [ψ i (x j )] ̸ = 0
⇛ The polynomial of the order n has n + 1 base functions, so system has only one solution.
⇛ We can choose the base functions on many different ways. This depends on the input points and our own judgement. Often use are trigonometric functions and polynomials.
⇛ The proper selection of the base function can hugely simplify the problem by making the matrix triangular or even diagonal.
5
/
21Linear Interpolation, simplest base
. Theorem:
..
.
If all of the interpolation points x
0, x
1, ..., x
nare pair different there exists one and only one solution of the interpolation function of the order max. n.
The simplest base functions one can imagine to use is the natural base:
ψ 0 = 1, ψ 1 = x, ψ 2 = x 2 , ... , ψ n = x n
. Attention/Achtung:
..
.
This kind of interpolation is wrongly conditioned: small changes of input parameters will cause large output differences.
⇛ Because of the above it is used almost not used in practice.
⇛ For education purposes we will analyse it.
6
/
21Linear Interpolation, simplest base
⇛ When using natural base the interpolation polynomial will have the structure:
F(x) = a 0 + a 1 x + a 2 x 2 + a 3 x 3 ...a n x n
⇛ Now base matrix (so-called Vandermonde matrix) has the form:
V =
1 x 0 x 2 0 ... x n 0 1 x 1 x 2 1 ... x n 1
...
1 x n x 2 n ... x n n
⇛ Now the only thing one needs to do is to solve:
1 x 0 x 2 0 ... x n 0 1 x 1 x 2 1 ... x n 1
...
1 x n x 2 n ... x n n
a 0 a 1 ...
a n
=
y 0 y 1 ...
y n
⇛ Because the natural base is wrongly conditioned the increase of the interpolation points significantly increases the condition parameter.
7
/
21Lagrange interpolation
⇛ Now let’s choose a better base function that are given by:
l i (x) =
∏ n j=0, j
̸=i
x − x j x i − x j
= (x − x 0 )(x − x 1 )...(x − x n )
(x i − x 0 )(x i − x 1 )...(x i − x n ) (1)
⇛ It is easy to notice:
l i (x j ) = {
0 i ̸ = j 1 i = j
⇛ Now putting this into base matrix we get:
1 0 0 ... 0 0 1 0 ... 0
...
0 0 0 ... 1
a 0
a 1
...
a n
=
y 0
y 1
...
y n
⇛ Now this is something we like :) The interpolating function then takes the form:
F (x) = L n (x) =
∑ n i=0
y i l i (x)
8
/
21Lagrange interpolation - the algorithm
1. Get the input data (x i , y i ), i = 0, ..., n.
2. The coefficients of the polynomial are y i . 3. Calculate the base functions from Eq. 1.
4. Put everything together to get the L n Lagrange polynomial.
. Disadvantage:
..
.
The problem with this method is the fact that if you add another interpolation point you need to start over and recalculate everything.
. Advantege
..
.
It is very popular and simple method. It is used in many different places:
differential equations, numerical integrations, etc. Because of this the precision of the method depends on the precision of the interpolation is self:
ϵ(x) = f (x) − F(x) = f
(n+1(ξ) (n + 1)! ω
n(x)
9
/
21Lagrange interpolation - the algorithm
1. Get the input data (x i , y i ), i = 0, ..., n.
2. The coefficients of the polynomial are y i . 3. Calculate the base functions from Eq. 1.
4. Put everything together to get the L n Lagrange polynomial.
. Disadvantage:
..
.
The problem with this method is the fact that if you add another interpolation point you need to start over and recalculate everything.
. Advantege
..
. where:
ω
n(x) = (x − x
0)(x − x
1)...(x − x
n) f function of class C
n+1on [a, b]
ξ average values.h
9
/
21Newton interpolation
⇛ Sir Isaac Newton proposed a different base for the problem:
p 0 (x) = 1 p 1 (x) = (x − x 0 )
p 2 (x) = (x − x 0 )(x − x 1 ) ...
p n (x) = (x − x 0 )(x − x 1 )...(x − x n )
⇛ Using the above functions the basis matrix takes a nice form:
V =
1 0 0 ... 0
1 p 1 (x 1 ) 0 ... 0 ...
1 p 1 (x n ) p 2 (x n ) ... p n (x n )
10
/
21Newton interpolation
⇛ Sir Isaac Newton proposed a different base for the problem:
p 0 (x) = 1 p 1 (x) = (x − x 0 )
p 2 (x) = (x − x 0 )(x − x 1 ) ...
p n (x) = (x − x 0 )(x − x 1 )...(x − x n )
⇛ And the linear equation system:
1 0 0 ... 0
1 p 1 (x 1 ) 0 ... 0 ...
1 p 1 (x n ) p 2 (x n ) ... p n (x n )
a 0
a 1
...
a n
=
y 0
y 1
...
y n
10
/
21Newton interpolation
⇛ We can solve the aforementioned system by noticing that a 0 depends only on y 0 , a 1 depends on y 0 and y 1 , etc.
⇛ They can be calculated using difference quotient:
f [x 0 , x 1 ] = y 1 − y 0 x 1 − x 0
, f [x 1 , x 2 ] = y 2 − y 1 x 2 − x 1
, etc.
f [x 0 , x 1 , x 2 ] = f [x 1 , x 2 ] − f [x 0 , x 1 ] x 2 − x 1
⇛ In general:
f [x i , x i+1 , ..., x i+k ] = f [x i+1 , x i+2 , ..., x i+k ] − f [x i , x i+1 , ..., x i+k
−1 ] x i+k − x i
⇛ Using difference quotient we can calculate the coefficients:
a i = f [x 0 , x 1 , ..., x i ]
11
/
21Newton interpolation, the algorithm
• We read the input interpolating points: (x i , y i ), i = 0, 1, ...., n.
• Calculate the polynomial coefficients in an iterative way.
• Calculate the base functions.
• Put things together to get the Newton interpolation polynomial.
. Why do we bother?
..
.
⇛ Despite the procedure seems more complicated then the Lagrange interpolation it has a huge advantage: adding one more interpolation point does require to repeat the whole procedure, but we can reuse the previous steps.
12
/
21Runge oscillation
⇛ If one has many points one needs to use a high order polynomial.
⇛ The problem with high order polynomials is that they start oscilationg at the edges 7→ so-called Runge oscillations.
⇛ Danny will tell you more about this next lecture.
13
/
21First degree spline function
⇛ Till now using couple numerical methods we were able to calculate coefficients of a big single polynomial equation.
⇛ Now we are smarter then Newton and Lagrange were over 300 years ago!
⇛ We can use several interpolating functions and then glue them (splice) them together!
⇛ In order to do so we divide the domain of the F function [a, b] into:
a = x 0 < x 1 < x 2 < ... < x n = b
⇛ The simplest is to use linear functions:
s(x) =
s 0 (x) x ∈ [x 0 , x 1 ) s 1 (x) x ∈ [x 1 , x 2 ) ...
s n
−1 (x) x ∈ [x n
−1 , x n ) ,
where
s k (x) = a k,0 + a k,1 (x − x k ) x ∈ [x k , x k+1 )
14
/
21First degree spline function
⇛ The above defined functions s k (x) need to obey the conditions:
s k (x k ) = y k
s k (x k+1 ) = s k+1 (x k+1 )
⇛ From the above we can calculate the linear equation coefficients:
a k,0 = y k
a k,1 = y k+1 − y k
x k+1 − x k
. Advantage:
..
.
⇛ The biggest advantage of this method is it’s simplicity. If you have a large number of interpolating points that are ”densely ” packed you can get very good and fast aproximation with this function.
. Disadvantage:
..
. ⇛ The interpolating function is not smooth on the interpolation points.
15
/
21Third degree spline function
⇛ As mentioned before the first degree spline function suffer from the continuity problem.
⇛ A solution would be to use the higher order polynomials to aproximate the function. For example the 3rd order:
s k (x) = a k,0 + a k,1 (x − x k ) + a k,2 (x − x k ) 2 + a k,3 (x − x k ) 3 , x ∈ [x k , x k+1 )
⇛ Now this polynomial needs to obey conditions:
s k (x k ) = y k
s k (x k+1 ) = s k+1 (x k+1 ) s
′k (x k+1 ) = s
′k+1 (x k+1 ) s
′′k (x k+1 ) = s
′′k+1 (x k+1 )
⇛ In total we will have 4n coefficients to find. But we will only have 4n − 2 equations.
⇛ To solve the system we need to assume somethings about the endpoints (there are many options here).
⇛ The most popular ones are:
s
′(a) = f
′(a) s
′(b) = f
′(b) or
s
′′(a) = s
′′(b) = 0
16/
21Third degree spline function
⇛ Now to efficiently solve the system we will define temporary variables:
h k = x k+1 − x k
d k = y k+1 − y k
h k m k = s
′′(x k )
⇛ Now the second derivative can be interpolated using the Lagrange polynomial:
s
′′k (x) = s
′′(x) x − x k+1
x k − x k+1
+ s
′′(x) x − x k
x k+1 − x k
= m k
h k
(x k+1 − x) + m k+
h k
(x − x k )
⇛ Now if we integrate the above equation (two times):
s
′k (x) = − m k 2h k
(x k+1 − x) 2 + m k+
2h k
(x − x k ) 2 − p k + q k s k (x) = m k
6h k
(x k+1 − x) 3 + m k+1 3h k
(x − x k ) 3 + p k (x k+1 − x) + q k (x − x k )
17
/
21Third degree spline function
⇛ Now we need to calculate the s k (x) for x k and x k+1 using relations:
s k (x k ) = y k and s k (x k + 1) = y k+1 : x k (x k ) = y k = m k
6h k
(x k+1 − x k ) 3 + p k (x k+1 − x k ) x k (x k+1 ) = y k+1 = m k+1
3h k (x k+1 − x k ) 3 + q k (x k+1 − x k )
⇛ From which we get:
p k = y k
h k − m k h k 6 p k+1 = y k+1
h k+1 − m k+1 h k+1 6
⇛ Putting all the things together:
S k (x) = m k
6h k (x k − x) 3 + m k+1
6h k (x − x k ) 3 + ( y k
h k − m k h k
6 )
(x k+1 − x) (2) +
( y k+1
h k − m k+1 h k
6 )
(x − x k ) (3)
18
/
21Third degree spline function
⇛ The only unknown in the above equation are the m i variables. To get those we need to use the equation: S k
′−1 (x k ) = s
′(x k ):
− 1
3 m k h k − 1
6 m k+1 h k + d k = 1
3 m k h k
−1 + 1
6 m k
−1 h k
−1 + d k
−1
⇛ So the complete solution is:
a k,0 = y k
a k,1 = d k − h k
6 (2m k + m k+1 ) a k,2 = m k
2
a k,3 = m k+1 − m k
6h k
⇛ To get the m 0 and m n one needs to use one of the aforementioned conditions.
19
/
21Third degree spline function algorithm
• Get the input interpolation points.
• Calculate the temporary variables: d i , h i .
• Assume adequate conditions to get m 0 , m n .
• Solve linear equation system.
• Obtained values of the a i,j coefficients put in the interpolation equation.
20
/
21Summary
⇛ Interpolation playes essential role in almost all numerical methods!
⇛ Even very old algorithms like Lagrange are still used(we used it in splines).
⇛ There is a lot of algorithms on the market. From simple Lagrange and Newton interpolation algorithms up to modern splines.
⇛ The most common used today are spline. They require a bit of work but they work very effectively.
21
/
21Backup
22
