Numerical Analysis : [ MA214 ]

Lecture 3

Instructor: Prof. Tony J. Puthenpurakal

Department of Mathematics
Indian Institute of Technology Bombay
[email protected]

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

Recall
Let f : [a, b] −→ R be a function and let x0 , x1 , · · · , xn be distinct points in
[a, b] Then there exists a unique polynomial Pn (x) such that

Pn (xi ) = f (xi ) ,∀i = 0, 1, · · · , n

Pn (x) → Interpolating polynomial of f ( w .r .t. x0 , x1 , · · · , xn )

we had two forms of Pn (x)

Lagrange’s form
n
Y x − xi
lk (x) = , k = 0, 1, · · · , n
xk − xi
i=0,i6=k
Xn
Pn (x) = f (xi )li (x)
i=0

Advantage of Lagrange’s form:- It is simple.

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3
Disvantage of Lagrange’s Form:-
There is no way of using Lagrange form of Pn (x) to determine Pn+1 (x).
Newtons form of interpolating polynomial
Pn (x) = f (x0 ) + f [x0 , x1 ](x − x0 )
+ f [x0 , x1 , x2 ](x − x0 )(x − x1 )
+ ···
..
.
n−1
Y
+ f [x0 , x1 , x2 , · · · , xn ] (x − xi )
i=0
n−1
Y
Pn (x) = Pn−1 (x) + f [x0 , x1 , · · · , xn ] (x − xi )
i=0
f (x1 ) − f (x0 )
Here f [x0 , x1 ] =
x1 − x0
f [x1 , x2 , · · · , xk ] − f [x0 , x1 , · · · , xk−1 ]
f [x0 , x1 , x2 , · · · , xk ] =
xk − x0
Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3
Recall:- Error of interpolating polynomial

en (x) = f (x) − Pn (x)

n
Y
en (x) = f [x0 , x1 , · · · , xn , x] (x − xi )
i=0
f (n+1)(ξ)
f [x0 , x1 , · · · , xn , x] =
(n + 1)!
Yn
Ψn+1 (x) = (x − xj )
j=0
Not in syllabus on how to
minimize Ψn+1 (x)

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

Osculatory interpolation
Sometimes we have
x0 , x1 , · · · , xk
f (x0 ), f (x1 ), f (x2 ), · · · , f (xk )
f 0 (x0 ), f 0 (x1 ), f 0 (x2 ), · · · , f 0 (xk )

We need a polynomial p(x) such that

p(xi ) = f (xi ) and p 0 (xi ) = f 0 (xi ) , for i = 0, 1, · · · , n

Note:- deg p(x) ≤ 2n + 1

Example:-
dy
= g (x, y ) , y (x0 ) = y0
dx
Note:- y 0 (xi ) = g (xi , yi ) can be readily computed.

Last time I had given an example to calculate P(x)

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3
Algorithm for computing P(x)

Set z0 = x0 , and z1 = x0
z2 = x1 , and z3 = x1
..
.
z2i = xi , and z2i+1 = xi
..
.
z2n = xn , and z2n+1 = xn
( )
f 0 (zi ) if zi = zi+1
Set f [zi , zi+1 ] = f (zi+1 )−f (zi )
zi+1 −zi if zi 6= zi+1
For higher order computation
f [zi+1 , · · · , zi+k ] − f [zi , · · · , zi+k−1 ]
f [zi , zi+1 , · · · , zi+k ] =
zi+k − zi
P2n+1 (x) = f [z0 ] + f [z0 , z1 ](x − z0 ) + f [z0 , z1 , z2 ](x − z0 )(x − z1 )
2n
Y
+ · · · + f [z0 , · · · , z2n+1 ] (x − zj )
Instructor: Prof. Tony J. Puthenpurakal j=0 Analysis : [ MA214 ] Lecture 3
Numerical
Problem

f (0) = 1, and f (1) = 1.9, f 0 (0) = 2, and f 0 (1) = 2.5

Approximate f (0.4)

Set z0 = 0, z1 = 0, z2 = 1, z3 = 1

z f (z) f [, ] f [, , ] f [, , , ]
0 1 2 −1.1 2.7
0 1 0.9 1.6
1 1.9 2.5
1 1.9

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

Here

