Interpolation with high degree polynomials

Numerical Analysis 2, lecture 4: is risky

f
Interpolating with splines
x0 f0
x1 f1
! ! x
(textbook sections 5.9–11)
x n fn

P(xi ) = fi ! P(x) " f (x)

• linear spline interpolation polynomials

• cubic spline • Vandermonde, Lagrange, Newton, Neville
• error formula f ( x̂) ! P( x̂) = 1 ( x̂ ! x0 )( x̂ ! x1 )!( x̂ ! xn ) f (n+1) (" )
(n + 1)!
• Runge’s example
2
si 1 1
s1 si+1 sn 0 f (x) =
1 + 25x 2
2

0
x0 x1 xi-1 xi xi+1 xn-1 xn
2

0
!1 0 1
Numerical Analysis 2, lecture 4, slide! 2

A linear spline is a continuous piecewise The error of linear spline interpolation

degree-1 polynomial is O( h2 )

s1 si
si+1 sn error in ith piece
f (x) ! si (x) " 18 (xi ! xi!1 ) 2 max f $$(# )
!#"# $ xi!1 "# " xi
hi
x0 x1 xi-1 xi xi+1 xn-1 xn knots f !(x) " si! (x) # hi max f !!($ )
xi"1 #$ # xi

Each si := s [ x is a polynomial of degree " 1 and s is continuous

i!1 , xi )
global error
max f (x) " s(x) ! 18 (max hi )2 max f ##(x)
x0! x! xn
!"# x0 ! x! xn
h

interpolating linear spline max

x0! x! xn
f "(x) # s"(x) ! h max
x0 ! x! xn
f ""(x)

si (x) = bi ! (x " xi"1 ) + ai

& ai = fi!1 variational characterization

si (xi!1 ) = fi!1 " (
# % ' fi ! fi!1 (i = 1, 2,…, n) xn
si (xi ) = fi $ ( bi = x ! x " ( g!(x))
2
Among all interpolants g, dx is minimized by g = s
) i i!1
x0

Numerical Analysis 2, lecture 4, slide! 3 Numerical Analysis 2, lecture 4, slide! 4

 Linear interpolating splines A cubic spline has continuous 2nd derivative
are easily computed in Matlab and degree-3 polynomial pieces
n = 10;
si
x = linspace(-1,1,n+1)';
" ai = fi!1 s1 si+1 sn
f = 1./(1+25*x.^2); $
coef = [diff(f)./diff(x) f(1:n)]; coefficients of spline pieces # fi ! fi!1
s = mkpp(x,coef); $ bi = x ! x
xx = linspace(-1,1,101); % i i!1
ff = 1./(1+25*xx.^2);
find the interval and x0 x1 xi-1 xi xi+1 xn-1 xn
ss = ppval(s,xx);
plot(xx,ss,'r',x,f,'ko') evaluate the polynomial (Horner)
plot(xx,ss-ff,'r',x,f-f,'ko') Each si := s [ x is a polynomial of degree " 3 and s## is continuous
i!1 , xi )
1 1

this doesn’t uniquely define it:

0 0 x ! xi!1
!1 0 1 !1 0 1 si (x) = di u 3 + ci u 2 + bi u + ai where u = (i = 1, 2,…, n) 4n coefficients
xi ! xi!1
si (x) = bi ! (x " xi"1 ) + ai
si (xi!1 ) = fi!1, si (xi ) = fi (i = 1, 2,…, n) 2n equations
0.1 0.1

si! (xi ) = si+1

! (xi ), si!!(xi ) = si+1
!! (xi ) (i = 1, 2, …, n " 1) 2n-2 equations
0 0

!0.1 !0.1 …we need 2 more conditions

!1 0 1 !1 0 1
Numerical Analysis 2, lecture 4, slide! 5 Numerical Analysis 2, lecture 4, slide! 6

The two extra conditions for cubic splines Recall the cardinal functions
can be chosen in many ways for Hermite interpolation (lecture 6)

“natural” boundary conditions

2 2
s1!!(x0 ) = 0, sn!!(xn ) = 0 " x ! x0 % " x ! x1 % " x ! x1 % " x ! x0 %
P(x) = $ 1 ! 2 f0 + $ 1 ! 2 f1
xn # x0 ! x1 &' #$ x0 ! x1 &' # x1 ! x0 &' #$ x1 ! x0 &'
theorem 5.11.3
# ( s""(x))
2
! among all interpolants g, dx is minimized by g = s 2 2
" x ! x1 % " x ! x0 %
+(x ! x0 ) $ f0( + (x ! x1 ) $ f1(
# x0 ! x1 '& # x1 ! x0 '&
x0

correct (clamped) b.c.

A0 B0
s1! (x0 ) = f !(x0 ), sn! (xn ) = f !(xn )
x0 x1 x0 x1

not-a-knot b.c.
s1!!!= s2!!!, sn"1
!!! = sn!!! A1 B1

periodic b.c.
s1! (x0 ) = sn! (xn ), s1!!(x0 ) = sn!!(xn )

Numerical Analysis 2, lecture 4, slide! 7 Numerical Analysis 2, lecture 4, slide! 8

 The cubic spline’s coefficients can be found The cubic spline’s coefficients can be found
by solving a tridiagonal linear system (I) by solving a tridiagonal linear system (II)

Hermite interpolation formula natural boundary conditions

si (xi!1 ) = fi!1, si (xi ) = fi , si" (xi!1 ) = gi!1, si" (xi ) = gi (i = 1, 2,…, n)
s1!!(x0 ) = 0 " # 6 f0 + 6 f1 # 4h1g0 # 2h1g1 = 0 " 2g0 + g1 = 3k1
si (x) = (1 + 2u)(1 ! u)2 fi!1 + (3 ! 2u)u 2 fi + u (1 ! u ) hi gi!1 + (u ! 1)u 2 hi gi
2 sn!!(xn ) = 0 " 6 fn#1 # 6 fn + 2hn gn#1 + 4hn gn = 0 " gn#1 + 2gn = 3kn

continuity of 2nd derivative

tridiagonal linear system
hi2 si!!(x) = (12u " 6) fi"1 + (6 " 12u) fi + (6u " 4)hi gi"1 + (6u " 2)hi gi
"2 1 % " g0 % " k1 %
hi2 si!!(xi ) = 6 fi"1 " 6 fi + 2hi gi"1 + 4hi gi $ h 2(h + h ) h ' $g ' $h k + h k '
$ 2 1 2 1 '$ 1 ' $ 2 1 1 2 '
hi2 si!!(xi"1 ) = "6 fi"1 + 6 fi " 4hi gi"1 " 2hi gi $ ! ! ! ' $ " ' = 3$ " '
$ hn 2(hn!1 + hn ) hn!1 ' $ gn!1 ' $h k + h k '
2
hi+1 !! (xi ) = "6 fi + 6 fi+1 " 4hi+1gi " 2hi+1gi+1
si+1 $ '$ ' $ n n!1 n!1 k '
$# 1 2 '& $# gn '& $# kn '&
si!!(xi ) = si+1
!! (xi ) # hi+1gi"1 + 2(hi + hi+1 )gi + hi gi+1 = 3hi+1ki + 3hi ki+1 (i = 1, 2,…, n " 1)
fi " fi"1 diagonally dominant nonsingular
where ki =
hi gaussian elimination requires O(n) operations

Numerical Analysis 2, lecture 4, slide! 9 Numerical Analysis 2, lecture 4, slide! 10

Cubic spline interpolants The boundary condition’s effect

are easily computed in Matlab diminishes with distance from boundary

x = [-0.2 0:6];
function pp = nspline(x,y) f = [1 0 0 0 0 0 0 0];
>> s = nspline(x,f); sk = spline(x,f);
>> ss = ppval(s,xx); x = x(:); y = y(:);
n = length(x)-1; sn = nspline(x,f);
>> plot(xx,ss,'r',x,f,'ko') xx = linspace(x(1),x(end));
>> plot(xx,ss-ff,’r’,x,f-f,’ko’) if n == 1
pp=mkpp(x',[diff(y)./diff(x) y(1)]); plot(xx,ppval(sk,xx),'r',xx,ppval(sn,xx),'b',x,f,'o')
else legend('not-a-knot','natural')
1 h = diff(x);
A = [[h(2:n);1;NaN],...
[2;2*(h(1:n-1)+h(2:n));2],... 1
[NaN;1;h(1:n-1)]];
T = spdiags(A,-1:1,n+1,n+1); not!a!knot
dy = diff(y); natural
k = dy./h;
0
!1 0 1 r = 3*[k(1);h(2:n).*k(1:n-1)...
+h(1:n-1).*k(2:n);k(n)]; 0
g = T\r;
0.02
d = (g(1:n) + g(2:n+1) - 2*k)./h.^2;
c = (3*k - 2*g(1:n) - g(2:n+1))./h;
b = g(1:n);
0
pp = mkpp(x',[ d c b y(1:n) ]);
end
0 6
!0.02
!1 0 1

Numerical Analysis 2, lecture 4, slide! 11 Numerical Analysis 2, lecture 4, slide! 12

The error of cubic spline interpolation what happened, what’s next
is O(h4)

theorem (p.133) • splines have degree ! k pieces and

max f (x) " s(x) ! 5
384
h 4 max f ####(x) continuous (k-1)th derivatives
x0 ! x! xn
x0 ! x! xn h = max hi

max f "(x) # s"(x) ! 9+ 3 3

216
h max
x0 ! x! xn
f """"(x) • linear spline = continuous piecewise linear
• cubic spline
x0 ! x! xn

1 + 3$ 2
max f ""(x) # s""(x) ! h max f """"(x)
x0 ! x! xn 12 x0 ! x! xn
!=
max hi
min hi
• boundary conditions: natural, correct, not-a-knot, periodic

max f """(x) # s"""(x) !

1+ $ 2
h max f """"(x) • tridiagonal equations
2 x0 ! x! xn

• splines minimize an integral

x0 ! x! xn

rules of thumb
• use small knot spacing where |f!!| is large Next: Chebyshev approximation (§9.7, 9.9)
• avoid rapidly changing knot spacing
Numerical Analysis 2, lecture 4, slide! 13 Numerical Analysis 2, lecture 4, slide! 14