Cubic splines are a popular choice for curve fitting for ease of data interpolation, integration, differentiation, and they are normally very smooth. This tutorial will describe two methods of constructing joined cubic splines through known data points. The first method (Part A) is for easy understanding of the mathematical concepts involved, but is rather computationally inefficient. The 2nd method (Part B) is more involved, but is computationally efficient and reflects the approach most likely used in practice.
Consider the problem of constructing 2 joined cubic splines to fit 3 data points (x1,y1), (x2,y2), (x3,y3). This is the simplest case of cubic spline interpolation that will illustrate the methods used in more normal cases where more points are present. The key characteristics of cubic spline interpolation are:
1. The curves pass through all specified data points
2. 1st derivative (slope) continuity at interior points
3. 2nd derivative continuity at interior points
4. boundary conditions specified at the free ends
Cubic Spline Part A
We begin with the equations of the two splines:
spline #1
y = a1x3 + b1x2 + c1x + d1
y' = 3a1x2 + 2b1x + c1
y" = 6a1x + 2b1
spline #2
y = a2x3 + b2x2 + c2x + d2
y' = 3a2x2 + 2b2x + c2
y" = 6a2x + 2b2
where the 1st derivative of the function is designated by y' and the 2nd derivative by y". The 2 splines share x2 as a common point. In other words, the end of spline #1 is the beginning of spline #2. Since there are a total of 8 unknown coefficients, 8 equations are needed:
1. 0 = 6a1x1 + 2b1 ==> y1" = 0 (boundary condition at free end)
2. y1 = a1x13 + b1x12 + c1x1 + d1 ==> y1 = f(x1) in spline #1
3. y2 = a1x23 + b1x22 + c1x2 + d1 ==> y2 = f(x2) in spline #1
4. y2 = a2x23 + b2x22 + c2x2 + d2 ==> y2 = f(x2) in spline #2
5. 3a1x22 + 2b1x2 + c1 = 3a2x22 + 2b2x2 + c2 ==> 1st derivative at x2 matches in both splines
6. 6a1x2 + 2b1 = 6a2x2 + 2b2 ==> 2nd derivative at x2 matches in both splines
7. y3 = a2x33 + b2x32 + c2x3 + d2 ==> y3 = f(x3) in spline #2
8. 0 = 6a2x3 + 2b2 ==> y3" = 0 (boundary condition at free end)
The boundary conditions (equations 1 and 8) reflect what are called ‘natural’ boundary conditions, where the 2nd derivative (and curvature) is zero. Other boundary conditions can be used, for example, so called ‘parabolic runout’ (y1" = y2" and/or y3" = y2") or zero slope at either end (0 = 3a1x12 + 2b1x1 + c1 and/or 0 = 3a2x32 + 2b2x3 + c2), or any specified value of 1st or 2nd derivative.
The above 8 equations in 8 unknowns can be represented in matrix form Ax = b (Note: . means zero):
Where the solution is x = A-1b. Let’s illustrate with a specific problem: fit 2 cubic splines to the function y = x3 in the range of x = 0 to 1. Thus, x1 = 0, y1 = 0, x3 = 1, y3 = 1. We’ll pick x2 = 0.5, thus y2 = 0.125. The matrix equation to be solved is (Note: . means zero):Code:┌ ┐ ┌ ┐ ┌ ┐ │ 6x1 2 . . . . . .│ │a1│ │ 0│ │ x1^3 x1^2 x1 1 . . . .│ │b1│ │y1│ │ x2^3 x2^2 x2 1 . . . .│ │c1│ │y2│ │ . . . . x2^3 x2^2 x2 1│ │d1│ = │y2│ │3x2^2 2x2 1 . -3x2^2 -2x2 -1 .│ │a2│ │ 0│ │ 6x2 2 . . -6x2 -2 . .│ │b2│ │ 0│ │ . . . . x3^3 x3^2 x3 1│ │c2│ │y3│ │ . . . . 6x3 2 . .│ │d2│ │ 0│ └ ┘ └ ┘ └ ┘
Resulting in:Code:┌ ┐ ┌ ┐ ┌ ┐ │ . 2 . . . . . .│ │a1│ │ 0│ │ . . . 1 . . . .│ │b1│ │ 0│ │0.125 0.25 0.5 1 . . . .│ │c1│ │0.125│ │ . . . . 0.125 0.25 0.5 1│ │d1│ = │0.125│ │ 0.75 1 1 . -0.75 -1 -1 .│ │a2│ │ 0│ │ 3 2 . . -3 -2 . .│ │b2│ │ 0│ │ . . . . 1 1 1 1│ │c2│ │ 1│ │ . . . . 6 2 . .│ │d2│ │ 0│ └ ┘ └ ┘ └ ┘
As can be seen in the plot, the cubic spline interpolation doesn’t fit the function very well. Wait a minute. How can 2 cubic splines not fit a cubic polynomial very well? It should be a perfect fit, especially since it only takes 1 cubic spline to represent the function y = x3. The answer is that we forced the 2nd derivative of the spline to be zero at each free end. This works fine at x = 0 for the function y = x3 because the 2nd derivative of this function is indeed 0 at x = 0. However, it isn’t a good choice at x = 1 because the 2nd derivative of y = x3 at x = 1 is 6x = 6. If we change the boundary condition of x3 (i.e. equation 8) to reflect a value of 6 instead of 0, the fit is perfect. This highlights the importance of choosing appropriate boundary conditions for the problem at hand.Code:┌ ┐ ┌ ┐ │a1│ │ 1.5│ │b1│ │ 0│ │c1│ │-0.125│ │d1│ = │ 0│ │a2│ │ -1.5│ │b2│ │ 4.5│ │c2│ │-2.375│ │d2│ │ 0.375│ └ ┘ └ ┘
The methodology illustrated above can be extended to larger problems with more points. The problem is that the dimension of the matrix equation to be solved is 4(n-1), where n is the number of points. There is another formulation that only requires 1 unknown per point – the 2nd derivative, and the solution to the matrix equation is much simpler. See Cubic Spline Part B....
.
Featured
Writing your app was the hard part, but it’s not the last part. You still need to sell your app. Here's how to make sure your description puts your app in the best light. |
|
Reply With Quote
for you already last week 
Originally Posted by thegaunlet
