# -*- 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