ラグランジュ補間でルンゲ現象(Runge's phenomenon）の再現

RungesPhenomenon.py

#coding:utf-8

import numpy as np
import matplotlib.pyplot as plt

def approx(x,x_data,y_data):
	""" x_data,y_dataを利用してy(x)の値をラグランジュ補間で予測する """

	y = 0						# 予想されるf(x)の近似値
	n_data = len(x_data)		# サンプルデータの数

	for i in range(n_data):        
	    # indxにはi以外の0からn_data-1が代入される
	    indx = np.arange(n_data)[ np.arange(n_data) != i ]
	    
	    # N(i,x)に当たる．np.productは全配列データの積を返す
	    N = np.product(x-x_data[indx])/np.product(x_data[i]-x_data[indx])
	    
	    # f(x_data[i])*N(i,x)を足していく
	    y += y_data[i]*N
	return y

for n in [10,15,20,25]:
	# プロットデータを作成(0から2*piからn個)
	x_data = np.linspace(-1.0,1.0,num=n)
	y_data = 1/(1+25*x_data**2)

	# 補間
	x = np.linspace(-1.0,1.0,num=100)
	y = np.array( [ approx(xi,x_data,y_data) for xi in x ] )

	# プロット
	plt.plot(x,y,label="approx,n="+str(n))

plt.plot(x,1/(1+25.0*x**2),label="1/(1+25*x**2)")
plt.title("Lagrange interpolation of 1/(1+25*x**2)")
plt.ylim(-0.5,1.5)
plt.legend()
plt.show()