[Python] 牛頓插值

插值公式為：

 

差商遞歸公式為：

# -*- coding: utf-8 -*-
#Program 0.4 Newton Interpolation

import numpy as np
import matplotlib.pyplot as plt

#遞歸求差商
def get_diff_quo(xi, fi):
	if len(xi) > 2 and len(fi) > 2:
		return (get_diff_quo(xi[:len(xi)-1], fi[:len(fi)-1]) - get_diff_quo(xi[1:len(xi)], fi[1:len(fi)])) / float(xi[0] - xi[-1])
	return (fi[0]-fi[1]) / float(xi[0]-xi[1])

#求w，使用閉包函數
def get_w(i, xi):
	def wi(x):
		result = 1.0
		for j in range(i):
			result *= (x - xi[j])
		return result
	return wi

#做插值
def get_Newton(xi, fi):
	def Newton(x):
		result = fi[0]
		for i in range(2, len(xi)):
			result += (get_diff_quo(xi[:i], fi[:i]) * get_w(i-1, xi)(x))
		return result
	return Newton

#已知結點
xn = [i for i in range(-50, 50, 10)]
fn = [i**2 for i in xn]

#插值函數
Nx = get_Newton(xn, fn)

#測試用例
tmp_x = [i for i in range(-50, 51)]
tmp_y = [Nx(i) for i in tmp_x]

#作圖
plt.plot(xn, fn, 'r*')
plt.plot(tmp_x, tmp_y, 'b-')
plt.title('Newton Interpolation')
plt.xlabel('x')
plt.ylabel('y')
plt.show()