もりとにーのブログ: 01/30/16

テスト期間中だけどラグランジュ補間のプログラムを作成した．で，試しにルンゲ現象を再現してみた．

ベッ……別に単位落とした時の言い訳なんかじゃないんだからねっ(｡>﹏<｡)

1.前知識

ラグランジュ補間は関数の近似を行うアルゴリズム手法の一つ．例えば値がわかっているデータが10個あったとする(例えばこんな感じ．x_data = { x0,x1,x2,...,x9}の時のy_data={y0,y1,y2,…,y9})

この時，f(x)を近似したい式として，近似式p(x)は次のように計算される．

この時，f(x)を近似したい式として，近似式p(x)は次のように計算される．
N(i,x)はxがi番目のデータに近い数字なら1を返す．もしxがi番目のサンプルデータに近ければ，i番目の時だけN(i,x)が1になるように出来てる．

こいつの利点としては雨の雨量みたいな離散的なデータを積分したいときに，そのサンプリングした奴の間のデータを出せたりする．

2.プログラム

とりあえずプログラムはこんな感じ．例としてsin関数の近似をやってみた．

	#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

	# プロットデータを作成(0から2*piから6個)
	x_data = np.linspace(0,2*np.pi,num=6)
	y_data = np.sin(x_data)

	# 試しに0から2*piまでのサイン関数を補間させてみる
	x = np.linspace(0,2*np.pi,num=100)
	y = np.array( [ approx(xi,x_data,y_data) for xi in x ] )

	# 表示
	plt.plot(x,y,label="approximation")
	plt.plot(x,np.sin(x),label="sin(x)")
	plt.title("Lagrange interpolation of sin(x)")
	plt.legend()
	plt.show()

view raw LagrangeInterpolation.py hosted with ❤ by GitHub

で，結果はこんな感じ．

6個のデータしか与えてないのにだいぶいい感じに近似出来てる．すごい．まぁサンプルデータが偏ったりしているとうまくいかない．あとルンゲ現象なんてものもある．ちょっとやってみますか．

3.ルンゲ現象

Wikipediaに載ってる式を試しにやってみます．この式です．

プログラムはこいつ．

	#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+25x_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.0x2),label="1/(1+25x**2)")
	plt.title("Lagrange interpolation of 1/(1+25x*2)")
	plt.ylim(-0.5,1.5)
	plt.legend()
	plt.show()