一往計算できるようにはなったが、k次まで計算するときに毎回0次から微分し直しているので無駄が多い。
import numpy as np import matplotlib.pyplot as plt def differential(func, start, stop, dx, order=1): x_list = np.arange(start, stop + dx / 2, dx) # ここでは中心差分で微分を計算するので微分するたびにxの範囲が両端からdx/2ずつ狭くなる。 # その狭くなる分を最初に補っておく。 MARGIN = order * dx / 2 START = start - MARGIN STOP = stop + MARGIN + dx / 2 # stopパラメーターまで確実にリスト化されるよう少しだけ拡げておく。 x_list_tmp = np.arange(START, STOP, dx) y_list = func(x_list_tmp) for _ in range(order): y_list = np.diff(y_list) / dx return x_list, y_list def idx_nearest(vect, val): return np.argmin(np.abs(np.array(vect) - val)) def a_k(func, start, stop, dx, order, point=0): x_list, y_list = differential(func, start, stop, dx, order) f_k_b = y_list[idx_nearest(x_list, point)] return f_k_b / np.math.factorial(order) def taylor(func, start, stop, dx, order, point=0): x_list = np.arange(start, stop + dx / 2, dx) taylor_sum = np.zeros(int((stop - start) / dx + 1)) for k in range(order+1): _a_k = a_k(func, start, stop, dx, k, point) taylor_sum += _a_k * ((x_list - point)**k) return x_list, taylor_sum ########################################################## func = np.sin start = -2 * np.pi stop = 2 * np.pi dx = (stop - start) / 360 order = 5 point = -np.pi/2 x_list, y_list = taylor(func, start, stop, dx, order, point) plt.plot(x_list, func(x_list)) plt.plot(x_list, y_list) plt.xlim(start, stop) plt.ylim(-1.1, 1.1) plt.grid() plt.show()
実行結果:
↓ sin(x)
を、-pi/2
を中心に5次まで計算した。
↓ 同じ条件でNspireの組み込み函数で計算した結果: