Cubic Spline Tutorial

1 Attachment(s)

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 (x₁,y₁), (x₂,y₂), (x₃,y₃). 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 = a₁x³ + b₁x² + c₁x + d₁
y' = 3a₁x² + 2b₁x + c₁
y" = 6a₁x + 2b₁

spline #2
y = a₂x³ + b₂x² + c₂x + d₂
y' = 3a₂x² + 2b₂x + c₂
y" = 6a₂x + 2b₂

where the 1st derivative of the function is designated by y' and the 2nd derivative by y". The 2 splines share x₂ 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 = 6a₁x₁ + 2b₁ ==> y₁" = 0 (boundary condition at free end)
2. y₁ = a₁x₁³ + b₁x₁² + c₁x₁ + d₁ ==> y₁ = f(x₁) in spline #1
3. y₂ = a₁x₂³ + b₁x₂² + c₁x₂ + d₁ ==> y₂ = f(x₂) in spline #1
4. y₂ = a₂x₂³ + b₂x₂² + c₂x₂ + d₂ ==> y₂ = f(x₂) in spline #2
5. 3a₁x₂² + 2b₁x₂ + c₁ = 3a₂x₂² + 2b₂x₂ + c₂ ==> 1st derivative at x₂ matches in both splines
6. 6a₁x₂ + 2b₁ = 6a₂x₂ + 2b₂ ==> 2nd derivative at x₂ matches in both splines
7. y₃ = a₂x₃³ + b₂x₃² + c₂x₃ + d₂ ==> y₃ = f(x₃) in spline #2
8. 0 = 6a₂x₃ + 2b₂ ==> y₃" = 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’ (y₁" = y₂" and/or y₃" = y₂") or zero slope at either end (0 = 3a₁x₁² + 2b₁x₁ + c₁ and/or 0 = 3a₂x₃² + 2b₂x₃ + c₂), 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):

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│ └ ┘ └ ┘ └ ┘

Where the solution is x = A^-1b. Let’s illustrate with a specific problem: fit 2 cubic splines to the function y = x³ in the range of x = 0 to 1. Thus, x₁ = 0, y₁ = 0, x₃ = 1, y₃ = 1. We’ll pick x₂ = 0.5, thus y₂ = 0.125. The matrix equation to be solved is (Note: . means zero):

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│ └ ┘ └ ┘ └ ┘

Resulting in:

Code:

