Cubic Spline (natural) with Sage

A spline is a interpolating function. When data has a cubic shape (i.e. if plotted resembles the graph of a cubic polynomial) one can try to model the data using a cubic spline. A natural cubic spline must have a knot at every point and if x is an endpoint, the f''(x) = 0. The following shows how to find a cubic spline in Sage with the five points given by: X = (0, 1.5, 3, 4.5, 6), Y = (1, 4, .5, 0, 3.5). I had some help from here:

1:  import numpy as np  
2:    
3:  # Array Array Array -> Array Array Array  
4:  # interp. helper function whose input is three parameters (from data_curvtr function) and output is three parameters  
5:  # to give to the first3solv function  
6:    
7:  def first3(c,d,e):  
8:    n = d.shape[0]  
9:    for i in xrange(1,n):  
10:      kappa = c[i-1]/d[i-1]  
11:      d[i] = d[i] - kappa*e[i-1]  
12:      c[i-1] = kappa  
13:    return c,d,e  
14:    
15:  # Array Array Array Array -> Array  
16:  # interp. helper function whose input is four arrays (three from the first3 function, and one from the data_curvtr  # function) whose output is a single array  
17:    
18:  def first3solv(c,d,e,b):  
19:    n = d.shape[0]  
20:    for i in xrange(1,n):  
21:      b[i] = b[i] - c[i-1]*b[i-1]  
22:    
23:    b[n-1] = b[n-1]/d[n-1]  
24:    for i in xrange(n-2,-1,-1):  
25:      b[i] = (b[i] - e[i]*b[i+1])/d[i]  
26:    
27:    return b  
28:    
29:  # Array Array -> Array  
30:  # interp. consumes a list of x-values and a list of y-values and returns the curvatures to be used in the cubic spline  
31:  # algorithm  
32:          
33:  def data_curvtr(xArray,yArray):  
34:    n = len(xArray) - 1  
35:    c = np.zeros(n)  
36:    d = np.ones(n+1)  
37:    e = np.zeros(n)  
38:    k = np.zeros(n+1)  
39:    c[0:n-1]= xArray[0:n-1] - xArray[1:n]  
40:    d[1:n] = 2.0*(xArray[0:n-1] - xArray[2:n+1])  
41:    e[1:n] = xArray[1:n] - xArray[2:n+1]  
42:    k[1:n] = 6.0*(yArray[0:n-1] - yArray[1:n]) \  
43:      /(xArray[0:n-1] - xArray[1:n]) \  
44:      -6.0*(yArray[1:n] - yArray[2:n+1]) \  
45:      /(xArray[1:n] - xArray[2:n+1])  
46:    first3(c,d,e)  
47:    print "c =",c,"d =",d,"e =",e  
48:    first3solv(c,d,e,k)  
49:    return k  
50:    
51:  # Array Array Array Number -> Number  
52:  # interp. consumes the data lists, the curvatures, and an x-value, and carries out the cubic spline algorithm   
53:  # to produce the height of the spline at the x-value  
54:    
55:    
56:  def get_spline(xArray,yArray,k,x):  
57:    def get_segment(xArray,x):  
58:      i_left = 0  
59:      i_right = len(xArray)- 1  
60:      while 1:  
61:        if (i_right-i_left)<=1: return i_left  
62:        i=(i_left+i_right)//2  
63:        if x < xArray[i]: i_right=i  
64:        else: i_left=i  
65:    i=get_segment(xArray,x)  
66:    h = xArray[i]-xArray[i+1]   
67:    y = ((x-xArray[i+1])**3/h-(x-xArray[i+1])*h)*k[i]/6.0 \  
68:      - ((x-xArray[i])**3/h-(x-xArray[i])*h)*k[i+1]/6.0 \  
69:      + (yArray[i]*(x-xArray[i+1])           \  
70:      - yArray[i+1]*(x-xArray[i]))/h  
71:    return y   
72:    
73:  xA1=np.array([0.,1.5,3.,4.5,6.])  
74:  yA1=np.array([1.,4.,.5,0.,3.5])  
75:  print "X =",xA1,"Y =",yA1   
76:    
77:  k= data_curvtr(xA1,yA1)  
78:  print "k =",k  
79:    
80:  y= get_spline(xA1,yA1,k,.5)  
81:  print "y =",y  
82:    
83:  var ( 'x t' )  
84:  n=len(xA1)  
85:  P=list_plot(zip(xA1,yA1), color="green", size=30)  
86:  C=sum(point((t,get_spline(xA1,yA1,k,t))) for t in [xA1[0]..xA1[n-1], step=0.1])  
87:  show(P+C)  

In the above code I have found a cubic spline for the following five points: X = (0, 1.5, 3, 4.5, 6), Y = (1, 4, .5, 0, 3.5) (one can see these entered in lines 73 and 74 above). When run in sage you get:

 X = [ 0.  1.5 3.  4.5 6. ] Y = [ 1.  4.  0.5 0.  3.5]  
 c = [-1.5     0.25    0.26666667 -0.    ] d = [ 1.  -6.  -5.625 -5.6  1.  ] e = [ 0. -1.5 -1.5 -1.5]  
 k = [ 0.     -5.02380952 2.76190476 1.97619048 0.    ]  
 y = 2.5582010582