参考: Interface(インターフェース) 2017年 09 月号, pp.109-111
import matplotlib.pyplot as plt import sympy as sym import numpy as np def newCoeff(c0List, x0List, ε0List): dimX = len(x0List) dimC = len(c0List) Δc = np.array([sym.Symbol("Δc" + str(n)) for n in np.arange(dimC)]) ε = sym.Symbol("ε") # 聯立方程式を解くための行列を組み立てる。 A = [np.concatenate(([x0List[k]**n for n in np.arange(dimC)], [np.sign(ε0List[k]) for i in range(1+(dimC-1)%2)])) # 正方行列にする手段として、求める必要のない未知項を 1 個または 2 個追加する。 for k in np.arange(dimX)] B = ε0List # その解を求めて、多項式の係数の補正値だけを抜き出す。 sol = np.linalg.solve(A, B) ΔcList = np.round(sol[0:dimC], 15) # その補正値で補正した新たな係数を返す。 return c0List + ΔcList ######## # test # ######## if __name__ == "__main__": # 函数 sin(pi * x / 2) を近似することにする。 x = sym.Symbol("x") f = sym.sin(sym.pi * x / 2) sinOrigin = lambda x : np.sin(np.pi * x / 2) # チェビシェフ補間で求めた近似多項式の係数 (0 次 ～ N 次) が下のとおりであったする。 c0ListOld = np.array([0.0, 1.57065736, -0.0, -0.64345777, 0.0, 0.07293465]) # その近似多項式は下のとおりである。 dim = len(c0ListOld) gOld_= sum([c0ListOld[n] * x**n for n in np.arange(dim)]) gOld = lambda x : sum([c0ListOld[n] * x**n for n in np.arange(dim)]) # 元の函数と近似式との誤差が極値を取るときの x の値が下のとおりであったとする。 x0ListOld = np.array([-1.,-0.8732605,-0.53591919,-0.14245605,0.14245605,0.53591919,0.8732605,1.]) # その x における誤差が下のとおりであったとする。 ε0ListOld = np.array([1.34240000e-04,-1.17721091e-04,7.14099653e-05,-1.29418676e-05,1.29418676e-05,-7.14099653e-05,1.17721091e-04,-1.34240000e-04]) # 補正された新たな係数 (0 次 ～ N 次) を求める。 c0ListNew = newCoeff(c0ListOld, x0ListOld, ε0ListOld) print(c0ListNew) # 新たな近似多項式は下のとおりとなる。 dim = len(c0ListNew) gNew_= sum([c0ListNew[n] * x**n for n in np.arange(dim)]) gNew = lambda x : sum([c0ListNew[n] * x**n for n in np.arange(dim)]) #""" # どれくらいミニマックスに近づいたのかを見てみる。 x_numbers = np.linspace(-1, 1, 2**8+1) sinOrigin = sinOrigin(x_numbers) # 元の sin(pi * x / 2) sinCheby = gOld(x_numbers) # チェビシェフ補間による近似多項式 sinMinMax = gNew(x_numbers) # 1 回だけ係数を補正した多項式 εCheby = sinOrigin - sinCheby # チェビシェフ近似との差 εMinMax = sinOrigin - sinMinMax # 係数を補正した多項式との差 plt.plot(x_numbers, εCheby) plt.plot(x_numbers, εMinMax) plt.legend(["εCheby", "εMinMax"]) plt.show() #""" """ # どれくらいミニマックスに近づいたのかを見てみる。 εOld = f - gOld_ # これが、チェビシェフ補間による近似式との差 εNew = f - gNew_ # これが、係数を 1 回だけ補正した近似式との差 sym.plot(εOld, εNew, (x,-1,1)) sym.plot(f, gOld, gNew, (x,-2.5,2.5)) #"""
実行結果 (0 次 ～ N 次の係数を 1 回だけ補正した係数): [ 0. 1.57031991 0. -0.64204565 0. 0.07178273]
真値と近似式との差は下のとおりである。
1 回補正しただけでほぼミニマックス (誤差の極値の大きさが全部同じである状態) になっている (オレンジ)。
最初のチェビシェフ補間だけで十分な気もする (水色)。
