スプライン補間（3次）

スプライン補間をやってみた。隣り合う2点を3次式で結ぶやつ。詳しくは次のサイトが参考になる。
・http://next1.cc.it-hiroshima.ac.jp/MULTIMEDIA/numeanal1/node16.html
・http://www.akita-nct.ac.jp/~yamamoto/lecture/2004/5E/interpolation/text/html/node3.html

def spline(x, y, a):
    l = len(x)
    matrix = []
    p = []
　　#まず最初に連立方程式を解くための行列を作る
    for i in range(l + 1):
        if i == 0:
            p.append(1.0)
        else:
            p.append(0.0)
    matrix.append(p)
    for i in range(1, l - 1):
        p = []
        for j in range(l + 1):
            p.append(0.0)
        p[i - 1] = x[i] - x[i - 1]
        p[i] = 2 * (x[i + 1] - x[i - 1])
        p[i + 1] = x[i + 1] - x[i]
        p[l] = 6*((y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]))
        matrix.append(p)
    p = []
    for i in range(l + 1):
        if i == l - 1:
            p.append(1.0)
        else:
            p.append(0.0)
    matrix.append(p)

    hakidashi(matrix, l)    #連立方程式を解く

    for i in range(l + 1):    #aがどの区間にあるかを調べる
        if a >= x[i] and a <= x[i + 1]:
            break
    h = x[i + 1] - x[i]
    s = (matrix[i+1][l]*(a-x[i])**3 - matrix[i][l]*(a-x[i+1])**3)/(6*h)+(y[i+1]/h - h*matrix[i+1][l]/6)*(a-x[i])-(y[i]/h - h*matrix[i][l]/6)*(a-x[i+1])
    return s

def hakidashi(matrix, dim):
    for i in range(dim):
        num = matrix[i][i]
        for j in range(dim + 1):
            matrix[i][j] = matrix[i][j] / num
        for j in range(dim):
            if i == j:
                pass
            else:
                a = matrix[j][i]
                for k in range(i, dim + 1):
                    matrix[j][k] = matrix[j][k]- a * matrix[i][k]

def main():
    x = [0.9,1.3,1.9,2.1,2.6,3.0,3.9,4.4,4.7,5.0,6.0,7.0,8.0,9.2]
    y = [1.3,1.5,1.85,2.1,2.6,2.7,2.4,2.15,2.05,2.1,2.25,2.3,2.25,1.95]
    f = open('data.dat', 'w')
    i = x[0]    #開始地点
    while i <= x[-1]:
        res = spline(x, y, i)    #iに対する補間値を求める
        f.write('%f %f\n' % (i, res))
        i += 0.01    #刻み幅
    f.close()

if __name__ == '__main__':
    main()

これで作られたdata.datをgnuplotで描いてみるとこんな感じになる。

赤い線がdata.datで描かれたもの。青い円形の点は元々のデータを表す。