自然样条计算

理论

#a,b,c,d are list of float that have same length
#b0 c0
#a1 b1 c1
#   a2 b2 c2
#      a3 b3
def Thomas(a, b, c, d):
    length = len(a)
    u = [b[0]]
    l = [0.]
    for i in range(1, length):
        l.append(a[i]/u[i-1])
        u.append(b[i]-l[i]*c[i-1])
    y = [d[0]]
    for i in range(1, length):
        y.append(d[i] - l[i]*y[i-1])
    x = [0.] * length
    x[length-1] = y[length-1] / u[length-1]
    for i in range(length-2, -1, -1):
        x[i] = (y[i] - c[i]*x[i+1])/u[i]
    return x
'''
x =

   0.65909
   1.63636
   0.90909
'''
def test_Thomas():
    a = [0., -4., -1.]
    b = [4., 4., 4.]
    c = [-1., -1., 0.]
    d = [1., 3., 2.]
    print(Thomas(a, b, c, d))

#x,y list of float, xq(x_query) float
def spline_natural(x, y, xq):
    assert(x[0] <= xq and xq <= x[-1])
    length = len(x)

    d = [0.] * length
    mu = [0.] * length
    lamda = [0.] * length
    h = [0.] * (length-1)
    f = [0.] * (length-1)
    for i in range(0, length-1):
        h[i] = x[i+1] - x[i]
    for i in range(0, length-1):
        f[i] = (y[i+1] - y[i])/h[i]
    for i in range(1, length-2):
        lamda[i] = h[i]/(h[i-1]+h[i])
    for i in range(2, length-1):
        mu[i] = h[i-1]/(h[i-1]+h[i])
    for i in range(1, length-1):
        d[i] = 6 /(h[i-1]+h[i]) * (f[i] - f[i-1])
    m = Thomas(mu[1:-1], [2]*(length-2), lamda[1:-1], d[1:-1])
    m = [0] + m + [0]
    idx = 0
    for i in range(0, length-1):
        if xq <= x[i+1]:
            idx = i
            break
    i = idx
    s = m[i] * (x[i+1] - xq)**3 /(6*h[i])
    s += m[i+1] * (xq - x[i])**3 /(6*h[i])
    s += (y[i] - m[i]*h[i]*h[i]/6) * (x[i+1] - xq)/h[i]
    s += (y[i+1] - m[i+1]*h[i]*h[i]/6) * (xq -x[i])/h[i]
    return s




def floatrange(start,stop,steps):
    ''' Computes a range of floating value.

        Input:
            start (float)  : Start value.
            end   (float)  : End value
            steps (integer): Number of values

        Output:
            A list of floats

        Example:
            >>> print floatrange(0.25, 1.3, 5)
            [0.25, 0.51249999999999996, 0.77500000000000002, 1.0375000000000001, 1.3]
    '''
    return [start+float(i)*(stop-start)/(float(steps)-1) for i in range(steps)]


def test_spline():
    x = [1, 2, 4, 5]
    y = [1, 3, 4, 2]
    x_test = floatrange(1, 5, 100)
    y_pre = []
    for x_value in x_test:
        y_pre.append(spline_natural(x,y,x_value))

    import numpy as np
    import matplotlib.pyplot as plt

    plt.figure()
    plt.plot(x,y,'r',x_test,y_pre,'b')
    plt.show()


test_spline()