アークタンジェントによる多項式補間法のPython実装

参考文献

数値計算の基礎と応用[新訂版]
数値解析学への入門
杉浦　洋(南山大学教授)　著
発行日 2009/12/10

準備

オンラインコンパイラを使用します。

ソースコード

sample.py

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt  # グラフ描画のために追加

def p(x, n, xi, b):
    y = b[n]
    for l in range(n-1, -1, -1):
        y = (x - xi[l]) * y + b[l]
    return y

def f(x):
    return np.exp(x)

def NewtonCoef(xi, m, b):
    for n in range(m+1):
        b[n] = f(xi[n])
        for l in range(n):
            b[n] = (b[n] - b[l]) / (xi[n] - xi[l])
    return 0

def main():
    m = 7
    Pi = np.pi
    dt = Pi / (m + 1)
    xi = [0.5 * np.atan((i + 0.5) * dt) for i in range(m+1)]
    b = [0] * (m+1)
    NewtonCoef(xi, m, b)
    print("degree={}".format(m))
    npoints = 100  # グラフ描画のために点数を増やす
    dx = 1.0 / npoints
    x_vals = []
    y_vals = []
    y_actual = []
    for i in range(npoints+1):
        x = -0.5 + i * dx
        y = p(x, m, xi, b)
        x_vals.append(x)
        y_vals.append(y)
        y_actual.append(f(x))
        print("p({:4.1f})={:17.10e} error={:9.2e}".format(x, y, y - np.exp(x)))
    
    # グラフ描画
    plt.plot(x_vals, y_vals, label='Polynomial Approximation')
    plt.plot(x_vals, y_actual, label='Actual exp(x)', linestyle='dashed')
    plt.legend()
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Polynomial Approximation vs Actual exp(x)')
    plt.show()

main()

実行結果

console

