Transcript Document

Chapter 15
Curve Fitting : Splines
Gab Byung Chae
Oscillation in a higher order
interpolation
The alternative : spline
Notation used to derive splines
n-1 intervals in n data points
Splines in
Linear spline
Quadratic spline
Cubic spline
st
1 ,
nd
2 ,
rd
3
order
15.2 Linear Splines
linear splines passing thru
( x1 , f 1 ),...., ( x n , f n ),
x1  x 2  ...  x n 1  x n

f 2  f1
s
(
x
)

f

( x  x1 ), x1  x  x 2
 1
1
x 2  x1




f i 1  f i

s ( x )   si ( x )  f i 
( x  x i ), x i  x  x i  1
x i 1  x i




f n  f n 1

s ( x )  f n 1 
( x  x n 1 ), x n 1  x  x n
 n 1
x n  x n 1

Example : (3,2.5),(4.5,1),(7,2.5),(9,.5) -> (5,?)
15.3 Quadratic Splines
The objective in quadratic splines is to derive
a second-order polynomial for each interval
between data points
The polynomial :
si(x) =ai + bi(x-xi)+ ci(x-xi)2
For n data points, there are n-1 intervals and,
consequently, 3(n-1) unknown constants to
evaluate.  3(n-1) equations or conditions
are required to evaluate the unknowns.
 s1 ( x )  a1  b1 ( x  x1 )  c1 ( x  x1 ) 2 , x1  x  x 2




2
s ( x )   s i ( x )  a i  bi ( x  x i )  c i ( x  x i ) , x i  x  x i 1



 s n 1 ( x )  a n 1  b n 1 ( x  x n 1 )  c n 1 ( x  x n 1 ) 2 , x n 1  x  x n

Let hi  x i 1  x i , for all i=1,,..,n-1
1. Continuity condition : the function must pass
through all the points.
si ( xi )  f i  ai  f i
--- (15.6)
for all i=1,,..,n-1 : (n-1) equations
2. The function values of adjacent polynomials
must be equal at the knots.
s i ( x i 1 )  f i 1  f i  bi ( x i 1  x i )  c i ( x i 1  x i )  f i 1
2

2
f i  bi hi  c i hi  f i 1
--- (15.8)
for all i=1,,..,n-1 : (n-1) equations
Let hi  x i 1  x i , for all i=1,,..,n-1
3. The first derivatives at the interior nodes must
be equal.
(n-2) eq.s
s i ' ( x i 1 )  s i 1 ' ( x i 1 )  bi  2 c i hi  bi 1
--- (15.9)
4. Assume that the second derivative is zero at
the first point.
s1 " ( x1 )  0  c1  0
1 equation
Example : (3,2.5),(4.5,1),(7,2.5),(9,.5) -> (5,?)
Example 15.2
Problem] Fit quadratic splines to the data in
Table 15.1. Use the results to estimate the
value at x=5.
Solution] Equation (15.8) is written for i=1
through 3 (with c1=0) to give
f1+ b1h1 = f2
f2+ b2h2 + c2h22 = f3
f3+ b3h3 + c3h32 = f4
Equation (15.9) creates 2 conditions with
c1=0 :
b1 = b2
b2+ 2c2h2 = b3
f1
f2
f3
f4
=
=
=
=
2.5
1.0
2.5
0.5
h1 = 4.5-3.0 = 1.5
h2 = 7.0-4.5 = 2.5
h3 = 9.0-7.0 = 2.0
2.5 + 1.5 b1
= 1.0
1.0
+ 2.5 b2
+ 6.25 c2
= 2.5
2.5
+ 2.0 b3
+ 4c3 = 0.5
b1 - b2
=0
b2 - b3
+ 5c2
=0
b1 = -1
b2 = -1
b3= 2.2
c2 = 0.64
c3 = -1.6
s1(x) = 2.5 –(x-3)
s2(x) = 1.0 –(x-4.5)+0.64(x-4.5)2
s3(x) = 2.5 + 2.2(x-7.0)-1.6(5-7.0)2
x=5 이므로
s2(5) = 1.0 –(5-4.5)+0.64(5-4.5)2 =0.66
15.4 Cubic Splines
 s1 ( x )  a 1  b1 ( x  x1 )  c1 ( x  x1 ) 2  d 1 ( x  x1 ) 3 , x1  x  x 2