f [0, 0] = f 0 (0) = 2
f [1, 1] = f 0 (1) = 2.5
f (1) − f (0) 1.9 − 1
f [0, 1] = = = 0.9
1−0 1
f [0, 1] − f [0, 0]
f [0, 0, 1] = = 0.9 − 2 = −1.1
1−0
f [1, 1] − f [0, 1]
f [0, 1, 1] = = 2.5 − 0.9 = 1.6
1−0
finally
f [0, 1, 1] − f [0, 0, 1]
f [0, 0, 1, 1] = = 1.6 − (−1.1) = 2.7
1−0
P3 (x) = 1 + 2(x − 0) + (−1.1)(x − 0)2 + 2.7(x − 0)2 (x − 1)
= 1 + 2x − 1.1x 2 + 2.7x 2 (x − 1)
f (0.4) ≈ P3 (0.4) = 1.365

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

Theory of Osculatory Interpolation

Convention:-

Let x0 , x1 , · · · , xm be not necessarily distinct points.

We say two function f (x) and g (x) agree at the points x0 , x1 , · · · , xm if

f (j) (z) = g (j) (z) , for j = 0, 1, · · · , k − 1

for every points z which occurs k times in the sequence x0 , x1 , · · · , xm .

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

Example:-
x = 1, 2, 2, 1, 3, 2
f (x) and g (x) agree at the points 1, 2, 2, 1, 3, 2. If

f (1) = g (1)
f 0 (1) = g 0 (1)
f (2) = g (2)
f 0 (2) = g 0 (2)
f 00 (2) = g 00 (2)
f (3) = g (3)

Problem:-
Given x0 , x1 , · · · , xm not necessarily distinct points and f : [a, b] → R.

We need a polynomial p(x) of degree ≤ m such that p(x) and f (x) agree
at x0 , x1 , · · · , xm .
Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3
Exercise:- Show that if two polynomials of degree ≤ m agree at
x0 , x1 , · · · , xm , then they are equal.

So it makes sense to talk about the polynomial of degree ≤ m which

agrees with f (x) at the m + 1 points x0 , x1 , · · · , xm .
Theorem
If f (x) has r continuous derivatives and no point in the x0 , x1 , · · · , xm
occurs more than r times then there exists exactly one polynomial Pm (x)
of degree ≤ m which agree with f (x) at x0 , x1 , · · · , xm .

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

Proof

Uniqueness-Exercise

Existence:-
Assume x0 ≤ x1 ≤ · · · ≤ xm
for m = 0 there is nothing to show.

Assume the statement is true for m = k − 1 and consider it for m = k.

There are two cases

Case 1 x0 = xk . Then x0 = x1 = x2 = · · · = xk .
So r ≥ k. By assumption f has at least k continuous derivatives.

Then the Taylor polynomial Pk (x) for f (x) around the center c = x0 does
the job.

f (k) (x0 )
Remark:- Note that its leading coefficient is k!

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

Continue proof

Case 2 x0 < xk

Then by induction hypothesis we can find polynomial Pk−1 (x) of degree

≤ k − 1 which agree with f (x) at x0 , x1 , · · · , xk and a polynomial qk−1 (x)
of degree ≤ k − 1 which agree with f (x) at x1 , x2 , · · · , xk .

Verify
x − x0 xk − x
Pk (x) = qk−1 (x) + Pk−1 (x)
xk − x0 xk − x0
does the job.

[ Slightly tricky to show. See textbook Conte and de Boor page 64. ]

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

Convention
x0 , x1 , · · · , xm not necessarily distinct points

Pm (x) = unique polynomial which agree with f (x) at x0 , x1 , · · · , xm .

f [x0 , x1 , · · · , xm ] = leading coefficient of Pm (x) = coefficient of x m in

Pm (x).

We have
m−1
Y
Pm (x) = Pm−1 (x) + f [x0 , x1 , · · · , xm ] (x − xi )
i=0

Proof:-
m−1
Y
Pm (x) − f [x0 , x1 , · · · , xm ] (x − xi )
i=0
has degree ≤ m − 1 and agrees with f (x) at x0 , x1 , · · · , xm .

So by uniqueness of interpolating polynomial the result follows.

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3
Thus we can write Pm (x) as
m−1
Y
Pm (x) = f (x0 ) + f [x0 , x1 ](x − x0 ) + · · · + f [x0 , x1 , · · · , xm ] (x − xi )
i=0

Note we are assuming x0 ≤ x1 ≤ · · · ≤ xm .

Case 1 x0 = x1 = x2 = · · · = xm .

Pm (x) = Taylor polynomial with center x0

f 00 (x0 )
= f (x0 ) + f 0 (x0 )(x − x0 ) + (x − x0 )2
2!
f (m) (x0 )
+··· + (x − x0 )m
m!
f (m) (x0 )
So f [x0 , x0 , · · · , x0 ( m + 1 times ) ] =
m!

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

