gistfile1.py

#!/usr/bin/python
# -*- coding: utf-8 -*-

'''
    補間(Interpolation)
        ・ラグランジュ補間
        ・ニュートン補間
    通る点
        (0, 0.8), (1, 3.1), (3, 4.5), (6, 3.9), (7, 2.8)
'''


def newton(pt, x):
    '''
        ニュートン補間
            n 個の点を通る曲線の式を導出する
            x[0~n-1], y[0~n-1] までの差分を再帰的に求め
           各段階の1つ目の差分を次の式に用いる
            差分から得た値を用いて補間を行う

        補間式(a は差分から得た値)
            f(x) = a0 + a1(x-x0) + a2(x-x0)(x-x1) + a3(x-x0)(x-x1)(x-x2) +
                   ... + an-1(x-x0)(x-x1)...(x-x[n-2])
    '''
    # xl is xlist, yl is ylist of points
    xl, yl = zip(*pt)

    # 再帰的に差分を求める
    w = []
    w.append(yl[:])
    for i in xrange(1, len(w[0])):
        w.append([])
        for j in xrange(0, len(w[i - 1]) - 1):
            wdst, wsrc = w[i - 1][j + 1], w[i - 1][j]
            xdst, xsrc = xl[i + j], xl[j]
            w[i].append((wdst - wsrc) / (xdst - xsrc))
    a = [_[0] for _ in w]

    # オーナーの方法を用いて乗法の重複を省く
    result = a[-1]
    for i in xrange(len(xl) - 2, -1, -1):
        result = result * (x - xl[i]) + a[i]
    return result


def lagrange(pt, x):
    '''
        ラグランジュ補間
            n 個の点を通る曲線の式を導出する

        ラグランジュ補間の式:
        f(x) = (x-x[1])(x-x[2])...(x-x[n-1])
                / (x[0]-x[1])(x[0]-x[2])...(x[0]-x[n-1]) * y[0]
             + (x-x[0])(x-x[2])...(x-x[n-1])
                / (x[1]-x[0])(x[1]-x[2])...(x[0]-x[n-1]) * y[1]
             + ...
             + (x-x[0])(x-x[2])...(x-x[n-2])
                / (x[n-1]-x[0])(x[n-1]-x[1])...(x[n-1]-x[n-2]) * y[n-1]
    '''
    def product(l):
        if l:
            return l[0] * product(l[1:])
        else:
            return 1
    result = 0
    xl, yl = zip(*pt)
    for i in xrange(len(pt)):
        p = 1
        xi, yi = xl[i], yl[i]
        for j in xrange(len(pt)):
            if i == j:
                continue
            p *= (x - xl[j]) / float(xl[i] - xl[j])
        result += p * yi
    return result


def main():
    points = [(0, 0.8), (1, 3.1), (3, 4.5), (6, 3.9), (7, 2.8)]

    print 'lagrange'
    print '      x       y'
    for i in xrange(0, 15):
        print '%7.2f %7.2f' % (i * 0.5, lagrange(points, i * 0.5))

    print 'newton'
    print '      x       y'
    for i in xrange(0, 15):
        print '%7.2f %7.2f' % (i * 0.5, newton(points, i * 0.5))


if __name__ == '__main__':
    main()