テイラー展開 / 数値計算をする

一往計算できるようにはなったが、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の組み込み函数で計算した結果: