Cubic Spline Interpolation

Cubic Spline Interpolation Sky McKinley and Megan Levine Math 45: 45: Linear Linear Algebra Algebra Abstract. An introduction into the theory and application of cubic splines with accompanying Matlab m-file cspline.m

Introduction Real world numerical data is usually difficult to analyze. Any function which would effectively correlate the data would be difficult to obtain and highly unwieldy. To this end, the idea of the cubic spline was developed. Using this process, a series of unique cubic polynomials are fitted between each of the data points, with the stipulation that the curve obtained be continuous and appear smooth. These cubic splines can then be used to determine rates of change and cumulative change over an interval. In this brief introduction, we will only discuss splines which interpolate equally spaced data points, although a more robust form could encompass unequally spaced points.

Theory The fundamental idea behind cubic spline interpolation is based on the engineer’s tool used to draw smooth curves through a number of points. This spline consists of weights attached to a flat surface at the points to be connected. A flexible strip is then bent across each of these weights, resulting in a pleasingly smooth curve. The mathematical spline is similar similar in principle. The points, in this case, are numerical data. The weights are the coefficients on the cubic polynomials used to interpolate the data. These coefficients ’bend’ the line so that it passes through each of the data points without any erratic behavior or breaks in continuity.

Process The essential idea is to fit a piecewise function of the form

SÝ xÞ

=

s 1 Ý xÞ

if x 1

²

x

<

x2

s 2 Ý xÞ

if x 2

²

x

<

x3

²

x

_ s n?1 Ý xÞ if

x n?1

<

(1)

xn

where s i is a third degree polynomial defined by s i Ý xÞ

for i

=

=

a i Ý x ? x i Þ 3

+

b i Ý x ? x i Þ 2

+

c i Ý x ? x i Þ + d i

(2)

1,2,..., n ? 1.

The first and second derivatives of these n ? 1 equations are fundamental to this process, and they are

1

Titles you can't find anywhere else

Try Scribd FREE for 30 days to access over 125 million titles without ads or interruptions! Start Free Trial Cancel Anytime.







2

for i

=

s iv Ý xÞ

=

3a i Ý x ? x i Þ 2

s vvi Ý xÞ

=

6a i Ý x ? x i Þ + 2b i

+

2b i Ý x ? x i Þ + c i

(3) (4)

1,2,..., n ? 1.

The Four Properties of Cubic Splines Our spline will need to conform to the following stipulations. 1. The piecewi piecewise se function function SÝ xÞ will interpolate all data points. 2. SÝ xÞ will be continuous on the interval ß x 1 , x n à 3. S v Ý xÞ will be continuous on the interval ß x 1 , x n à 4. S vv Ý xÞ will be continuous on the interval ß x 1 , x n à Since the piecewiece function SÝ xÞ will interpolate all of the data points, we can conclude that

for i = 1,2,..., n ? 1. Since x i produce

5

SÝ x i Þ

=

yi

ß x i , x i+1 à, SÝ x i Þ

=

s i Ý x i Þ and we can use equation (2) to

y i

=

s i Ý x i Þ

y i

=

a i Ý x i

y i

=

d i

(5)

(6) xiÞ3

?

b i Ý x i

+

xiÞ2

?

+

c i Ý x i

?

x i Þ + d i

for each i = 1,2,..., n ? 1. Since the curve SÝ xÞ must be continuous across its entire interval, it can be concluded that each sub-function must join at the data points, so s i Ý x i Þ

=

s i?1 Ý x i Þ

(7)

for i = 2,3,..., n. From equation (2), s i Ý x i Þ

=

d i

and s i?1 Ý x i Þ

=

a i?1 Ý x i

?

x i?1 Þ 3

+

b i?1 Ý x i

(8) ?

x i?1 Þ 2

+

c i?1 Ý x i

?

x i?1 Þ + d i?1

so d i

for i

=

=

a i?1 Ý x i

?

2,3,..., n ? 1. Letting h

x i?1 Þ 3 =

d i

xi =

?

+

b i?1 Ý x i

?

x i?1 Þ 2

+

c i?1 Ý x i

?

x i?1 Þ + d i?1

(9)

x i?1 in equation (8), we have

a i?1 h 3

+

b i?1 h 2

+

c i?1 h + d i?1

(10)







3 However, by equation (3), s iv Ý x i Þ

=

ci

and s iv?1 Ý x i Þ

=

3a i?1 Ý x i

x i?1 Þ 2

?

+

2b i?1 Ý x i

?

x i?1 Þ + c i?1

so ci

Again, letting h

=

xi

?

=

3a i?1 Ý x i

for i = 2,3, u, n ? 1. From equation (4), s vvi Ý xÞ

i

=

=

x i?1 Þ 2

+

2b i?1 Ý x i

+

2b i?1 h + c i?1

?

x i?1 Þ + c i?1 .

(12)

x i?1 , we arrive at ci

for i

?

=

=

3a i?1 h 2

(13)

6a i Ý x ? x i Þ + 2b i , so s vvi Ý xÞ

=

6a i Ý x ? x i Þ + 2b i

s vvi Ý x i Þ

=

6a i Ý x i

s vvi Ý x i Þ

=

2b i

?

(14)

x i Þ + 2b i

2,3,..., n ? 2.

Lastly, since s vvi Ý xÞ has to be continuous across the interval, s vvi Ý x i Þ 1,2,3, `, n ? 1. This and equation (14) lead us to to the equation

and, letting h

=

x i+1

?

=

s vvi+1 Ý x i Þ for

s vvi Ý x i+1 Þ

=

6a i Ý x i+1

?

x i Þ + 2b i

(15)

s vvi+1 Ý x i+1 Þ

=

6a i Ý x i+1

?

x i Þ + 2b i

(16)

x i and using the conclusion from equations (14) and (16), s vvi+1 Ý x i+1 Þ

=

6a i Ý x i+1

2b i+1

=

6a i h + 2b i

?

x i Þ + 2b i

(17) (18)

These equations can be much simplified by substituting M i for s vvi Ý x i Þ and expressing the above equations in terms of M i and y i . This makes the determination of the weights a i, b i , c i , and d i a much easier easier task. Each b i can be represented by

and d i has already been determined to be

s vvi Ý x i Þ

=

2b i

M i

=

bi

=

2b i M i 2

(19)







4

2b i+1

=

6a i h + 2b i

6a i h

=

2b i+1

ai

=

?

(21)

2b i

2b i+1 ? 2b i 6h M M 2Ý 2i 1 Þ ? 2Ý 2 i Þ 6h M i+1 ? M i 6h +

ai

=

ai

=

and c i can be re-written as d i+1 cih

=

aih3

= ?a i h

ci

=

ci

=

ci

=

ci

=

ci

=

ci

=

ci

=

+ 3

bih2

?

+

bih2

?a i h 3 ?

bih2 h ?a i h 3 ? b i h 2 h

(22)

c i h + d i

?

d i

?

+

d i

+

d i+1

+

d i+1

?d i +

d i+1

h d ? d i+1 Ý?a i h 2 ? b i hÞ ? i h y ? y i+1 M i+1 ? M i 2 M i h + hÞ ? i ?Ý 2 6h h y i+1 ? y i M i+1 ? M i 3 M i h+ hÞ ?Ý 6 6 h y i+1 ? y i M i+1 ? M i + 3 M i Þh ?Ý 6 h y i+1 ? y i M i+1 + 2 M i Þh. ?Ý 6 h

We now have our equations for determining the weights for our n ? 1 equations

ai

=

bi

=

ci

=

d i

=

M i+1 ? M i 6h M i 2 y i+1 ? y i M i+1 + 2 M i ?Ý Þh 6 h yi

These systems can be handled more conveniently by putting them into matrix form as follows

(23)







5

c i+1

=

3a i h 2

+

2b i h + c i (24)

y ? yi M i+1 ? M i 2 M M i+1 + 2 M i Þh + 2Ý i Þh + i+1 Þh ?Ý 2 6 6h h M ? M i 2 M M + 2 M i M + 2 M i+1 3Ý i+1 Þh + 2Ý i Þh ? Ý i+1 Þh + Ý i+2 Þh 2 6 6 6h 3 M i+1 ? 3 M i 6 M i M i+1 + 2 M i M + 2 M i+1 hÝ Þ + Ý i+2 ÞÞ + ?Ý 6 6 6 6 h Ý M + 4 M + M Þ i i+1 i+2 6

3Ý

M i

for i

=

+

4 M i+1

+

M i+2

=

y i+2

?

y i+1

Ý

?

h

M i+2

y i+1 ? y i y Þ + i+2 h y i ? 2 y i+1 + y i+2 h y i ? 2 y i+1 + y i+2 h y ? 2 y i+1 + y i+2 6Ý i Þ h2

= ?Ý = = =

+

6 ? y i+1 h

1,2,3, `, n ? 1

which leads to the matrix equation M 1

1

4

1

0

` 0

0

0

0

M 2

y 1

?

2 y 2

+

y3

0

1

4

1

` 0

0

0

0

M 3

y 2

?

2 y 3

+

y4

0

0

1

4

` 0

0

0

0

M 4

y 3

?

2 y 4

+

y5

_ _ _ _ b _ _ _ _

_

=

6 h2

_

0

0

0

0

` 4

1

0

0

M n?3

y n?4

?

2 y n?3

+

y n?2

0

0

0

0

` 1

4

1

0

M n?2

y n?3

?

2 y n?2

+

y n?1

0

0

0

0

` 0

1

4

1

M n?1

y n?2

?

2 y n?1

+

yn

M n

(25) Note that this system has n ? 2 rows and n columns, and is therefore therefore under-determined. In order to generate a unique cubic spline, two other conditions must be imposed upon the system. This leads us to our next section.

Three types of Splines Natural splines This first spline type includes the stipulation that the second derivative be equal to zero at the endpoints.

2 M i+1

Þh







6

0 1

0

0

0

` 0

0

0

0

M 2

y 1

?

2 y 2

+

y3

0

1

4

1

` 0

0

0

0

M 3

y 2

?

2 y 3

+

y4

0

0

1

4

` 0

0

0

0

M 4

y 3

?

2 y 4

+

y5

_ _ _ _ b _ _ _ _

_

=

6 h2

_

0

0

0

0

` 4

1

0

0

M n?3

y n?4

?

2 y n?3

+

y n?2

0

0

0

0

` 1

4

1

0

M n?2

y n?3

?

2 y n?2

+

y n?1

0

0

0

0

` 0

0

0

1

M n?1

y n?2

2 y n?1

?

+

yn

0 (27) For reasons of convenience, the first and last columns of this matrix can be eliminated, as they correspond to the M 1 and M n values, which are both 0.

4

1

0

` 0

0

0

M 2

y 1

?

2 y 2

+

y3

1

4

1

` 0

0

0

M 3

y 2

?

2 y 3

+

y4

0

1

4

` 0

0

0

M 4

y 3

?

2 y 4

+

y5

_ _ _ b _ _ _

_

=

6 h2

_

0

0

0

` 4

1

0

M n?3

y n?4

?

2 y n?3

+

y n?2

0

0

0

` 1

4

1

M n?2

y n?3

?

2 y n?2

+

y n?1

0

0

0

` 0

1

4

M n?1

y n?2

?

2 y n?1

+

yn

(28) This results in an n ? 2 by n ? 2 matrix, which will determine the remaining solutions for M 2 through M n?1 . The spline spline is now unique. unique.









7

1 0

5

0

-5

-1 0

-1 5

-2 0 0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

Figure 1. Natural interpolating curve

Parabolic Runout Spline The parabolic spline imposes the condition that the second derivative at the endpoints, M 1 and M n , be equal to M 2 and M n?1 respectively. M 1

=

M 2

M n

=

M n?1

(29)

The result of this condition is the curve becomes a parabolic curve at the endpoint. This type of cubic spline is useful for periodic and exponential data. The matrix equation for this type of spline is







8

1 0

0

-1 0

-2 0

-3 0

-4 0

-5 0

-6 0 0

1

2

3

4

5

6

Figure 2. Parabolic Runout curve Note that the endpoint behavior is a bit more extreme than with the natural spline option.

Cubic Runout Spline This last type of spline has the most extreme endpoint behavior. It assigns M 1 to be 2 M 2 ? M 3 and M n to be 2 M n?1 ? M n?2 . This causes the curve to degrade to a single single cubic curve over the last two intervals, rather than two separate functions. The matrix equation for this type is







9

2 0

0

-2 0

-4 0

-6 0

-8 0

-1 0 0

-1 2 0

-1 4 0

-1 6 0 0

1

2

3

4

5

6

Figure 3. Cubic Runout Runout curve curve Note the pronounced curvature at the endpoints. Keep in mind that there are many other types of interpolating spline curves, such as the Periodic Spline and the Clamped Spline. The three discussed in this work are simply the ones which we have chosen to examine; they are not intrinsically superior to, or more widely used than these other types of splines.

Curve Fitting with Splines The main application of cubic spline interpolation techniques is, of course, curve fitting. To this end, the consistency and efficiency of the spline as a data correlation tool will be demonstrated.







10

1

0. 8

0. 6

0. 4

0. 2

0

-0 .2

-0 .4

-0 .6

-0 .8

-1 -4

-3

-2

-1

0

Figure 4. The graph of y

1

=

2

3

4

sinÝ xÞ.

This next figure was generated by connecting six data points along the above line with a cubic spline function. 1

0 .8

0 .6

0 .4

0 .2

0

-0 .2







11

1

0 .8

0 .6

0 .4

0 .2

0

-0 .2

-0 .4

-0 .6

-0 .8

-1 -3

-2

-1

0

1

2

3

Figure 6. Superimposition of a spline curve on a sine function. Clearly, the cubic curve closely imitates imitates the sine curve. It has no extreme behavior between data points, and it effectively correlates the points. This characteristic also works with more erratic functions. Take for instance the function below.









12

1

0.8

0.6

0.4

0.2

0

-0.2

-0.4

-0.6

-0.8 -6

-4

-2

0

2

4

6

3

Figure 7. A cubic spline approximation of ßsinÝ xÞ + cosÝ xÞà 4 .







13

1

0. 9

0. 8

0. 7

0. 6

0. 5

0. 4

0. 3

0. 2

0. 1

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figure 8. Random data points Clearly there is no relation between these data points. A spline, however, can interpolate all 100 points without the drastic behavior that the necessary 99th degree polynomial would exhibit.







14

1

0. 9

0. 8

0. 7

0. 6

0. 5

0. 4

0. 3

0. 2

0. 1

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figure 10. Approximating cubic curve The above figure demonstrates the general trends in the data (of which there are none) without necessarily connecting all data points. Its curvature is much less severe than that of fig 9.







15 Rorres, Chris and Howard Anton. Applications of Linear Algebra 3ed . 1984 New York: John Wiley and Sons. Van Loan, Charles F. Introduction to Scientific Computing . 1997 New Jersey: Prentice Hall

Cubic Spline Interpolation

Recommend Documents