1 - ESM Intranet Site



Chapter 4. Common Splines

A spline function is a piecewise polynomial function joined together with certain continuity conditions satisfied.

4.1. Linear Spline

For a set of data points usually termed knots: [pic], the linear spline is

[pic] (4.1.1)

The C0 conditions, [pic] and [pic] yield [pic] equations at interior points. These along with one condition at each end point give a total of [pic] equations to match the [pic] unknowns: [pic]. By applying these, we get

[pic] (4.1.2)

which results in straight lines joining neighboring points.

Clearly, [pic] is the Lagrange interpolation formula for the data set consisting of the following two points: [pic] and [pic]. Hence it is the solution approximation for linear finite elements in one dimension.

Example 4.1.1. Find the linear splines for the following data set

|i |1 |2 |3 |4 |5 |

|x |0 |5 |7 |8 |10 |

|y |0 |2 |-1 |-2 |20 |

Linear spline fit

[pic]

[pic]

[pic]

[pic]

[pic]

Figure 4.1.1 Linear Splines for data set 1.

A MATLAB code “linear_spline_examples.m” is used to generate Figure 4.1.1. Here is a brief algorithm for the code

(1) Input the data points x(i) ,y(i), i = 1,n; Plot them as symbols using plot;

(2) For i = 1, n-1, do

A = linspace(x(i),x(i+1),nint) % Divide each data interval into nint points

Call c1=linear_spline % Call the linear_spline.m function

Call c2=interp1 % Call the MATLAB function (default option: linear)

Plot c1

Plot c2

Enddo

The “linear_spline.m” function is built according to the following algorithm:

(1) Input x where the interpolation value y is needed;

Check to see if x is out of range; if yes, abort the program; if no, proceed;

(2) Find the data interval in which x resides;

(3) Calculate y according to Equation 4.1.2.

Example 4.1.2 Finding the linear splines for the following data set

|i |1 |2 |3 |4 |5 |6 |

|x |0 |1 |2 |3 |4 |5 |

|y |1 |1 |1 |-1 |-1 |-1 |

We will use the same MATLAB code for this example. The only change needed in “linear_spline_examples.m” is to replace the data set x and y with that in the table above. The splines are plotted in Figure 4.1.2.

[pic]

Figure 4.1.2 Linear Splines for data set 2.

Clearly the splines in both examples are just straight lines. We will use the two data point sets above throughout this chapter and compare the splines due to different methods.

4.2. Quadratic Spline

A quadratic spline has a quadratic function for each data interval

[pic] (4.2.1)

which is constrained to satisfy the C0 and C1 conditions.

For the C0 condition, we get

[pic] (4.2.2)

For the C1 condition, we get

[pic] (4.2.3)

From the above three requirements, there are 3n-4 constraint conditions. But the splines require a total of 3n-3 conditions, so we are one condition short and the problem is not completely defined. Usually, s(1(x1) = 0 is used for the additional condition. This results in the so-called natural quadratic spline. Certainly other conditions can be used, such as s(1(x1) = s(n-1(xn ).

To find the formula of the splines, first denote di = s(i (xi ), then since s((x) is a linear spline for the data set (xi, di, i = 1,2, (, n), we have

[pic]

[pic]

which is of course is the linear equation for the slope.

Integrating over x, we get

[pic] (4.2.4)

after solving for the constant of integration using the C0 boundary condition: [pic]. So the quadratic spline is defined once we have the [pic]’s.

Letting x = xi+1, we get

[pic]

which leads to

[pic] (4.2.5)

which is a recursive scheme. If d1 is known, then di can be derived from the above equation.

Various conditions can be imposed to obtain d1:

• For the natural spline, [pic].

• If d1 = d2, then d1 can be calculated by

[pic] (4.2.6)

• If d1 = dn, then d1 can be calculated by

[pic](4.2.7)

This condition can not be applied when n is an odd number.

The algorithm of the MATLAB code “quadratic_spline_examples.m” is as follows

(1) Input the data points x(i) ,y(i), i = 1,n; Plot them as symbols using plot;

(2) For i = 1, n-1, do

A = linspace(x(i),x(i+1),nint) % Divide each data interval into nint points

Call c1=quadratic_spline with d1 = 0; % Natural quadratic spline

Call c2=quadratic_spline with d1 = d2;

Plot c1

Plot c2

Enddo

The “quadratic_spline.m” function is built according to the following algorithm

1) Input x where the interpolation value y is needed;

Check to see if x is out of range; if yes, abort the program; if no, proceed;

2) Find the data interval where x belongs;

3) Calculate d1; Calculate di (Equation 4.2.5);

4) Calculate y according to Equation 4.2.4.

Example 4.2.1. Given the same data points in Example 4.1.1, derive the quadratic splines.

Using the natural quadratic splines, where d1 = 0, then we have

[pic]

[pic]

[pic]

[pic]

[pic]

[pic]

The results and the splines due to d1 = d2, are plotted in Figure 4.2.1. Near x = 0, the natural spline has a flatter curve due to zero slope at x = 0. The other spline has a linear slope in the first interval. Throughout all the data intervals, these two splines are different, although the differences are small for the most part.

[pic]

Figure 4.2.1. Quadratic Splines for data set 1.

Example 4.2.2. Given the same data points in Example 4.1.2, obtain the quadratic splines.

The only change needed in “quadratic_spline_examples.m” is to replace the data set x and y with the table in Example 4.1.2. The splines are plotted in Figure 4.2.2. Note that the natural splines and splines due to d1 = d2 are exactly the same. This is caused by the fact that both options produce the same result: d1 = 0.

Also, the splines in the 4th and 5th intervals have large oscillations. We would like a smooth transition in the 3rd interval ideally and the splines do have this feature. But because of the under-shoot at interval 4, there are oscillations in the 4th and 5th intervals. Generally, when the order of polynomial splines goes up, this problem can become very severe.

[pic]

Figure 4.2.2. Quadratic Splines for data set 2.

4.3. Natural Cubic Splines

Generally, a function s is called a spline of degree k on [pic] if

(1) [pic]

(2) [pic] are all continuous functions on [pic] where s(j) is the jth derivative;

(3) [pic] is a polynomial of degree [pic]on each interval [pic].

If [pic], the spline is called a cubic spline. Then the curvatures are linear in an interval and denoting values at the points [pic] as [pic], they comprise a linear spline for the data set [pic]. Therefore

[pic] (4.3.1)

Integrating over x twice, we get

[pic] (4.3.2)

where hi = xi+1 - xi and c and d are constants of integration.

Applying the C0 conditions, we get

[pic] (4.3.3)

which leads to

[pic](4.3.4)

So the cubic spline is defined once we have the ei‘s.

The values of ei can be derived from the C1 continuity condition. Differentiating the above equation, we get

[pic]

and

[pic]

where

[pic] (4.3.5)

Also

[pic]

By setting at all interior points, we have

[pic] (4.3.6)

which is a tridiagonal system of equations in terms of the unknowns [pic]’s. There are [pic] unknowns and [pic] equations. Two additional conditions are needed to determine the cubic splines. When they are [pic], the resulting splines are called the natural cubic splines. These equations can be rearranged as

[pic] (4.3.7)

where

[pic] (4.3.8)

The algorithm of the code “cubic_spline_examples.m” is very similar to “linear_spline_examples.m”, so only the algorithm for the function “cubic_spline.m” will be presented here.

1) Input x where the interpolation value y is needed; % s = y here.

2) Check to see if x is out of range; if yes, abort the program; if no, proceed;