Case 2 xm 6= x0 . Then
f [x1 , x2 , · · · , xm ] − f [x0 , x1 , · · · , xm−1 ]
f [x0 , x1 , x2 , · · · , xm ] =
xm − x0

Example:- f (0) = 1, f 0 (0) = 0, f 00 (0) = 1, f (0.1) = 9.95E −1.

Approximate f (0.05).

z f (z) f [, ] f [, , ] f [, , , ]
0 1 0 0.5 −10
0 1 0 −5E −1
0 1 −5E −2
0.1 9.95E −1

P(x) = 1 + 0(x − 0) + (0.5)(x − 0)2 + (−10)(x − 0)3

= 1 + 0.5x 2 − 10x 3
P(0.05) ≈ 1 + 0.5(0.05)2 − 10(0.05)3 = 1

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

Calculations

z0 = 0, z1 = 0, z2 = 0, z3 = 0.1

f [z0 , z1 ] = f 0 (0) = 0
f [z1 , z2 ] = f 0 (0) = 0
f (z3 ) − f (z2 )
f [z2 , z3 ] = = −5E −2
z3 − z2
f 00 (z0 )
f [z0 , z1 , z2 ] = = 0.5
2
f [z2 , z3 ] − f [z1 , z2 ] −5E −2 − 0
f [z1 , z2 , z3 ] = = = −5E −1
z3 − z1 0.1
f [z1 , z2 , z3 ] − f [z0 , z1 , z2 ]
f [z0 , z1 , z2 , z3 ] = = −10
z3 − z0

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

Examples where Osculatory interpolation is useful
dy
1. dx = f (x, y ), and y (x0 ) = y0

Note that numerical methods for solving O.D.Es yield

x y
x1 y (x1 )
x2 y (x2 )
.. ..
. .
xn y (xn )

Note that y 0 (xi ) = f (xi , yi ) can be calculated easily.

Rx
2. f (x) = 0 g (t)dt

Numerical integration techniques yield f (x0 ), f (x1 ), · · · , f (xn ).

Note f 0 (x) = g (x), f 0 (xi ) = g (xi ) can be calculated easily.

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

As the following example will be used often I give a direct proof.

Example:- a, b distinct points. We know f (a), f (b), f 0 (a), f 0 (b).

P3 (x) = f (a) + f [a, a](x − a) + f [a, a, b](x − a)2 + f [a, a, b, b](x − a)2 (x − b)

We prove by direct computation that P3 (x) agree with f (x) at a, a, b, b

P3 (a) = f (a)
P30 (a) = f [a, a] = f 0 (a)
f (b) − f (a)
f [a, b] =
b−a
f [a, b] − f [a, a]
f [a, a, b] =
b−a
f (b)−f (a)
b−a − f 0 (a)
=
b−a
f (b) − f (a) − (b − a)f 0 (a)
=
(b − a)2

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

 f (b) − f (a) − (b − a)f 0 (a)
P3 (b) = f (a) + f 0 (a)(b − a) + (b − a)2
(b − a)2
= f (b)

Verify

(b − a)f 0 (b) − (f (b) − f (a))

f [a, b, b] =
(b − a)2
(b − a)(f 0 (b) + f 0 (a)) − 2(f (b) − f (a))
f [a, a, b, b] =
(b − a)3
P30 (b) = f 0 (a) + 2f [a, a, b](b − a) + f [a, a, b, b](b − a)2
= f 0 (b) ( Check ? )

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

Continuity of divided differences

Theorem
f [x0 , x1 , · · · , xn ] is a continuous function of x0 , x1 , · · · , xn . (Assume f has
n continuous derivatives )
(r ) (r )
i.e. for each r , x0 , · · · , xn are n + 1 points in [a, b] and
(r )
lim x = yi for i = 0, 1, · · · , n
r →∞ i

Then
(r ) (r ) (r )
lim f [x0 , x1 , · · · , xn ] = f [y0 , y1 , · · · , yn ]
r →∞

Proof.
See textbook Conte and de Boor page 65.

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

Disadvantages of Interpolation
Note that if x0 , x1 , · · · , xk are points in [a, b] then the interpolating
polynomial has degree k.

In practice k is large. Furthermore a polynomial of degree k with k large

Oscillates a lot

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

For example if there are 101 points then it is not advisable to work with a
degree 100 interpolating polynomial as this also creates lot of round-off
error.

Strategy:-
Use piecewise-polynomial approximation

Simplest case Piecewise-liner interpolation:-

a = x1 < x2 < x3 < · · · < xn = b

f (x) is approximated at a point x by first locating the interval [xk , xk+1 ]
containing x and then taking