degree=7
p(-0.5)= 6.0651199761e-01 error=-1.87e-05
p(-0.5)= 6.1260930349e-01 error=-1.71e-05
p(-0.5)= 6.1876775669e-01 error=-1.56e-05
p(-0.5)= 6.2498798021e-01 error=-1.43e-05
p(-0.5)= 6.3127060289e-01 error=-1.30e-05
p(-0.5)= 6.3761625943e-01 error=-1.19e-05
p(-0.4)= 6.4402559049e-01 error=-1.08e-05
p(-0.4)= 6.5049924280e-01 error=-9.85e-06
p(-0.4)= 6.5703786915e-01 error=-8.95e-06
p(-0.4)= 6.6364212859e-01 error=-8.12e-06
p(-0.4)= 6.7031268640e-01 error=-7.36e-06
p(-0.4)= 6.7705021423e-01 error=-6.66e-06
p(-0.4)= 6.8385539018e-01 error=-6.02e-06
p(-0.4)= 6.9072889883e-01 error=-5.43e-06
p(-0.4)= 6.9767143138e-01 error=-4.89e-06
p(-0.3)= 7.0468368570e-01 error=-4.40e-06
p(-0.3)= 7.1176636640e-01 error=-3.96e-06
p(-0.3)= 7.1892018495e-01 error=-3.55e-06
p(-0.3)= 7.2614585972e-01 error=-3.18e-06
p(-0.3)= 7.3344411608e-01 error=-2.84e-06
p(-0.3)= 7.4081568647e-01 error=-2.53e-06
p(-0.3)= 7.4826131052e-01 error=-2.26e-06
p(-0.3)= 7.5578173508e-01 error=-2.01e-06
p(-0.3)= 7.6337771434e-01 error=-1.78e-06
p(-0.3)= 7.7105000989e-01 error=-1.58e-06
p(-0.2)= 7.7879939083e-01 error=-1.39e-06
p(-0.2)= 7.8662663383e-01 error=-1.23e-06
p(-0.2)= 7.9453252322e-01 error=-1.08e-06
p(-0.2)= 8.0251785109e-01 error=-9.47e-07
p(-0.2)= 8.1058341737e-01 error=-8.29e-07
p(-0.2)= 8.1873002988e-01 error=-7.23e-07
p(-0.2)= 8.2695850449e-01 error=-6.29e-07
p(-0.2)= 8.3526966514e-01 error=-5.46e-07
p(-0.2)= 8.4366434396e-01 error=-4.73e-07
p(-0.2)= 8.5214338135e-01 error=-4.08e-07
p(-0.1)= 8.6070762607e-01 error=-3.50e-07
p(-0.1)= 8.6935793534e-01 error=-3.00e-07
p(-0.1)= 8.7809517491e-01 error=-2.56e-07
p(-0.1)= 8.8692021917e-01 error=-2.18e-07
p(-0.1)= 8.9583395122e-01 error=-1.84e-07
p(-0.1)= 9.0483726300e-01 error=-1.55e-07
p(-0.1)= 9.1393105533e-01 error=-1.30e-07
p(-0.1)= 9.2311623806e-01 error=-1.08e-07
p(-0.1)= 9.3239373011e-01 error=-8.98e-08
p(-0.1)= 9.4176445962e-01 error=-7.40e-08
p(-0.0)= 9.5122936399e-01 error=-6.05e-08
p(-0.0)= 9.6078939002e-01 error=-4.91e-08
p(-0.0)= 9.7044549399e-01 error=-3.96e-08
p(-0.0)= 9.8019864175e-01 error=-3.16e-08
p(-0.0)= 9.9004980885e-01 error=-2.49e-08
p( 0.0)= 9.9999998059e-01 error=-1.94e-08
p( 0.0)= 1.0100501522e+00 error=-1.49e-08
p( 0.0)= 1.0202013288e+00 error=-1.13e-08
p( 0.0)= 1.0304545256e+00 error=-8.32e-09
p( 0.0)= 1.0408107682e+00 error=-5.99e-09
p( 0.1)= 1.0512710922e+00 error=-4.16e-09
p( 0.1)= 1.0618365438e+00 error=-2.74e-09
p( 0.1)= 1.0725081796e+00 error=-1.67e-09
p( 0.1)= 1.0832870668e+00 error=-8.68e-10
p( 0.1)= 1.0941742834e+00 error=-2.93e-10
p( 0.1)= 1.1051709182e+00 error= 1.05e-10
p( 0.1)= 1.1162780708e+00 error= 3.66e-10
p( 0.1)= 1.1274968521e+00 error= 5.20e-10
p( 0.1)= 1.1388283839e+00 error= 5.95e-10
p( 0.1)= 1.1502737995e+00 error= 6.13e-10
p( 0.2)= 1.1618342433e+00 error= 5.90e-10
p( 0.2)= 1.1735108715e+00 error= 5.42e-10
p( 0.2)= 1.1853048518e+00 error= 4.78e-10
p( 0.2)= 1.1972173635e+00 error= 4.07e-10
p( 0.2)= 1.2092495980e+00 error= 3.35e-10
p( 0.2)= 1.2214027584e+00 error= 2.66e-10
p( 0.2)= 1.2336780602e+00 error= 2.03e-10
p( 0.2)= 1.2460767307e+00 error= 1.48e-10
p( 0.2)= 1.2586000100e+00 error= 1.01e-10
p( 0.2)= 1.2712491504e+00 error= 6.27e-11
p( 0.2)= 1.2840254167e+00 error= 3.28e-11
p( 0.3)= 1.2969300867e+00 error= 1.04e-11
p( 0.3)= 1.3099644507e+00 error=-5.32e-12
p( 0.3)= 1.3231298123e+00 error=-1.55e-11
p( 0.3)= 1.3364274880e+00 error=-2.11e-11
p( 0.3)= 1.3498588076e+00 error=-2.32e-11
p( 0.3)= 1.3634251141e+00 error=-2.27e-11
p( 0.3)= 1.3771277643e+00 error=-2.05e-11
p( 0.3)= 1.3909681284e+00 error=-1.72e-11
p( 0.3)= 1.4049475906e+00 error=-1.35e-11
p( 0.3)= 1.4190675486e+00 error=-9.83e-12
p( 0.4)= 1.4333294146e+00 error=-6.45e-12
p( 0.4)= 1.4477346147e+00 error=-3.58e-12
p( 0.4)= 1.4622845894e+00 error=-1.34e-12
p( 0.4)= 1.4769807939e+00 error= 2.56e-13
p( 0.4)= 1.4918246976e+00 error= 1.25e-12
p( 0.4)= 1.5068177851e+00 error= 1.73e-12
p( 0.4)= 1.5219615556e+00 error= 1.80e-12
p( 0.4)= 1.5372575235e+00 error= 1.58e-12
p( 0.4)= 1.5527072185e+00 error= 1.21e-12
p( 0.5)= 1.5683121855e+00 error= 7.67e-13
p( 0.5)= 1.5840739850e+00 error= 3.53e-13
p( 0.5)= 1.5999941932e+00 error= 2.42e-14
p( 0.5)= 1.6160744022e+00 error=-1.90e-13
p( 0.5)= 1.6323162200e+00 error=-2.84e-13
p( 0.5)= 1.6487212707e+00 error=-2.75e-13

[Execution complete with exit code 0]