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

data = [(0.0,0.0),(0.5,1.6),(1.0,2.0),(6.0,2.0),(7.0,1.5),(9.0,0.0)]
x = var('x')
R = PolynomialRing(QQ, 'x')
f = R.lagrange_polynomial(data)
print "7.5"
print
print "a."
print
print f
print
plot(f(x), (x, 0, 9)).show()
print
print "b."
print
natural_spline(X=[0.0,0.5,1.0,6.0,7.0,9.0],DATA=[0.0,1.6,2.0,2.0,1.5,0.0],usefunct=false,f=0,approx=false,p=0,prec=15)
print
print "c."
print
print "Cubic spline interpolation seems to produce more reasonable values, as compared to Lagrange interpolation, since Lagrange interpolation tends to produce polynomials that oscillate above and below the correct function, while spline (and Chebyshev) interpolation tend to minimize (if not remove) such discrepancies."
print
print "d."
print
print "Piecewise linear interpolation is not as good as cubic spline interpolation, since the latter interpolation shows a smooth curve, something that would not be as visible if piecewise linear interpolation is used."

7.5 a. 0.00568777627601157*x^5 - 0.134787108316520*x^4 + 1.12075734722794*x^3 - 3.85592316180552*x^2 + 4.86426514661809*x

b.

S_0(x) = -2.400*x^3 + 3.800*x

S_1(x) = 2.398*(x - 0.500000000000000)^3 - 3.599*(x - 0.500000000000000)^2 + 2.000*x + 0.5999

S_2(x) = -0.007541*(x - 1.00000000000000)^3 - 0.002171*(x - 1.00000000000000)^2 + 0.1994*x + 1.801

S_3(x) = 0.003165*(x - 6.00000000000000)^3 - 0.1153*(x - 6.00000000000000)^2 - 0.3879*x + 4.327

S_4(x) = 0.01763*(x - 7.00000000000000)^3 - 0.1058*(x - 7.00000000000000)^2 - 0.6089*x + 5.763

c. Cubic spline interpolation seems to produce more reasonable values, as compared to Lagrange interpolation, since Lagrange interpolation tends to produce polynomials that oscillate above and below the correct function, while spline (and Chebyshev) interpolation tend to minimize (if not remove) such discrepancies. d. Piecewise linear interpolation is not as good as cubic spline interpolation, since the latter interpolation shows a smooth curve, something that would not be as visible if piecewise linear interpolation is used.