2340.sagews

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

print
natural_spline(X=[1,2,3,4,5],DATA=[28,210,1120,4760,17136],usefunct=false,f=0,approx=false,p=0,prec=15)
natural_spline(X=[12,24,36,48,60,72],DATA=[10400660,7124934600,5.25194*10^11,1.313244*10^13,1.725416*10^14,1.4776*10^15],usefunct=false,f=0,approx=false,p=0,prec=15)
natural_spline(X=[1,2,3,4,5,12,24,36,48,60,72],DATA=[28,210,1120,4760,17136,10400660,7124934600,5.25194*10^11,1.313244*10^13,1.725416*10^14,1.4776*10^15],usefunct=false,f=0,approx=false,p=0,prec=15)
print

S_0(x) = 156.0*(x - 1)^3 + 26.00*x + 2.000

S_1(x) = -52.00*(x - 2)^3 + 468.0*(x - 2)^2 + 494.0*x - 778.0

S_2(x) = 2054.*(x - 3)^3 + 312.0*(x - 3)^2 + 1274.*x - 2702.

S_3(x) = -2158.*(x - 4)^3 + 6474.*(x - 4)^2 + 8060.*x - 27480.

S_0(x) = -1.969e9*(x - 12)^3 + 2.841e11*x - 3.409e12

S_1(x) = 1.014e10*(x - 24)^3 - 7.088e10*(x - 24)^2 - 5.665e11*x + 1.360e13

S_2(x) = -3.190e10*(x - 36)^3 + 2.942e11*(x - 36)^2 + 2.113e12*x - 7.555e13

S_3(x) = 1.954e11*(x - 48)^3 - 8.540e11*(x - 48)^2 - 4.605e12*x + 2.342e14

S_4(x) = -1.717e11*(x - 60)^3 + 6.180e12*(x - 60)^2 + 5.931e13*x - 3.386e15

S_0(x) = 7.128e7*(x - 1)^3 - 7.128e7*x + 7.128e7

S_1(x) = -3.564e8*(x - 2)^3 + 2.138e8*(x - 2)^2 + 1.426e8*x - 2.851e8

S_2(x) = 1.354e9*(x - 3)^3 - 8.553e8*(x - 3)^2 - 4.990e8*x + 1.497e9

S_3(x) = -5.061e9*(x - 4)^3 + 3.208e9*(x - 4)^2 + 1.853e9*x - 7.413e9

S_4(x) = 1.852e9*(x - 5)^3 - 1.198e10*(x - 5)^2 - 6.915e9*x + 3.458e10

S_5(x) = -2.917e9*(x - 12)^3 + 2.691e10*(x - 12)^2 + 9.766e10*x - 1.172e12

S_6(x) = 1.040e10*(x - 24)^3 - 7.810e10*(x - 24)^2 - 5.165e11*x + 1.240e13

S_7(x) = -3.196e10*(x - 36)^3 + 2.961e11*(x - 36)^2 + 2.100e12*x - 7.506e13

S_8(x) = 1.954e11*(x - 48)^3 - 8.545e11*(x - 48)^2 - 4.602e12*x + 2.340e14

S_9(x) = -1.717e11*(x - 60)^3 + 6.180e12*(x - 60)^2 + 5.931e13*x - 3.386e15