3) Calculate hi and bi (Equation 4.3.5);

4) Calculate ui and vi (Equations 4.3.8);

5) Assign e1 = 0, then calculate ei (Equation 4.3.7)

6) Find the data interval in which x resides;

7) Calculate y according to Equation 4.3.4.

Example 4.3.1. Given the data points in Example 4.1.1, obtain the natural cubic splines.

To calculate the 2nd derivatives, first determine

h1 = 5, h2 = 2, h3 = 1, h4 = 2

b1 = 0.4, b2 = -1.5, b3 = -1, b4 = 1

Using above we can calculate ui’s and vi’s as

u2 = 14, u3 = 6, u4= 6

v2 = -11.4, v3 = 3, v4 = 12

Then from equation 4.3.7, we get

[pic]

The values of e2, e3, and e4 can be obtained by using the tridiagonal solver. These values along with e1= e5 = 0 can then be used to determine the splines.

The natural cubic splines are plotted in Figure 4.3.1. MATLAB functions “spline” and “interp1” with option set to spline produce exactly the same splines. The additional two constraints built into these two MATLAB functions are as follows

[pic]

[pic]

Figure 4.3.1. Cubic Splines for data set 1

Example 4.3.2. Given the data points in Example 4.1.2, obtain the cubic splines.

The splines are plotted in Figure 4.3.2. All of them did better than the quadratic splines. The oscillations are smaller. The natural splines here did even better than the ones due to “spline” and “interp1: spline” in the 1st and 5th intervals.

[pic]

Figure 4.3.2. Cubic Splines for data set 2.

4.4. Conclusion

Splines of increasing order are obtained by increasing the order of continuity in the matching conditions across intervals at interior points: linear splines match the data, quadratic splines match first derivatives, cubic splines match second derivatives, etc. We stop there because the concept is set and because higher order splines are rarely used in applications. But is an increase in order worth the extra mathematical effort required to generate it?

Observation of the graphs reveals that splines of order greater one oscillate through the data. A comparison of splines is shown in Figures 4.4.1 and 4.4.2. We discover that the quadratic spline, in general, oscillates more than the cubic about the linear spline which we use as a reference. So increasing the order does not always increase the oscillations. So which spline do we choose?

Selection of an appropriate spline depends on the application. If we want to interpolate the data itself, perhaps the linear is the best choice. If we have data that we know is cubic in nature, then a cubic spline may be a better choice. If we want derivatives of data, a linear spline creates jumps at data points, hence a degree of smoothness may be necessary. So again, one must make a choice based upon the application and carefully weigh the results. In design, splines allow local smoothing through data and the cubic spline is the usual choice.

How do splines compare to Hermite interpolations which also treat derivatives? It is very important to understand the following concept: splines match derivative conditions at intersections of splines, not derivatives in the data itself. In other words, splines are piecewise functions which approximate the data by filling in the voids between data: their matching conditions insure smooth connections between themselves, not the data. On the other hand, Hermite interpolations require derivatives of data, that is, the derivatives are part of the data set, e.g., position, velocity, acceleration. The only relation between spline and Hermite interpolation is that they are both polynomials.

Figure 4.4.1 Spline comparison for dataset 1

Figure 4.4.2 Spline comparison for dataset 2

-----------------------

[pic]

................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download