┌ ┐ ┌ ┐ │a1│ │ 1.5│ │b1│ │ 0│ │c1│ │-0.125│ │d1│ = │ 0│ │a2│ │ -1.5│ │b2│ │ 4.5│ │c2│ │-2.375│ │d2│ │ 0.375│ └ ┘ └ ┘

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 = x³. 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 = x³ 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 = x³ at x = 1 is 6x = 6. If we change the boundary condition of x₃ (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.

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....
.

Re: Cubic Spline Tutorial

Cubic Spline Part B

We begin with a slightly different form of the equations for the 2 cubic splines:

spline #1
y = a₁(x-x₁)³ + b₁(x-x₁)² + c₁(x-x₁) + d₁
y' = 3a₁(x-x₁)² + 2b₁(x-x₁) + c₁
y" = 6a₁(x-x₁) + 2b₁

spline #2
y = a₂(x-x₂)³ + b₂(x-x₂)² + c₂(x-x₂) + d₂
y' = 3a₂(x-x₂)² + 2b₂(x-x₂) + c₂
y" = 6a₂(x-x₂) + 2b₂

For now, we’ll focus on spline #1. We start with the 2nd derivative. Imposing the compatibility constraints that y" = y₁" at x = x₁ and y" = y₂" at x = x₂, and calling x₂-x₁ = h₁:

y₁" = 6a₁(x₁-x₁) + 2b₁ = 0 + 2b₁ = 2b₁ ==> b₁ = y₁"/2
y₂" = 6a₁(x₂-x₁) + 2b₁ = 6a₁h₁ + y₁" ==> a₁ = (y₂"-y₁")/6h₁

This results in the following equation for the 2nd derivative:

y" = (x-x₁)(y₂"-y₁")/(x₂-x₁) + y₁"

which can be verified to be correct (i.e. y" = y₁" at x = x₁ and y" = y₂" at x = x₂)

Next, we apply the conditions that the spline pass though the points, in other words y₁ = f(x₁) and y₂ = f(x₂):

y₁ = 0 + 0 + 0 + d₁
d₁ = y₁

y₂ = (x₂-x₁)³(y₂"-y₁")/6h₁ + y₁"(x₂-x₁)²/2 + c₁(x₂-x₁) + y₁
y₂ = h₁³(y₂"-y₁")/6h₁ + y₁"h₁²/2 + c₁h₁ + y₁
y₂ = h₁²(y₂"-y₁")/6 + y₁"h₁²/2 + c₁h₁ + y₁
y₂-y₁ = y₂"h₁²/6 - y₁"h₁²/6 + y₁"h₁²/2 + c₁h₁
(y₂-y₁)/h₁ = y₂"h₁/6 - y₁"h₁/6 + y₁"h₁/2 + c₁
(y₂-y₁)/h₁ = y₂"h₁/6 - y₁"h₁/6 + 3y₁"h₁/6 + c₁
(y₂-y₁)/h₁ = y₂"h₁/6 + y₁"h₁/3 + c₁
c₁ = (y₂-y₁)/h₁ – y₂"h₁/6 - y₁"h₁/3

Finally, we impose the compatibility condition that y₂' in spline #1 must equal y₂' in spline #2:

3a₁(x₂-x₁)² + 2b₁(x₂-x₁) + c₁ = 3a₂(x₂-x₂)² + 2b₂(x₂-x₂) + c₂
3a₁h₁² + 2b₁h₁ + c₁ = c₂
h₁(y₂"-y₁")/2 + y₁"h₁ + (y₂-y₁)/h₁ - y₂"h₁/6 - y₁"h₁/3 = (y₃-y₂)/h₂ – y₃"h₂/6 – y₂"h₂/3
h₁(y₂"-y₁")/2 + y₁"h₁ - y₂"h₁/6 - y₁"h₁/3 + y₃"h₂/6 + y₂"h₂/3 = (y₃-y₂)/h₂ - (y₂-y₁)/h₁
3h₁(y₂"-y₁") + 6y₁"h₁ - y₂"h₁ - 2y₁"h₁ + y₃"h₂ + 2y₂"h₂ = 6(y₃-y₂)/h₂ - 6(y₂-y₁)/h₁
3h₁y₂" – 3h₁y₁" + 6y₁"h₁ - y₂"h₁ - 2y₁"h₁ + y₃"h₂ + 2y₂"h₂ = 6(y₃-y₂)/h₂ - 6(y₂-y₁)/h₁
y₁"(6h₁ – 3h₁ – 2h₁) + y₂"(2h₁ + 2h₂) + y₃"h₂ = 6(y₃-y₂)/h₂ - 6(y₂-y₁)/h₁

y₁"h₁ + y₂"(2h₁ + 2h₂) + y₃"h₂ = 6(y₃-y₂)/h₂ - 6(y₂-y₁)/h₁

Generalizing, this equation results in a tri-diagonal set of linear equations, the unknowns being the 2nd derivatives of the interior points. The end points require special consideration as before.

Code:

┌ ┐ ┌ ┐ ┌ ┐ │ ? │ │ y1" │ │ ? │ │h1 2(h1+h2) h2 │ │ y2" │ │ (y3 - y2)/h2 - (y2 - y1)/h1 │ │ h2 2(h2+h3) h3 │ │ y3" │ │ (y4 - y3)/h3 - (y3 - y2)/h2 │ │ ... │ │ ... │ = 6 │ ... │ │ ... │ │ ... │ │ ... │ │ hn-2 2(hn-2+hn-1) hn-1 │ │yn-1"│ │(yn-yn-1)/hn-1 - (yn-1-yn-2)/hn-2│ │ ? │ │ yn"│ │ ? │ └ ┘ └ ┘ └ ┘

Even more simplification can be done if equal point spacing is used (i.e. h₁ = h₂ = ... = h_n-1 = h):

Code:

┌ ┐ ┌ ┐ ┌ ┐ │? │ │y1" │ │ ? │ │1 4 1 │ │y2" │ │ y3 - 2y2 + y1 │ │ 1 4 1 │ │y3" │ │ y4 - 2y3 + y2 │ │ ... │ │... │ = 6/(h*h) │ ... │ │ 1 4 1 │ │... │ │ ... │ │ 1 4 1 │ │yn-1"│ │yn - 2yn-1 + yn-2│ │ ? │ │yn" │ │ ? │ └ ┘ └ ┘ └ ┘

1 Attachment(s)

Re: Cubic Spline Tutorial

Cubic Spline Part B continued...

Let’s solve the same problem using the same points (x = [0,0.5, 1]^T, y = [0, 0.125, 1]^T) and boundary conditions (y₁" = y₃" = 0 that we did using the 1st method. Because the only unknowns are the 2nd derivative at each point, we have a 3 x 3 matrix to solve (vs 8 x 8 using the first method). Also since h₁ = h₂ = 0.5, we can use the simplified version:

Code:

┌ ┐ ┌ ┐ ┌ ┐ ┌ ┐ ┌ ┐ ┌ ┐ │1 0 0 │ │y1"│ │ 0 │ │ 0 │ │0 │ │ 0│ │1 4 1 │ │y2"│ = 6/(h*h)│y3-2y2+y1│ = 6/0.25 │1-2(0.125)+0│ = 24│0.75│ = │18│ │0 0 1 │ │y3"│ │ 0 │ │ 0 │ │0 │ │ 0│ └ ┘ └ ┘ └ ┘ └ ┘ └ ┘ └ ┘

The solution is [y₁", y₂", y₃"]^T = [0, 4.5, 0]^T. From this we can calculate the coefficients of the cubic spline segments:

a₁ = (y₂"-y₁")/6h₁ = 1.5
b₁ = y₁"/2 = 0
c₁ = (y₂-y₁)/h₁ – y₂"h₁/6 – y₁"h₁/3 = -0.125
d₁ = y₁ = 0

a₂ = (y₃"-y₂")/6h₂ = -1.5
b₂ = y₂"/2 = 2.25
c₂ = (y₃-y₂)/h₂ – y₃"h₂/6 – y₂"h₂/3 = 0
d₂ = y₂ = 0.125

The resulting plot is identical to that in the first example. In summary, using this formulation, the dimension of the matrix problem is n instead of 4(n-1). Also, tri-diagonal matrices require fewer computational resources. To be sure, additional computation is needed to find the spline coefficients, which are the direct result of the matrix solution in the first method, but the total computation is far less than the first method.

As a final illustration, we will show how to enforce a slope at either end. Recall the equation of the 1st derivative:

At x = x₁:

y₁' = 3a₁(x₁-x₁)² + 2b₁(x₁-x₁) + c₁ = 0 + 0 + c₁ = c₁ = (y₂-y₁)/h₁ – y₂"h₁/6 – y₁"h₁/3
y₁"(2h₁) + y₂"(h₁) = 6[(y₂-y₁)/h₁ - y₁']

At x = x₃:

y₃' = 3a₂(x₃-x₂)² + 2b₂(x₃-x₂) + c₂ = 3a₂h₂² + 2b₂h₂ + c₂
y₃' = 3h₂²(y₃"-y₂")/6h₂ + 2h₂y₂"/2 + (y₃-y₂)/h₂ – y₃"h₂/6 – y₂"h₂/3
y₃' = 3h₂(y₃"-y₂")/6 + h₂y₂" + (y₃-y₂)/h₂ – y₃"h₂/6 – y₂"h₂/3
h₂y₃"/2 - h₂y₂"/2 + h₂y₂" – y₃"h₂/6 – y₂"h₂/3 = y₃' - (y₃-y₂)/h₂
3h₂y₃" - 3h₂y₂" + 6h₂y₂" – y₃"h₂ – 2y₂"h₂ = 6(y₃' - (y₃-y₂)/h₂)
y₃"(3h₂ – h₂) - y₂"(6h₂ – 3h₂ – 2h₂) = 6(y₃' - (y₃-y₂)/h₂)
y₂"(h₂) + y₃"(2h₂) = 6[y₃' - (y₃-y₂)/h₂]

If we force the slope to be zero at both ends, the matrix equation would be (again, using the simplified version):

Code:

┌ ┐ ┌ ┐ ┌ ┐ ┌ ┐ ┌ ┐ ┌ ┐ │2 1 0 │ │y1"│ │ y2-y1 │ │ 0.125-0 │ │ 0.125│ │ 3│ │1 4 1 │ │y2"│ = 6/(h*h)│y3-2y2+y1│ = 6/0.25 │1-2(0.125)+0│ = 24│ 0.75 │ = │ 18│ │0 1 2 │ │y3"│ │ -(y3-y2)│ │ -(1-0.125)│ │-0.875│ │-21│ └ ┘ └ ┘ └ ┘ └ ┘ └ ┘ └ ┘

Conclusion:

We have demonstrated two different methods of formulating cubic splines to interpolate a given set of points. Using the first method, the solution of the matrix equation directly yields the coefficients of the spline segments, but has the disadvantage of a matrix dimension of 4(n-1). The second method is much more computationally efficient (matrix dimension n) but requires some calculation to determine the spline segment coefficients. We also demonstrated how to implement various free end boundary conditions.

To extend the discussion to parametric cubic splines, see:http://www.vbforums.com/showthread.php?t=481546
.
:)
.

Re: Cubic Spline Tutorial

VBAhack thanks for the excellent posts and job well done!

Unfortunately I am not allowed by the administrator to rate your post as I have already done it the last time.

Re: Cubic Spline Tutorial

Thanks for this

Re: Cubic Spline Tutorial

Nice job! I too would have rated it... but I can't...

I always wondered how splines were actually calculated. I never thought of the 2nd derivative -> 0 boundary conditions, or the fact that cubic fits would yield such nice results. Thanks!

{Space for edit after I get the chance to read it over more carefully}

Re: Cubic Spline Tutorial

Rassis and jemidiah, don't worry, I rated him up :thumb: for you already last week ;)

Re: Cubic Spline Tutorial

Thanks all for the positive comments. :thumb: If that wasn't enough punishment, I'm finalizing a follow-on dealing with parametric cubic splines. Stay tuned...

Re: Cubic Spline Tutorial

This is one of the best cubic spline interpolation tutorials available online! Thanks to the simple representation of Cubic Spline Part A in the form of a matrix, I was able to use excel to help me with interpolating my data. Thanks alot!

Is there a way to edit the equation to set the third derivative of both splines to equal each other as I've read of "not-a-knot" boundary end conditions?

Re: Cubic Spline Tutorial

Glad it was helpful. The short answer is no. A cubic polynomial has 4 parameters to specify which means only 4 'things' can be fixed. By fixing both end points and the 1st & 2nd derivatives at the join point, that's 4 things. Can't do any more. The only way to match the 3rd derivative is to use a quartic polynomial instead of a cubic or not match the 1st or 2nd derivatives.

Re: Cubic Spline Tutorial

Thank you so much for the prompt reply! So it's a "quartic spline interpolation"? Many say MATLAB uses the "not-a-knot" boundary end condition for cubic spline interpolation, is it erroneous? So if I match the 4 basic conditions together with the 3rd derivative of the quartic polynomial, it'll have a "better" fit than a "natural cubic spline"?

Re: Cubic Spline Tutorial

Great article, thanks.
If I want to use this in 3D (like in an image) can that be done with a 2D algoritm?

Re: Cubic Spline Tutorial

Quote:

Originally Posted by thegaunlet

Thank you so much for the prompt reply! So it's a "quartic spline interpolation"? Many say MATLAB uses the "not-a-knot" boundary end condition for cubic spline interpolation, is it erroneous? So if I match the 4 basic conditions together with the 3rd derivative of the quartic polynomial, it'll have a "better" fit than a "natural cubic spline"?

I'd never heard of "not-a-knot" boundary end condition, so I did a little research: http://www.mathworks.com/products/sp.../csapidem.html

From what I can tell, "not-a-knot" is just one of several ways to specify the end conditions of a cubic spline. If you look at the details of cubic spline formulation, the end conditions at both ends of the spline must be specified in order to solve the problem. I view them as constraints. Thus each end spline needs an additional contraint (in the form of a derivative) to be specified.

There are several options. One is what's called "natural" boundary conditions, which is zero second derivative at the end points. Another option is to specify the first derivative at the end points. The "not-a-knot" condition just means you are forcing the 3rd derivative at the 1st and last interior point to be continuous. It's just another option for end conditions of a cubic spline. No magic.

Is "not-a-knot" boundary end condition better than other choices? I doubt there is a universal answer that can cover all possible situations where cubic spline interpolation is used. Seems to me it depends entirely on what you are trying to fit. It could very well be that it is better for a certain class of problems, I don't know.

Re: Cubic Spline Tutorial

Thanks alot VBAhack (= You've been a great help! I'm interpolating biomechanical data, mainly distance, time and angular parameters of gait, hence will require interpolation to ensure the intervals are the same to compare the patterns between individuals/gender. I would have never been able to understand any of this so quickly without your tutorial. Kudos to you (=

Re: Cubic Spline Tutorial

There are infinitely many ways to interpolate. Some of those (also infinitely many) are sane. A small fraction of those (finite, but decently large) are regularly used by humans. Arbitrarily choosing to use cubic spline interpolation out of all the possible choices might cause some statistical problems--why is that interpolation valid? At the least, you should try to justify that assumption somehow. One way is by taking highly accurate data (if that's possible) and seeing how close the interpolation of coarse data matches with fine data.

Re: Cubic Spline Tutorial

Wow, thanks.

I have been trying to do some 4D polynomial interpolation without much success. Kindly help --- how does one stretch the cubic spline to 4 dimensions ? By 4 dimensions, I mean 4 variables. Thanks.

Re: Cubic Spline Tutorial

I'm sure there are a ton of ways you could use the above 2D cubic splines in higher dimensions. Check out VBAhack's parametric cubic spline tutorial for the start of one. Instead of just

x = f₁(s) = a_x(s-s₀)³ + b_x(s-s₀)² + c_x(s-s₀) + d_x
y = f₂(s) = a_y(s-s₀)³ + b_y(s-s₀)² + c_y(s-s₀) + d_y

you would have

w = f₁(s) = a_w(s-s₀)³ + b_w(s-s₀)² + c_w(s-s₀) + d_w
x = f₂(s) = a_x(s-s₀)³ + b_x(s-s₀)² + c_x(s-s₀) + d_x
y = f₃(s) = a_y(s-s₀)³ + b_y(s-s₀)² + c_y(s-s₀) + d_y
z = f₄(s) = a_z(s-s₀)³ + b_z(s-s₀)² + c_z(s-s₀) + d_z

The rest of the algebra is similar, but with more columns in your matrices and more boundary conditions.

Re: Cubic Spline Tutorial

Quote:

Originally Posted by jemidiah

I'm sure there are a ton of ways you could use the above 2D cubic splines in higher dimensions. Check out VBAhack's parametric cubic spline tutorial for the start of one. Instead of just

x = f₁(s) = a_x(s-s₀)³ + b_x(s-s₀)² + c_x(s-s₀) + d_x
y = f₂(s) = a_y(s-s₀)³ + b_y(s-s₀)² + c_y(s-s₀) + d_y

you would have

w = f₁(s) = a_w(s-s₀)³ + b_w(s-s₀)² + c_w(s-s₀) + d_w
x = f₂(s) = a_x(s-s₀)³ + b_x(s-s₀)² + c_x(s-s₀) + d_x
y = f₃(s) = a_y(s-s₀)³ + b_y(s-s₀)² + c_y(s-s₀) + d_y
z = f₄(s) = a_z(s-s₀)³ + b_z(s-s₀)² + c_z(s-s₀) + d_z

The rest of the algebra is similar, but with more columns in your matrices and more boundary conditions.

Urm ... are you sure that would work in 4 dimensions ? Looks to easy .....

Re: Cubic Spline Tutorial

Something looking too easy is no real reason to think it won't work. The obvious generalization is often the correct one. If you pick the right boundary conditions and spacing, it'll "work" in that (f1, f2, f3, f4) can be determined uniquely to pass through the points you give it. That function will be continuous by definition. I dunno what other properties you want, since you've been vague.

Re: Cubic Spline Tutorial

Quote:

Originally Posted by jemidiah

Something looking too easy is no real reason to think it won't work. The obvious generalization is often the correct one. If you pick the right boundary conditions and spacing, it'll "work" in that (f1, f2, f3, f4) can be determined uniquely to pass through the points you give it. That function will be continuous by definition. I dunno what other properties you want, since you've been vague.

Ok. To be precise, I have a value, call it gamma, that is a function of 4 variables, call them (w, x, y and z), i.e. gamma =fn(w,x,y,z). Each of these 4 variables are independent and have their own ranges. E.g, w ranges from 0 to 0.1, x from 200 to 3000, y from 0.2 to 200 and z from 0 to 0.08. (The ranges have their own independent steps).
Thus, one can form 4 spline curves for each variable and each variable will have it's own 4 points for interpolation. So, a user will give inputs of w, x, y and z and the user will be seeking the value of gamma. Now, can the parametric method be used for this interpolation in 4 variables ?
Perhaps another way of putting the question is, "How does one combine the 4 independent curves" ?
Thanks.

Re: Cubic Spline Tutorial

Quote:

Originally Posted by frankly747

Perhaps another way of putting the question is, "How does one combine the 4 independent curves" ?
Thanks.

Oh, you want to fit a surface and not a line. That's a different question than the one I was thinking--I was imagining you wanted a cubic spline (which is by definition a curve, which is by definition 1-dimensional) that simply lived in 4D space, which passed through each of your points in turn. I'll have to think about surface fitting, though you might want to look it up elsewhere, since I'm sure people have thought about it a lot.

Re: Cubic Spline Tutorial

Quote:

Originally Posted by frankly747

Thus, one can form 4 spline curves for each variable and each variable will have it's own 4 points for interpolation.

How are you doing this? I can imagine interpolating along some coordinate lines, but you'd need a bunch of sample data points first. In any case, sure you can use cubic splines to help in surface fitting in 4 variables, and in a ton of different ways. Should you? It depends entirely on your application, which again you haven't been specific about.

You might be able to get by just fine with linear interpolations. To be somewhat more specific, suppose you know gamma for any integer values from 1 to 5 in each variable [I know this isn't quite your situation, but it's clearly analogous by scaling and shifting your ranges]. To compute, say, gamma(1.1, 2.0, 3.5, 4.9) you could take several weighted averages:

Code:

gamma(1.1, 2.0, 3.5, 4.9) = [0.9*gamma(1, 2.0, 3.5, 4.9) + 0.1*gamma(2, 2.0, 3.5, 4.9)] gamma(1, 2.0, 3.5, 4.9) = [1*gamma(1, 2, 3.5, 4.9) + 0*gamma(1, 3, 3.5, 4.9)] gamma(2, 2.0, 3.5, 4.9) = [1*gamma(2, 2, 3.5, 4.9) + 0*gamma(2, 3, 3.5, 4.9)] gamma(1, 2, 3.5, 4.9) = [0.5*gamma(1, 2, 3, 4.9) + 0.5*gamma(1, 2, 4, 4.9)] gamma(1, 3, 3.5, 4.9) = [0.5*gamma(1, 3, 3, 4.9) + 0.5*gamma(1, 3, 4, 4.9)] gamma(2, 2, 3.5, 4.9) = [0.5*gamma(2, 2, 3, 4.9) + 0.5*gamma(2, 2, 4, 4.9)] gamma(2, 3, 3.5, 4.9) = [0.5*gamma(2, 3, 3, 4.9) + 0.5*gamma(2, 3, 4, 4.9)] gamma(1, 2, 3, 4.9) = [0.1*gamma(1, 2, 3, 4) + 0.9*gamma(1, 2, 3, 5)] gamma(1, 2, 4, 4.9) = [0.1*gamma(1, 2, 4, 4) + 0.9*gamma(1, 2, 4, 5)] gamma(1, 3, 3, 4.9) = [0.1*gamma(1, 3, 3, 4) + 0.9*gamma(1, 3, 3, 5)] gamma(1, 3, 4, 4.9) = [0.1*gamma(1, 3, 4, 4) + 0.9*gamma(1, 3, 4, 5)] gamma(2, 2, 3, 4.9) = [0.1*gamma(2, 2, 3, 4) + 0.9*gamma(2, 2, 3, 5)] gamma(2, 2, 4, 4.9) = [0.1*gamma(2, 2, 4, 4) + 0.9*gamma(2, 2, 4, 5)] gamma(2, 3, 3, 4.9) = [0.1*gamma(2, 3, 3, 4) + 0.9*gamma(2, 3, 3, 5)] gamma(2, 3, 4, 4.9) = [0.1*gamma(2, 3, 4, 4) + 0.9*gamma(2, 3, 4, 5)]

or in all

Code:

gamma(1.1, 2.0, 3.5, 4.9) = sum over i=1,2; j=2,3; k=3,4; m=4,5 of interp(1.1, i)*interp(2.0, j)*interp(3.5, k)*interp(4.9, m)*gamma(i, j, k, m) interp(x, x0) = 1 - abs(x - x0)

If this discussion gets more involved, though, we should start another thread and not hijack this one.

Re: Cubic Spline Tutorial

Quote:

Originally Posted by jemidiah

How are you doing this? I can imagine interpolating along some coordinate lines, but you'd need a bunch of sample data points first. In any case, sure you can use cubic splines to help in surface fitting in 4 variables, and in a ton of different ways. Should you? It depends entirely on your application, which again you haven't been specific about.

You might be able to get by just fine with linear interpolations. To be somewhat more specific, suppose you know gamma for any integer values from 1 to 5 in each variable [I know this isn't quite your situation, but it's clearly analogous by scaling and shifting your ranges]. To compute, say, gamma(1.1, 2.0, 3.5, 4.9) you could take several weighted averages:

Code:

gamma(1.1, 2.0, 3.5, 4.9) = [0.9*gamma(1, 2.0, 3.5, 4.9) + 0.1*gamma(2, 2.0, 3.5, 4.9)] gamma(1, 2.0, 3.5, 4.9) = [1*gamma(1, 2, 3.5, 4.9) + 0*gamma(1, 3, 3.5, 4.9)] gamma(2, 2.0, 3.5, 4.9) = [1*gamma(2, 2, 3.5, 4.9) + 0*gamma(2, 3, 3.5, 4.9)] gamma(1, 2, 3.5, 4.9) = [0.5*gamma(1, 2, 3, 4.9) + 0.5*gamma(1, 2, 4, 4.9)] gamma(1, 3, 3.5, 4.9) = [0.5*gamma(1, 3, 3, 4.9) + 0.5*gamma(1, 3, 4, 4.9)] gamma(2, 2, 3.5, 4.9) = [0.5*gamma(2, 2, 3, 4.9) + 0.5*gamma(2, 2, 4, 4.9)] gamma(2, 3, 3.5, 4.9) = [0.5*gamma(2, 3, 3, 4.9) + 0.5*gamma(2, 3, 4, 4.9)] gamma(1, 2, 3, 4.9) = [0.1*gamma(1, 2, 3, 4) + 0.9*gamma(1, 2, 3, 5)] gamma(1, 2, 4, 4.9) = [0.1*gamma(1, 2, 4, 4) + 0.9*gamma(1, 2, 4, 5)] gamma(1, 3, 3, 4.9) = [0.1*gamma(1, 3, 3, 4) + 0.9*gamma(1, 3, 3, 5)] gamma(1, 3, 4, 4.9) = [0.1*gamma(1, 3, 4, 4) + 0.9*gamma(1, 3, 4, 5)] gamma(2, 2, 3, 4.9) = [0.1*gamma(2, 2, 3, 4) + 0.9*gamma(2, 2, 3, 5)] gamma(2, 2, 4, 4.9) = [0.1*gamma(2, 2, 4, 4) + 0.9*gamma(2, 2, 4, 5)] gamma(2, 3, 3, 4.9) = [0.1*gamma(2, 3, 3, 4) + 0.9*gamma(2, 3, 3, 5)] gamma(2, 3, 4, 4.9) = [0.1*gamma(2, 3, 4, 4) + 0.9*gamma(2, 3, 4, 5)]

or in all

Code:

gamma(1.1, 2.0, 3.5, 4.9) = sum over i=1,2; j=2,3; k=3,4; m=4,5 of interp(1.1, i)*interp(2.0, j)*interp(3.5, k)*interp(4.9, m)*gamma(i, j, k, m) interp(x, x0) = 1 - abs(x - x0)

If this discussion gets more involved, though, we should start another thread and not hijack this one.

Hi all.
Thanks for the replies. To the gentleman who suggested linear interpolation --- I have already done that and it works but the accuracy is not good, that's why I need to do non-linear interpolation for four variables.

Now, I have since learnt that one of the "easiest" ways to do it is by use of tensor product constructs.
Now, my question is, given, say, 4 univariate cubic polynomials (each of the 4 is interpolating its respective variable ), how does one multiply all four together ? (That is the definition of a tensor product construct in this context).
Thus, given the cubic polynomials,
var1 = aw^3 + bw^2 + bw + d
var2 = ex^3 + fx^2 + gx + h
var3 = iy^3 + jy^2 + ky + l
var4 = mz^3 + nz^2 + oz + p

F(w,x,y,z) = tensor product construct. How does one do that ?
Now please, no references to wikipedia articles -- or any other articles for that matter. I have read and re-read lots of articles online and at this point, they wont help --- what I need now is for someone who knows how to do a tensor product construct of the 4 univariate cases above to produce a multivariate case in 4 dimensions.
Thanks.

Frank.

Re: Cubic Spline Tutorial

Quote:

Originally Posted by frankly747

...what I need now is for someone who knows how to do a tensor product construct of the 4 univariate cases above to produce a multivariate case in 4 dimensions.

That would not be me (though the tensor product is probably just a convenient place to put a lot of indexes). However, I can generalize the linear method above to cubic splines, assuming you can evaluate gamma at lattice points, say (n1, n2, n3, n4) for integer ni's. The general idea is to start by interpolating along each coordinate axis--that is, allow only one of the four variables to vary. This would generate a lot of interpolations--for a 4x4x4x4 lattice, it would generate 4^4 = 256 interpolations, however only a few of these are used. Combine these interpolations by successively creating new interpolations based on intermediate points in these coordinate axis interpolations. After several new interpolations, you'll have your point. I can be more specific if it would help. I would bet that the method I just described is equivalent to an existing method, but I'm not into computer graphics enough to know specifically.

Re: Cubic Spline Tutorial

Quote:

Originally Posted by jemidiah

That would not be me (though the tensor product is probably just a convenient place to put a lot of indexes). However, I can generalize the linear method above to cubic splines, assuming you can evaluate gamma at lattice points, say (n1, n2, n3, n4) for integer ni's. The general idea is to start by interpolating along each coordinate axis--that is, allow only one of the four variables to vary. This would generate a lot of interpolations--for a 4x4x4x4 lattice, it would generate 4^4 = 256 interpolations, however only a few of these are used. Combine these interpolations by successively creating new interpolations based on intermediate points in these coordinate axis interpolations. After several new interpolations, you'll have your point. I can be more specific if it would help. I would bet that the method I just described is equivalent to an existing method, but I'm not into computer graphics enough to know specifically.

Hi Jeidiah.
Thanks for the help; however, I could read what you wrote for the next 2 weeks and I promise you, I will never understand it. You introduced new terms, new concepts and I wont even bother to ask for further clarification. Thing is, how do you generalise a multivariate case from univariate cases using a tensor product construct ? That's the one and only answer I seek, anything else wont help --- but thanks for your consideration and time, I appreciate it. And it will really help if someone could show an example using polynomials (multiplication of the univariate cases to produce multivariate cases).

Re: Cubic Spline Tutorial

I was interested in the algorithm I briefly described above and wrote it up here. It extends cubic spline interpolation (as described in this thread) to higher dimensions. No tensors were used, though I have to say it's inefficient unless you keep the lattice small. It may or may not be useful to you; I wrote it up mostly for my own interest.

Re: Cubic Spline Tutorial

Hi, I realise I'm coming to the party late, but I've just found this post and it's super helpful! I just have one question though, and apologies if it's quite obvious. I'm not clear on how you jump from the 8x8 matrix to the 3x3 matrix. You say that we have the same boundary conditions y1'' = y3'' = 0, but then say that the only unknowns are the 2nd derivatives at each point. Please could you kindly explain a bit further what it is I'm missing? Thank you so much.

Re: Cubic Spline Tutorial

He noted that if you rewrite the form of the spline a bit, you can determine the other parameters easily from just the second derivatives. The relevant equations are in red (a1 through d1 equaling stuff involving only the second derivative). It's a clever way to break the original problem up into smaller pieces. Step 1 is computing the second derivatives, which turns out to be relatively straightforward--solve the tridiagonal matrix equation. Step 2 is also easy--use the results of step 1 to determine the coefficients of the various splines, which were derived explicitly in the post.

Re: Cubic Spline Tutorial

Hi VBAhack, this information is great but unfortunately I am having trouble creating a VBA spline function. I have found some examples online, some which allow me to specify end point conditions such as natural or a specific value for slope or curvature, however I want to use parabolic and cubic runouts for the endpoints but cannot find out how to code this. Do you happen to have a VBA function where you can specify end point conditions which include these? In these cases the curvature is not simply a hardcoded value but is dependent upon the curvature of the portions of the spline next to the end points, so I am unsure how to code this. Any help would be greatly appreciated!

Re: Cubic Spline Tutorial

Quote:

Originally Posted by andibuck

Hi VBAhack, this information is great but unfortunately I am having trouble creating a VBA spline function. I have found some examples online, some which allow me to specify end point conditions such as natural or a specific value for slope or curvature, however I want to use parabolic and cubic runouts for the endpoints but cannot find out how to code this. Do you happen to have a VBA function where you can specify end point conditions which include these? In these cases the curvature is not simply a hardcoded value but is dependent upon the curvature of the portions of the spline next to the end points, so I am unsure how to code this. Any help would be greatly appreciated!

Hopefully you've found an answer by now, but if not, I think the best way to answer is by a simple example. Attached are the equations for (1) natural boundary conditions per explanation in my original post, and (2) parabolic runout. Note that in the case of parabolic runout, we want y1" = y2", which is the same as y1" - y2" = 0.

https://postimg.cc/jwwJ8zcx

Re: Cubic Spline Tutorial

Quote:

Originally Posted by VBAhack

Hopefully you've found an answer by now, but if not, I think the best way to answer is by a simple example. Attached are the equations for (1) natural boundary conditions per explanation in my original post, and (2) parabolic runout. Note that in the case of parabolic runout, we want y1" = y2", which is the same as y1" - y2" = 0.

https://postimg.cc/jwwJ8zcx

hi VBAHack user, this NaturalSpline is same thing NaturalInterpolation or not
do you knows nautral(neighbour)Interpolation algorithm (spatial interpolation)

VBAhack user i am need speak with you about Numerical algorithms in VBA excel, please PM to me.
is very usefull add more algorithms to AlgoLib VBA and create another library