Spline Interpolation using matrix math.

Started at 2 pm yesterday and finished up at at 1:50 am today.  Finished it after 3 hours today.  Total of 15 hours, with about 8 just implementing the matrix library.

https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/04_Spline.c
https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/matrix.c
https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/matrix.h

Compile with the following command:

gcc matrix.c 04_Spline.c

I did the matrix work earlier in support of this project.  That article is here: http://mystry-geek.blogspot.com/2014/02/a-small-c-matrix-library.html

The output and inline comments in the code step through the processing that it is doing.  I tend to keep the verbose output until I can get a good test harness in place so that I know I am not breaking things.

Output is:

time ./a.out 

h:     2.0000  2.0000  2.0000   1.0000   2.0000  

x: 1.0000   3.0000  5.0000  7.0000  8.0000  10.0000 
a: 6.0000   2.0000  5.0000 10.0000  4.0000   1.0000 
A = 
    1.000000000     0.000000000     0.000000000     0.000000000     0.000000000     0.000000000 
    2.000000000     8.000000000     2.000000000     0.000000000     0.000000000     0.000000000 
    0.000000000     2.000000000     8.000000000     2.000000000     0.000000000     0.000000000 
    0.000000000     0.000000000     2.000000000     6.000000000     1.000000000     0.000000000 
    0.000000000     0.000000000     0.000000000     1.000000000     6.000000000     2.000000000 
    0.000000000     0.000000000     0.000000000     0.000000000     0.000000000     1.000000000 

b_vector = 
    0.000000000 
   10.500000000 
    3.000000000 
  -25.500000000 
   13.500000000 
    0.000000000 

A^-1 = 
    1.000000000     0.000000000     0.000000000     0.000000000     0.000000000     0.000000000 
   -0.268343816     0.134171908    -0.036687631     0.012578616    -0.002096436     0.004192872 
    0.073375262    -0.036687631     0.146750524    -0.050314465     0.008385744    -0.016771488 
   -0.025157233     0.012578616    -0.050314465     0.188679245    -0.031446541     0.062893082 
    0.004192872    -0.002096436     0.008385744    -0.031446541     0.171907757    -0.343815514 
    0.000000000     0.000000000     0.000000000     0.000000000     0.000000000     1.000000000 

A^-1 * b_vector = c = 
    0.000000000 
    0.949685535 
    1.451257862 
   -5.254716981 
    3.125786164 
    0.000000000 

b = -2.6331  -0.7338   4.0681  -3.5388  -5.6677  
d =  0.1583   0.0836  -1.1177   2.7935  -0.5210  

Solution for x: -1.500000 with 6 points: nan 
First derivative at x: nan 

Solution for x: 1.500000 with 6 points: 4.703223 
First derivative at x: -2.514413 

Solution for x: 3.500000 with 6 points: 1.880994 
First derivative at x: 0.278629 

Solution for x: 4.500000 with 6 points: 3.318298 
First derivative at x: 2.679573 

Solution for x: 5.500000 with 6 points: 7.257174 
First derivative at x: 4.681145 

Solution for x: 6.500000 with 6 points: 10.595421 
First derivative at x: 0.877686 

Solution for x: 7.500000 with 6 points: 7.266116 
First derivative at x: -6.698375 

Solution for x: 8.500000 with 6 points: 1.882469 
First derivative at x: -2.932652 

Solution for x: 9.500000 with 6 points: 0.773192 
First derivative at x: 0.193134 

Solution for x: 15.000000 with 6 points: nan 
First derivative at x: nan 

real 0m0.009s
user 0m0.000s
sys 0m0.004s

As you can see from above the code does the following:

• Builds the (x, a) point list, with a= f(x).
• Builds the h list in points, this is the difference between 2 consecutive x values.
• Builds the A matrix from the h values.
• Builds the b vector from a formula using a and h
• Inverts the A matrix and multiplies it against the b vector to get the c_vector
• Copies the c_vector into points.
• Loads b into points using a formula involving h, j, and c.
• Loads d into points using a formula involving h and c.

Finished the following a few hours later, after sleeping:

• In the finished version we can release the matrices we generated and return the point value.
• We can reportedly generate interpolated values using the points table as is.
• Actually generate specific values.

To generate an f(x) for a given x, you have to:  

• figure out which interval between points you are in (j).
• if you are outside the bounds of the points, reject the request, return NaN.
• solve a third degree polynomial using the aj, bj, cj, dj, and xj and the provided x value.

Two ways to get the rate of change at any point

    Calculate a point on each side of the value you are looking at and find the slope of the line that connects the two points.

   You could easily program the computer to take the derivative of the simple formula for one of these lines, basically just pre-program a function that is the derivative.
  I went with calculating the derivative to get an absolutely exact estimated value. Yeah, exact estimate is a problem. :D 

---

We got a take home test with 21 points and the code I wrote calculated the result in .016 seconds right out of the box.  I was able to solve it before the class period was over.   It looks like a duck!  I also calculated the first derivative of the function so we can perform more analysis on the information. :D  Of course it took me about 20 hours to write the matrix code and this code to get to that point. :D