牛顿插值法

一、简介

插值法利用函数f(x)在某区间中若干点的函数值，作出适当的特定函数，在这些点上取已知值，在区间的其他点上用这特定函数的值作为函数f(x)的近似值。

插值法具有唯一性，需要求解各阶差商。如果这特定函数是多项式，就称它为插值多项式。利用插值基函数很容易得到拉格朗日插值多项式，公式结构紧凑，在理论分析中方便，但当插值节点增减时全部插值基函数均要随之变化，整个公式也将发生变化，这在实际计算中是很不方便的，为了克服这一缺点，提出了牛顿插值法。

二、实现

# -*- coding: utf-8 -*-
"""
Created on Fri Dec 16 21:57:44 2016

    Newton插值

@author: Administrator
"""

from numpy import *
import matplotlib.pyplot as plt

# x [0,1]
def f(x):
    return sin(pi * x)

#均差
def get_order_diff_quot(xi, fi):
    if len(xi) > 2 and len(fi) > 2:
        return (get_order_diff_quot(xi[:len(xi) - 1], fi[:len(fi) - 1]) - get_order_diff_quot(xi[1:len(xi)], fi[1:len(fi)])) / float(xi[0] - xi[-1])
    return (fi[0] - fi[1]) / float(xi[0] - xi[1])

#基底
def get_wi(i,xi):
    def wi(x):
        result = 1.0
        for each in range(i):
            result *= (x - xi[each])
        return result
    return wi

def get_newton_inter(xi, fi):
    def newton_inter(x):
        result = fi[0]
        for i in range(2, len(xi)+1):
            result += (get_order_diff_quot(xi[:i], fi[:i]) * get_wi(i-1, xi)(x))
        return result
    return newton_inter

if __name__ == '__main__':
    x = linspace(0,1,10)
    fx = f(x)
    nx = get_newton_inter(x,fx)            # 获得插值函数
    draw_x = linspace(0,1,100)    
    nx_fx = nx(draw_x)
    f_fx = f(draw_x)

    fig = plt.figure(figsize=(8,4))
    ax = fig.add_subplot(111)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.scatter(x, fx, marker='o', color='r',label='interpolation point')
    ax.plot(draw_x, f_fx,color='k',linestyle=':',label='f(x)')
    ax.plot(draw_x,nx_fx,color='g',linestyle='--',label='nx(x)')
    ax.legend(loc='lower center')
    fig.show()
    fig.savefig('a.png')