Cubic Spline Tutorial-VBForums
Results 1 to 28 of 28

Thread: Cubic Spline Tutorial

  1. #1

    Thread Starter
    Fanatic Member VBAhack's Avatar
    Join Date
    Dec 2004
    Location
    Sector 000
    Posts
    585

    Cubic Spline Tutorial

    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):

    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 = 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:
    ┌							      ┐	┌  ┐   ┌     ┐
    │    .	     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 = 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.

    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....
    .
    Attached Images Attached Images  
    Last edited by VBAhack; Jul 30th, 2007 at 05:39 PM.

  2. #2

    Thread Starter
    Fanatic Member VBAhack's Avatar
    Join Date
    Dec 2004
    Location
    Sector 000
    Posts
    585

    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 = a1(x-x1)3 + b1(x-x1)2 + c1(x-x1) + d1
    y' = 3a1(x-x1)2 + 2b1(x-x1) + c1
    y" = 6a1(x-x1) + 2b1

    spline #2
    y = a2(x-x2)3 + b2(x-x2)2 + c2(x-x2) + d2
    y' = 3a2(x-x2)2 + 2b2(x-x2) + c2
    y" = 6a2(x-x2) + 2b2

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

    y1" = 6a1(x1-x1) + 2b1 = 0 + 2b1 = 2b1 ==> b1 = y1"/2
    y2" = 6a1(x2-x1) + 2b1 = 6a1h1 + y1" ==> a1 = (y2"-y1")/6h1

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

    y" = (x-x1)(y2"-y1")/(x2-x1) + y1"

    which can be verified to be correct (i.e. y" = y1" at x = x1 and y" = y2" at x = x2)

    Next, we apply the conditions that the spline pass though the points, in other words y1 = f(x1) and y2 = f(x2):

    y1 = 0 + 0 + 0 + d1
    d1 = y1

    y2 = (x2-x1)3(y2"-y1")/6h1 + y1"(x2-x1)2/2 + c1(x2-x1) + y1
    y2 = h13(y2"-y1")/6h1 + y1"h12/2 + c1h1 + y1
    y2 = h12(y2"-y1")/6 + y1"h12/2 + c1h1 + y1
    y2-y1 = y2"h12/6 - y1"h12/6 + y1"h12/2 + c1h1
    (y2-y1)/h1 = y2"h1/6 - y1"h1/6 + y1"h1/2 + c1
    (y2-y1)/h1 = y2"h1/6 - y1"h1/6 + 3y1"h1/6 + c1
    (y2-y1)/h1 = y2"h1/6 + y1"h1/3 + c1
    c1 = (y2-y1)/h1 – y2"h1/6 - y1"h1/3

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

    3a1(x2-x1)2 + 2b1(x2-x1) + c1 = 3a2(x2-x2)2 + 2b2(x2-x2) + c2
    3a1h12 + 2b1h1 + c1 = c2
    h1(y2"-y1")/2 + y1"h1 + (y2-y1)/h1 - y2"h1/6 - y1"h1/3 = (y3-y2)/h2 – y3"h2/6 – y2"h2/3
    h1(y2"-y1")/2 + y1"h1 - y2"h1/6 - y1"h1/3 + y3"h2/6 + y2"h2/3 = (y3-y2)/h2 - (y2-y1)/h1
    3h1(y2"-y1") + 6y1"h1 - y2"h1 - 2y1"h1 + y3"h2 + 2y2"h2 = 6(y3-y2)/h2 - 6(y2-y1)/h1
    3h1y2" – 3h1y1" + 6y1"h1 - y2"h1 - 2y1"h1 + y3"h2 + 2y2"h2 = 6(y3-y2)/h2 - 6(y2-y1)/h1
    y1"(6h1 – 3h1 – 2h1) + y2"(2h1 + 2h2) + y3"h2 = 6(y3-y2)/h2 - 6(y2-y1)/h1

    y1"h1 + y2"(2h1 + 2h2) + y3"h2 = 6(y3-y2)/h2 - 6(y2-y1)/h1

    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. h1 = h2 = ... = hn-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"  │           │        ?        │
    └               		┘ └     ┘           └                 ┘
    Last edited by VBAhack; Jul 30th, 2007 at 04:28 PM.

  3. #3

    Thread Starter
    Fanatic Member VBAhack's Avatar
    Join Date
    Dec 2004
    Location
    Sector 000
    Posts
    585

    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 (y1" = y3" = 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 h1 = h2 = 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 [y1", y2", y3"]T = [0, 4.5, 0]T. From this we can calculate the coefficients of the cubic spline segments:

    a1 = (y2"-y1")/6h1 = 1.5
    b1 = y1"/2 = 0
    c1 = (y2-y1)/h1 – y2"h1/6 – y1"h1/3 = -0.125
    d1 = y1 = 0

    a2 = (y3"-y2")/6h2 = -1.5
    b2 = y2"/2 = 2.25
    c2 = (y3-y2)/h2 – y3"h2/6 – y2"h2/3 = 0
    d2 = y2 = 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 = x1:

    y1' = 3a1(x1-x1)2 + 2b1(x1-x1) + c1 = 0 + 0 + c1 = c1 = (y2-y1)/h1 – y2"h1/6 – y1"h1/3
    y1"(2h1) + y2"(h1) = 6[(y2-y1)/h1 - y1']

    At x = x3:

    y3' = 3a2(x3-x2)2 + 2b2(x3-x2) + c2 = 3a2h22 + 2b2h2 + c2
    y3' = 3h22(y3"-y2")/6h2 + 2h2y2"/2 + (y3-y2)/h2 – y3"h2/6 – y2"h2/3
    y3' = 3h2(y3"-y2")/6 + h2y2" + (y3-y2)/h2 – y3"h2/6 – y2"h2/3
    h2y3"/2 - h2y2"/2 + h2y2" – y3"h2/6 – y2"h2/3 = y3' - (y3-y2)/h2
    3h2y3" - 3h2y2" + 6h2y2" – y3"h2 – 2y2"h2 = 6(y3' - (y3-y2)/h2)
    y3"(3h2 – h2) - y2"(6h2 – 3h2 – 2h2) = 6(y3' - (y3-y2)/h2)
    y2"(h2) + y3"(2h2) = 6[y3' - (y3-y2)/h2]

    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
    .

    .
    Attached Images Attached Images  
    Last edited by VBAhack; Jul 30th, 2007 at 04:29 PM.

  4. #4
    Addicted Member Rassis's Avatar
    Join Date
    Jun 2004
    Location
    Lisbon
    Posts
    248

    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.
    ...este projecto dos Deuses que os homens teimam em arruinar...

  5. #5
    New Member Math_Girl's Avatar
    Join Date
    Jul 2007
    Location
    A place called home
    Posts
    1

    Re: Cubic Spline Tutorial

    Thanks for this
    Last edited by Math_Girl; Jul 29th, 2007 at 04:04 AM.

  6. #6
    Only Slightly Obsessive jemidiah's Avatar
    Join Date
    Apr 2002
    Posts
    2,365

    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}
    The time you enjoy wasting is not wasted time.
    Bertrand Russell

    <- Remember to rate posts you find helpful.

  7. #7
    vbuggy krtxmrtz's Avatar
    Join Date
    May 2002
    Location
    [41.5622N-2.0045E]
    Posts
    5,496

    Re: Cubic Spline Tutorial

    Rassis and jemidiah, don't worry, I rated him up for you already last week
    Lottery is a tax on people who are bad at maths
    If only mosquitoes sucked fat instead of blood...
    To do is to be (Descartes). To be is to do (Sartre). To be do be do (Sinatra)

  8. #8

    Thread Starter
    Fanatic Member VBAhack's Avatar
    Join Date
    Dec 2004
    Location
    Sector 000
    Posts
    585

    Re: Cubic Spline Tutorial

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

  9. #9
    New Member
    Join Date
    Oct 2009
    Posts
    3

    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?

  10. #10

    Thread Starter
    Fanatic Member VBAhack's Avatar
    Join Date
    Dec 2004
    Location
    Sector 000
    Posts
    585

    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.
    Last edited by VBAhack; Oct 26th, 2009 at 03:26 PM.

  11. #11
    New Member
    Join Date
    Oct 2009
    Posts
    3

    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"?

  12. #12
    Fanatic Member namrekka's Avatar
    Join Date
    Feb 2005
    Location
    Netherlands
    Posts
    548

    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?

  13. #13

    Thread Starter
    Fanatic Member VBAhack's Avatar
    Join Date
    Dec 2004
    Location
    Sector 000
    Posts
    585

    Re: Cubic Spline Tutorial

    Quote Originally Posted by thegaunlet View Post
    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.

  14. #14
    New Member
    Join Date
    Oct 2009
    Posts
    3

    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 (=

  15. #15
    Only Slightly Obsessive jemidiah's Avatar
    Join Date
    Apr 2002
    Posts
    2,365

    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.
    The time you enjoy wasting is not wasted time.
    Bertrand Russell

    <- Remember to rate posts you find helpful.

  16. #16
    Junior Member
    Join Date
    Aug 2008
    Posts
    16

    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.

  17. #17
    Only Slightly Obsessive jemidiah's Avatar
    Join Date
    Apr 2002
    Posts
    2,365

    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 = f1(s) = ax(s-s0)3 + bx(s-s0)2 + cx(s-s0) + dx
    y = f2(s) = ay(s-s0)3 + by(s-s0)2 + cy(s-s0) + dy

    you would have

    w = f1(s) = aw(s-s0)3 + bw(s-s0)2 + cw(s-s0) + dw
    x = f2(s) = ax(s-s0)3 + bx(s-s0)2 + cx(s-s0) + dx
    y = f3(s) = ay(s-s0)3 + by(s-s0)2 + cy(s-s0) + dy
    z = f4(s) = az(s-s0)3 + bz(s-s0)2 + cz(s-s0) + dz

    The rest of the algebra is similar, but with more columns in your matrices and more boundary conditions.
    The time you enjoy wasting is not wasted time.
    Bertrand Russell

    <- Remember to rate posts you find helpful.

  18. #18
    Junior Member
    Join Date
    Aug 2008
    Posts
    16

    Re: Cubic Spline Tutorial

    Quote Originally Posted by jemidiah View Post
    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 = f1(s) = ax(s-s0)3 + bx(s-s0)2 + cx(s-s0) + dx
    y = f2(s) = ay(s-s0)3 + by(s-s0)2 + cy(s-s0) + dy

    you would have

    w = f1(s) = aw(s-s0)3 + bw(s-s0)2 + cw(s-s0) + dw
    x = f2(s) = ax(s-s0)3 + bx(s-s0)2 + cx(s-s0) + dx
    y = f3(s) = ay(s-s0)3 + by(s-s0)2 + cy(s-s0) + dy
    z = f4(s) = az(s-s0)3 + bz(s-s0)2 + cz(s-s0) + dz

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

  19. #19
    Only Slightly Obsessive jemidiah's Avatar
    Join Date
    Apr 2002
    Posts
    2,365

    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.
    The time you enjoy wasting is not wasted time.
    Bertrand Russell

    <- Remember to rate posts you find helpful.

  20. #20
    Junior Member
    Join Date
    Aug 2008
    Posts
    16

    Re: Cubic Spline Tutorial

    Quote Originally Posted by jemidiah View Post
    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.

  21. #21
    Only Slightly Obsessive jemidiah's Avatar
    Join Date
    Apr 2002
    Posts
    2,365

    Re: Cubic Spline Tutorial

    Quote Originally Posted by frankly747 View Post
    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.
    The time you enjoy wasting is not wasted time.
    Bertrand Russell

    <- Remember to rate posts you find helpful.

  22. #22
    Only Slightly Obsessive jemidiah's Avatar
    Join Date
    Apr 2002
    Posts
    2,365

    Re: Cubic Spline Tutorial

    Quote Originally Posted by frankly747 View Post
    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.
    The time you enjoy wasting is not wasted time.
    Bertrand Russell

    <- Remember to rate posts you find helpful.

  23. #23
    Junior Member
    Join Date
    Aug 2008
    Posts
    16

    Re: Cubic Spline Tutorial

    Quote Originally Posted by jemidiah View Post
    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.

  24. #24
    Only Slightly Obsessive jemidiah's Avatar
    Join Date
    Apr 2002
    Posts
    2,365

    Re: Cubic Spline Tutorial

    Quote Originally Posted by frankly747 View Post
    ...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.
    The time you enjoy wasting is not wasted time.
    Bertrand Russell

    <- Remember to rate posts you find helpful.

  25. #25
    Junior Member
    Join Date
    Aug 2008
    Posts
    16

    Re: Cubic Spline Tutorial

    Quote Originally Posted by jemidiah View Post
    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).

  26. #26
    Only Slightly Obsessive jemidiah's Avatar
    Join Date
    Apr 2002
    Posts
    2,365

    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.
    The time you enjoy wasting is not wasted time.
    Bertrand Russell

    <- Remember to rate posts you find helpful.

  27. #27
    New Member
    Join Date
    Sep 2013
    Location
    London
    Posts
    5

    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.

  28. #28
    Only Slightly Obsessive jemidiah's Avatar
    Join Date
    Apr 2002
    Posts
    2,365

    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.
    The time you enjoy wasting is not wasted time.
    Bertrand Russell

    <- Remember to rate posts you find helpful.

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •  



Featured


Click Here to Expand Forum to Full Width

Survey posted by VBForums.