ラグランジュ補間

ラグランジュ補間をやってみた。次のサイトが参考になる。
・ʪÍý¤Î¤«¤®¤·¤Ã¤Ý¡§·×»»ÊªÍý³Ø¡§Êä´ÖË¡¡§¥é¥°¥é¥ó¥¸¥åÊä´°
・http://next1.cc.it-hiroshima.ac.jp/MULTIMEDIA/numeanal1/node13.html

def lagrange(point, x):
    l = 0.0
    for i in range(len(point)):
        c = 1.0
        p = 1.0
        for j in range(len(point)):
            if i != j:
                c *= (x - point[j][0])    #分子
                p *= (point[i][0] - point[j][0])    #分母
        l += point[i][1] * c / p
    return l

def main():
    p = [[1.6, 2.7], [3.2, 1.5], [5.0, 2.7]]    #y=5/12*x^2-11/4*x+181/30上の3点。某教科書に例として載っていたものそのまま。
    #p = [[1,0.5],[6.2,0.02535],[-3,0.1],[-2,0.2],[9,0.01219],[0,1.0]]
    f = open('plot.dat', 'w')
    x = -50.0
    while x < 50:    #xを-50から50まで0.1刻みで変えてそのときのyを求める
        y = lagrange(p, x)
        f.write('%f %f\n' % (x, y))
        x += 0.1
    f.close()

if __name__ == '__main__':
    main()

これで作られたplot.datをgnuplotで描画すると次のグラフが得られる。

拡大したものがこちら。

うまく元の関数を描いてくれていることが分かる。

それでは、次の場合はどうなるのかを試してみた。

p = [[1,0.5],[6.2,0.02535],[-3,0.1],[-2,0.2],[9,0.01219],[0,1.0]]    #y=1/(1+x^2)上の点

結果はこちら。

拡大したものがこちら。

元の関数と全然違う関数を描いてしまっていることが分かる。

ラグランジュの補間多項式はxを±∞に飛ばすと±∞に発散していくけど、y=1/(1+x^2)はxを±∞に飛ばすと0に収束するのでこのようなずれが生じる。
たくさんの点を1度に使うと極値もたくさんできてしまうので小分けにしてラグランジュ補間を適用させるといいかもしれない。