p(x) = f (xk ) + f [xk , xk+1 ](x − xk )

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

Graphical representation

Note that if N is large and so |xi − xi+1 | is small then this is a good
approximation of f (x).

If we use higher degree ( say cubic ) piecewise-polynomial approximation

then we get better approximation.
Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3
Construction of piecewise-cubic function at the points x1 , · · · , xN+1 where

a = x1 < x2 < x3 < · · · < xN+1 = b

On each [xi , xi+1 ] we construct g3 (x) as a cubic polynomial Pi (x),

i = 1, 2, · · · , N

Pi (x) = c1,i + c2,i (x − xi ) + c3,i (x − xi )2 + c4,i (x − xi )3

for i = 1, 2, · · · , N. Since g3 (xi ) = f (xi ), for i = 1, 2, · · · , N + 1

We have
Pi (xi ) = f (xi ) and Pi (xi+1 ) = f (xi+1 )
for i = 1, 2, · · · , N

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

In particular

Pi−1 (xi ) = Pi (xi ) = f (xi ), for i = 1, 2, · · · , N

So g3 (x) is continuous on [a, b].

Only constraints for Pi (x) is

Pi (xi ) = f (xi ) and Pi (xi+1 ) = f (xi+1 )

So we have some freedom in choosing Pi (x)

We study 2 cases
1 Piecewise-Cubic Hermite interpolation
2 Cubic-Spline interpolation

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

Piecewise-Cubic Hermite interpolation

One determine Pi (x) so as to interpolate f (x) at xi , xi , xi+1 , xi+1 ,

i.e. we also have

Pi0 (xi ) = f 0 (xi ) and Pi0 (xi+1 ) = f 0 (xi+1 )

By Newtone formula

Pi (x) = f (xi ) + f [xi , xi ](x − xi ) + f [xi , xi , xi+1 ](x − xi )2

+f [xi , xi , xi+1 , xi+1 ](x − xi )2 (x − xi+1 )
Write x − xi+1 = (x − xi ) + (xi − xi+1 )
Pi (x) = f (xi ) + f 0 (xi )(x − xi )
+ (f [xi , xi , xi+1 ] − f [xi , xi , xi+1 , xi+1 ]4xi )(x − xi )2
+ f [xi , xi , xi+1 , xi+1 ](x − xi )3

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

Algorithim

For i = 1, 2, · · · , N + 1

fi = f xi
4xi = xi+1 − xi
si = f 0 (xi ), c1,i = fi , c2,i = si

c3,i = f [xi , xi , xi+1 ] − f [xi , xi , xi+1 , xi+1 ]4xi

f [xi , xi+1 ] − si
= − c4,i 4xi
4xi
c4,i = f [xi , xi , xi+1 , xi+1 ]
f [xi , xi+1 , xi+1 ] − f [xi , xi , xi+1 ]
=
4xi
si+1 + si − 2f [xi , xi+1 ]
=
4xi

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

Example

dy
= y − x 2 + 1, and y (0) = 0.5, 0 ≤ x ≤ 1
dx
xi y (xi ) y 0 (xi )
0 0.5 1.5
0.2 0.826 1.786
0.4 1.207 2.047
0.6 1.637 2.277
0.8 2.110 2.470
1.0 2.618 2.618

Find y (0.7), y (0.9)

Note:- y (xi ) is found using a Numerical method

Remarks:- Usual oscillatory polynomial has degree 11.

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

So we use piecewise Hermite interpolation

x y (x) f [, ] f [, , ] f [, , , ]
0.6 1.637 2.277 4.4E −1 4.25E −1
0.6 1.637 2.365 5.25E −1
0.8 2.110 2.470
0.8 2.110

P4 (x) = 1.637 + 2.277(x − 0.6) + 4.4E −1(x − 0.6)2

+ 4.25E −1(x − 0.6)2 (x − 0.8)
P4 (0.7) = 1.869

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3

 x y (x) f [, ] f [, , ] f [, , , ]
0.8 2.110 2.47 0.35 0.2
0.8 2.110 2.54 0.39
1.0 2.618 2.618
1.0 2.2.618

P4 (x) = f [0.8] + f [0.8, 0.8](x − 0.8) + f [0.8, 0.8, 1](x − 0.8)2

+ f [0.8, 0.8, 1, 1](x − 0.8)2 (x − 1)
= 2.11 + 2.47(x − 0.8) + 0.35(x − 0.8)2 + 0.2(x − 0.8)2 (x − 1)
P5 (0.9) = 2.360 in 4 sig digits

Instructor: Prof. Tony J. Puthenpurakal Numerical Analysis : [ MA214 ] Lecture 3