下の聯立方程式を Δc0 ～ Δc5 について sympy.solve()、numpy.linalg.solve() の両方で解いてみる。
sympy.solve() は 2 秒ほど時間がかかった。
numpy.linalg.solve() は瞬時に解が求まった。
この聯立方程式を Δc0 ～ Δc5 について解いてみる。ε の符号は ε0List[0] ～ ε0List[7] に合わせる。 x0List[0]**0*Δc0 +x0List[0]**1 *Δc1 +x0List[0]**2 *Δc2 +x0List[0]**3 *Δc3 +x0List[0]**4 *Δc4 x0List[0]**5 *Δc5 ±ε = ε0List[0], x0List[1]**0*Δc0 +x0List[1]**1 *Δc1 +x0List[1]**2 *Δc2 +x0List[1]**3 *Δc3 +x0List[1]**4 *Δc4 x0List[1]**5 *Δc5 ±ε = ε0List[1], x0List[2]**0*Δc0 +x0List[2]**1 *Δc1 +x0List[2]**2 *Δc2 +x0List[2]**3 *Δc3 +x0List[2]**4 *Δc4 x0List[2]**5 *Δc5 ±ε = ε0List[2], x0List[3]**0*Δc0 +x0List[3]**1 *Δc1 +x0List[3]**2 *Δc2 +x0List[3]**3 *Δc3 +x0List[3]**4 *Δc4 x0List[3]**5 *Δc5 ±ε = ε0List[3], x0List[4]**0*Δc0 +x0List[4]**1 *Δc1 +x0List[4]**2 *Δc2 +x0List[4]**3 *Δc3 +x0List[4]**4 *Δc4 x0List[4]**5 *Δc5 ±ε = ε0List[4], x0List[5]**0*Δc0 +x0List[5]**1 *Δc1 +x0List[5]**2 *Δc2 +x0List[5]**3 *Δc3 +x0List[5]**4 *Δc4 x0List[5]**5 *Δc5 ±ε = ε0List[5], x0List[6]**0*Δc0 +x0List[6]**1 *Δc1 +x0List[6]**2 *Δc2 +x0List[6]**3 *Δc3 +x0List[6]**4 *Δc4 x0List[6]**5 *Δc5 ±ε = ε0List[6], x0List[7]**0*Δc0 +x0List[7]**1 *Δc1 +x0List[7]**2 *Δc2 +x0List[7]**3 *Δc3 +x0List[7]**4 *Δc4 x0List[7]**5 *Δc5 ±ε = ε0List[7] ただし x0List、ε0List は下のとおりとする。 x0List = [-1., -0.8732605, -0.53591919, -0.14245605, 0.14245605, 0.53591919, 0.8732605, 1.] ε0List = [1.34240000e-04, -1.17721091e-04, 7.14099653e-05, -1.29418676e-05, 1.29418676e-05, -7.14099653e-05, 1.17721091e-04, -1.34240000e-04]
import sympy as sym import numpy as np N = 5 # 近似多項式の次数 x0List = np.array([-1. , -0.8732605 , -0.53591919 , -0.14245605 , 0.14245605 , 0.53591919 , 0.8732605 , 1. ]) ε0List = 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]) dimX = len(x0List) ######################## # sympy.solve() で解く # ######################## x, ε = sym.symbols("x, ε") Δc = np.array([sym.Symbol("Δc" + str(n)) for n in np.arange(N+1)]) SysEq = np.array([sum(np.concatenate(([Δc[n] * x0List[k]**n for n in np.arange(N+1) ], [np.sign(ε0List[k])*ε ], [-ε0List[k] ]))) for k in range(dimX)]) sol = sym.solve(SysEq) ΔcList = np.array([sol[Δc[n]] for n in np.arange(N+1)]) print(ΔcList) ############################ # np.linalg.solve() で解く # ############################ A = [np.concatenate(([x0List[k]**n for n in np.arange(N+1)], [np.sign(ε0List[k]) for i in range(1+N%2)])) # 行列のサイズを揃えるための処理 for k in np.arange(dimX)] B = ε0List sol = np.linalg.solve(A, B) ΔcList = np.round(sol[0:N+1], 15) print(ΔcList)
実行結果:
[0.0 -0.000337445128266279 0.0 0.00141211995462523 0.0 -0.00115191666560199] [ 0. -0.00033745 0. 0.00141212 -0. -0.00115192]