2
3
s ( x )  s i ( x )  a i  bi ( x  x i )  c i ( x  x i )  d i ( x  x i ) , x i  x  x i 1



 s n 1 ( x )  a n 1  b n 1 ( x  x n 1 )  c n 1 ( x  x n 1 ) 2  d n 1 ( x  x n 1 ) 3 , x n 1  x  x n

si ( xi )  f i  a i  f i
a.
2
(n-1) eq.s
3
s i ( x i 1 )  f i 1  f i  bi hi  c i hi  d i hi  f i 1
2
s i ' ( x i 1 )  s i 1 ' ( x i 1 )  bi  2 c i hi  3 d i hi  bi 1
hi  x i  1  x i ,
c.
(n-1) eq.s
(n-2) eq.s
d. (n-2) eq.s
2 c n 1  6 d n 1 h n 1  0 e. 2 eq.s
s i " ( x i 1 )  s i 1 " ( x i 1 )  c i  3 d i hi  c i 1
s1 " ( x1 )  s n  1 " ( x n )  0  2 c1 
b.
Solving Cubic Splines
f. c i  3 d i h i  c i  1 ( Eq . d )  d i  ( c i  1  c i ) / 3 h i
2
3
g . f i  bi hi  c i hi  d i hi  f i  bi hi 
 bi 
f i 1  f i

hi
hi
3
hi
2
3
( 2 c i  c i  1 )  f i  1 ( Eqs . b & f )
( 2 c i  c i 1 )
2
h . b i  2 c i h i  3 d i h i  b i  h i ( c i  c i  1 )  b i  1 ( Eqs . c & f )
i . h i 1 c i 1  2 ( h i 1  h i ) c i  h i c i  1  3
j. 2 c1  0 ,
1

h
 1
0

0
0

 0
f i 1  f i
2 c n 1  6 d n 1 h n 1  2 c n  0
0
0
0
3
hi
f i  f i 1
( Eqs . g & h )
h i 1
(e, f )
0
2 ( h1  h 2 )
h2
0
0



0
0



0
0
hn  2
2 ( h n  2  h n 1 )
0
0
0
0

0   c1 



3

c2
0



0  :  


0  :  
 fn



h n 1 c n 1
3



xn
1   c n  

f3  f2
x3  x 2
0
3
f 2  f1
x 2  x1
:
 f n 1
 x n 1
:
3
0
f n 1  f n  2
x n 1  x n  2










Example 15.3
Fit cubic splines to the data . Utilize the
results to estimate the value at x =5
Solution : employ Eq.(15.27)
1

h
 1
0

0
f1
f2
f3
f4
0
0
2 ( h1  h 2 )
h2
h2
2 ( h1  h 2 )
0
0
=
=
=
=
2.5
1.0
2.5
0.5
0   c1  
0

  

0  c 2   3 ( f [ x 3 , x 2 ]  f [ x 2 , x1 ]) 
   


h 3  c 3   3 ( f [ x 4 , x 3 ]  f [ x 3 , x 2 ]) 
  

1  c4  
0
h1 = 4.5-3.0 = 1.5
h2 = 7.0-4.5 = 2.5
h3 = 9.0-7.0 = 2.0
So we have
 1

1 .5

 0

 0
0
0
8
2 .5
2 .5
9
0
0
0   c1   0 
  

0  c 2   4 .8 
   

2   c 3    4 .8 
  
