補間法(3次スプライン関数)

# -*- coding: UTF-8 -*-
from math import *
import numpy as np

############################################
# クラスSplineの定義
#      coded by Y.Suganuma
############################################

class Spline :

	#########################
	# コンストラクタ
	#      n1 : 区間の数
	#      x1, y1 : 点の座標
	#########################

	def __init__(self, n1, x1, y1) :
			# 設定
		self.n = n1
		self.h = np.empty(self.n, np.float)
		self.b = np.empty(self.n, np.float)
		self.d = np.empty(self.n, np.float)
		self.g = np.empty(self.n, np.float)
		self.u = np.empty(self.n, np.float)
		self.q = np.empty(self.n, np.float)
		self.s = np.empty(self.n, np.float)
		self.r = np.empty(self.n+1, np.float)
		self.x = np.empty(self.n+1, np.float)
		self.y = np.empty(self.n+1, np.float)

		for i1 in range(0, self.n+1) :
			self.x[i1] = x1[i1]
			self.y[i1] = y1[i1]
			# ステップ1
		for i1 in range(0, self.n) :
			self.h[i1] = self.x[i1+1] - self.x[i1]
		for i1 in range(1, self.n) :
			self.b[i1] = 2.0 * (self.h[i1] + self.h[i1-1])
			self.d[i1] = 3.0 * ((self.y[i1+1] - self.y[i1]) / self.h[i1] - (self.y[i1] - self.y[i1-1]) / self.h[i1-1])
			# ステップ2
		self.g[1] = self.h[1] / self.b[1]
		for i1 in range(2, self.n-1) :
			self.g[i1] = self.h[i1] / (self.b[i1] - self.h[i1-1] * self.g[i1-1])
		self.u[1] = self.d[1] / self.b[1]
		for i1 in range(2, self.n) :
			self.u[i1] = (self.d[i1] - self.h[i1-1] * self.u[i1-1]) / (self.b[i1] - self.h[i1-1] * self.g[i1-1])
			# ステップ3
		self.r[0]        = 0.0
		self.r[self.n]   = 0.0
		self.r[self.n-1] = self.u[self.n-1]
		for i1 in range(self.n-2, 0, -1) :
			self.r[i1] = self.u[i1] - self.g[i1] * self.r[i1+1]
			# ステップ4
		for i1 in range(0, self.n) :
			self.q[i1] = (self.y[i1+1] - self.y[i1]) / self.h[i1] - self.h[i1] * (self.r[i1+1] + 2.0 * self.r[i1]) / 3.0
			self.s[i1] = (self.r[i1+1] - self.r[i1]) / (3.0 * self.h[i1])

	#########################
	# 補間値の計算
	#      x1 : 補間値を求める値
	#      return : 補間値
	#########################

	def value(self, x1) :
			# 区間の決定
		i = -1
		for i1 in range(1, self.n) :
			if x1 < self.x[i1] :
				i = i1 - 1
				break
		if i < 0 :
			i = self.n - 1
			# 計算
		xx = x1 - self.x[i]
		y1 = self.y[i] + xx * (self.q[i] + xx * (self.r[i]  + self.s[i] * xx))

		return y1

----------------------------------

# -*- coding: UTF-8 -*-
import numpy as np
from math import *
from function import Spline

############################################
# 補間法(3次スプライン関数)
#      coded by Y.Suganuma
############################################

			# 設定
n = 5
x = np.array([0, 1, 2, 3, 4, 5], np.float)
y = np.array([0, 1, 4, 9, 16, 25], np.float)

spline = Spline(n, x, y)
			# 実行と出力
x1 = 0.0
sp = 0.1
while x1 <= 5.01 :
	y1 = spline.value(x1)
	print(x1, y1)
	x1 += sp