Project 1.5b_cubic_splines
system:sage


<h1 style="text-align: left;">Part II: Natural Cubic Spline</h1>
<p>&nbsp;</p>
<p>This program is used to construct the cubic spline interpolant $S$ for the function $f$, defined at the interval [x0,xn], satisfying $S"(x0)=S"(xn)$. If $f$ is unknown, the program will use a data set $f$. If a point $p$ is given, the program will approximate $f(p)$ using $S$.</p>
<p>User defined variables</p>
<ul>
<li>distinct $x-values$: the $n+1$ distinct $x$ values must be given&nbsp;</li>
<li>function&nbsp;data : you can input the values the function takes at the above $x-values$ or you can use the function to generate the data (check box)&nbsp;</li>
<li>$f$ : the function (optional)&nbsp;</li>
<li>Approximate $f(p)$ : input $x-value$ you wish use the spline to approximate the function at.&nbsp;</li>
<li>prec : the precision (in binary digits)</li>
</ul>
<p>Author : Crissy Ruffo</p>
<p>Contact : University of Texas at San Antonio, Department of Mathematics  Date : October 2009  Version : 2.0</p>

{{{id=0|
def natural_spline(X,DATA,usefunct,f,approx,p, prec):
   m=0
   N = len(X)
   a = mrange_iter([range(N)], sum)
   alpha = mrange_iter([range(N)], sum)
   b = mrange_iter([range(N)], sum)
   c = mrange_iter([range(N)], sum)
   d = mrange_iter([range(N)], sum)
   h = mrange_iter([range(N)], sum)
   l = mrange_iter([range(N)], sum)
   mu = mrange_iter([range(N)], sum)
   z = mrange_iter([range(N)], sum)
   S = mrange_iter([range(N)], sum)

   if (usefunct == true):
       for n in range(N):
           a[n] = f(x=X[n]).n(prec)
   else:
       a = DATA


   h[0] = X[1] - X[0] 
   l[0] = 1
   mu[0] = 0
   z[0] = 0

   for n in range(1, N-1):
       h[n]     = X[n+1] - X[n]
       alpha[n] = 3/h[n]*(a[n+1] - a[n]) - 3/h[n-1]*(a[n] - a[n-1])
       l[n]     = 2*(X[n+1] - X[n-1]) - h[n-1]*mu[n-1]
       mu[n]    = h[n]/l[n]
       z[n]     = (alpha[n] - h[n-1]*z[n-1])/l[n]

   n = N -1
   c[n] = 0
   l[n] = 1
   z[n] = 0
   n = n-1

   while n >= 0:
       c[n] = (z[n] - mu[n]*c[n+1]).n(prec)
       b[n] = ((a[n+1] - a[n])/h[n] - h[n]*(c[n+1] +2*c[n])/3).n(prec)
       d[n] = ((c[n+1] - c[n])/(3*h[n])).n(prec)
       S[n] = a[n] + b[n]*(x - X[n]) + c[n]*(x - X[n])^2 + d[n]*(x-X[n])^3
       n = n-1

   # Show spline function 
   for n in range(N-1):
       html("$$S_%s(x) = %s$$"%(n,S[n]))
       
   # Approximation        
   if (approx == true):
       for n in range(N-1):
           if (X[n] <= p < X[n+1]):
               m = n
               html("S(p) = %s"%S[n](x=p))
               n = N+1
            
   if (usefunct == true):
           err = abs(f(x=p) - S[m](x=p))
           html("f(p) = %s"%f(x=p))
           html("The absolute error is %s"%err)
               
   # Graph           
   if (usefunct == true):
       fig1 = plot(f, xmin = min(X), xmax = max(X), rgbcolor = (0,0,1))
       for n in range(N):
           fig1 = fig1 + point([X[n], a[n]], rgbcolor = (0,0,0), pointsize = 30)
   else:
       fig1 = point([X[0],a[0]], rgbcolor = (0,0,1), pointsize = 30)
       for n in range(1, N):
           fig1 = fig1 + point([X[n],a[n]],rgbcolor = (0,0,1), pointsize = 30)

   if(approx == true):   
       fig1 = fig1 + point([p,S[m](x=p)], rgbcolor = (1,0,0), pointsize = 30)

   for n in range(N-1):
       fig1 = fig1 +  plot(S[n], X[n], X[n+1], rgbcolor = (1,0,0), linestyle = '--', thickness = 1)

   fig1.show() 
       	
x = var('x')
@interact
def _(X = input_box(default=[-3, -2, -1, 0, 1, 2, 3], label = 'distinct x-values'),  DATA = input_box(default=[-19.95021293, -11.86466471, -5.632120558, -1.0, 2.718281828, 7.389056098, 18.08553692], label='function data'), usefunct=("Or use function to generate data",False), f = 3*x - x^2 - 2 + exp(x), approx=("Approximate f(p)", True), p = input_box(default=.5, label = 'estimate function value at p = '),prec = input_box(default=32, label = "precision")):
   natural_spline(X,DATA,usefunct,f,approx,p, prec)
///
}}}

{{{id=2|

///
}}}