1   c 4   0 
Using MATLAB the results are
c1 = 0
c2 = 0.839543726
c3 = -0.766539924
c4 = 0
Using Eq. (15.21) and
b1 = -1.419771863
b2 = -0.160456274
b3 = 0.022053232
(15.18),
d1 = 0.186565272
d2 = -0.214144487
d3 = 0.127756654
s1 ( x )  2 . 5  1 . 419771863 ( x  3 )  0 . 186565272 ( x  3 )
3
s 2 ( x )  1 . 0  0 . 160456274 ( x  4 . 5 )  0 . 839543726 ( x  4 . 5 )
 0 . 214144487 ( x  4 . 5 )
3
s 3 ( x )  2 . 5  0 . 022053232 ( x  7 . 0 )  0 . 766539924 ( x  7 . 0 )
 0 . 127756654 ( x  7 . 0 )
2
2
3
s 2 ( 5 )  1 . 0  0 . 160456274 ( 5  4 . 5 )  0 . 839543726 ( 5  4 . 5 )
 0 . 214144487 ( 5  4 . 5 )  1 . 102889734
3
2
15.4.2 End Conditions
Clamped End Condition : Specifying the first
derivatives at the first and last nodes.
Not-a-Knot End Condition : Force continuity
of the third derivative at the second and the
next-to-last knots. – same cubic functions
will apply to each of the first and last two
adjacent segments,
15.5 Piecewise interpolation in
Matlab
15.5.1 MATLAB function : spline
Syntax :
yy = spline(x, y, xx)
Example 15.4 Splines in MATLAB
Use MATLAB to fit 9 equally spaced data points
sampled from this function in the interval [-1, 1].
a) a not-a-knot spline b) a clamped spline with end
slopes of f1’=1 and f’n-1 = -4
f ( x) 
1
1  25 x
2
Solution :
a)
>> x = linspace(-1,1,9);
>> y = 1./(1+25*x.^2);
>> xx = linspace(-1,1);
>> yy = spline(x,y,xx);
>> yr = 1./(1+25*xx.^2);
>> plot(x,y,’o’,xx,yy,xx,yr,’—’)
Cubic spline of the Runge
function
f ( x) 
1
1  25 x
2
b)
>> yc =[1 y -4];
>> yyc =spline(x,yx,xx);
>> plot(x,y,’o’,xx,yyc,xx,yr,’—’)
15.5.2 MATLAB function : interp1
Implement a number of different types of
piecewise one-dimensional interpolation.
Syntax :
yi = interp1(x,y,xi, ‘method’)
yi = a vector containing the results of the
interpolation as evaluated at the points in the
vector xi , and ‘method’ = the desired method.
method
Nearest – nearest neighbor interpolation.
Linear – linear interpolation
Spline – piecewise cubic spline
interpolation(identical to the spline function)
pchip and cubic – piecewise cubic Hermite
interpolation.
The default is linear.
MATLAB’s methods
>> yy=spline(x,y,xx); cubic spilne
>> yy=interp1(x,y,xx): (piecewise) linear spline
>> yy=interp1(x,y,xx,’spline’): cubic spline
>> yy=interp1(x,y,xx,’nearest’) : nearest
neighbor method
>> yy=interp1(x,y,xx,’pchip’): piecewise cubic
Hermite interpolation
>> zz=interp2(x,y,z,xx,yy): 2 dimensional
interpolation
Example 15.5 Trade-offs Using
interpl
Linear, nearest, cubic spline, piecewise cubic
Hermite interpolation
>>
>>
>>
>>
>>
t = [ 0 20 40 56 68 80 84 96 104 110];
v = [0 20 20 38 80 80 100 100 125 125];
tt = linspace(0,110);
vl=interp1(t,v,tt); %(piecewise) linear spline
plot(t,v,’o’,tt,vl)
>> vn=interp1(t,v,tt,’nearest’); % nearest
>> plot(t,v,’o’,tt,vn)
>> vs=interp1(t,v,tt,’spline’); % cubic spline
>> plot(t,v,’o’,tt,vs)
>> vh=interp1(t,v,tt,’pchip’); % piecewise cubic
Hermite interpolation
>> plot(t,v,’o’,tt,vh)