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$.
User defined variables
Author : Crissy Ruffo
Contact : University of Texas at San Antonio, Department of Mathematics Date : October 2009 Version : 2.0
{{{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| /// }}}