ハイパボリックサインによる多項式補間法の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.sinh((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)= 5.9785719632e-01 error=-8.67e-03
p(-0.5)= 6.0444467816e-01 error=-8.18e-03
p(-0.5)= 6.1107043908e-01 error=-7.71e-03
p(-0.5)= 6.1773595645e-01 error=-7.27e-03
p(-0.5)= 6.2444268892e-01 error=-6.84e-03
p(-0.5)= 6.3119207693e-01 error=-6.44e-03
p(-0.4)= 6.3798554334e-01 error=-6.05e-03
p(-0.4)= 6.4482449399e-01 error=-5.68e-03
p(-0.4)= 6.5171031828e-01 error=-5.34e-03
p(-0.4)= 6.5864438969e-01 error=-5.01e-03
p(-0.4)= 6.6562806643e-01 error=-4.69e-03
p(-0.4)= 6.7266269187e-01 error=-4.39e-03
p(-0.4)= 6.7974959520e-01 error=-4.11e-03
p(-0.4)= 6.8689009188e-01 error=-3.84e-03
p(-0.4)= 6.9408548421e-01 error=-3.59e-03
p(-0.3)= 7.0133706185e-01 error=-3.35e-03
p(-0.3)= 7.0864610235e-01 error=-3.12e-03
p(-0.3)= 7.1601387163e-01 error=-2.91e-03
p(-0.3)= 7.2344162452e-01 error=-2.71e-03
p(-0.3)= 7.3093060525e-01 error=-2.52e-03
p(-0.3)= 7.3848204792e-01 error=-2.34e-03
p(-0.3)= 7.4609717705e-01 error=-2.17e-03
p(-0.3)= 7.5377720799e-01 error=-2.01e-03
p(-0.3)= 7.6152334746e-01 error=-1.86e-03
p(-0.3)= 7.6933679396e-01 error=-1.71e-03
p(-0.2)= 7.7721873831e-01 error=-1.58e-03
p(-0.2)= 7.8517036403e-01 error=-1.46e-03
p(-0.2)= 7.9319284787e-01 error=-1.34e-03
p(-0.2)= 8.0128736021e-01 error=-1.23e-03
p(-0.2)= 8.0945506552e-01 error=-1.13e-03
p(-0.2)= 8.1769712281e-01 error=-1.03e-03
p(-0.2)= 8.2601468603e-01 error=-9.44e-04
p(-0.2)= 8.3440890454e-01 error=-8.61e-04
p(-0.2)= 8.4288092352e-01 error=-7.84e-04
p(-0.2)= 8.5143188437e-01 error=-7.12e-04
p(-0.1)= 8.6006292515e-01 error=-6.45e-04
p(-0.1)= 8.6877518095e-01 error=-5.83e-04
p(-0.1)= 8.7756978436e-01 error=-5.26e-04
p(-0.1)= 8.8644786581e-01 error=-4.73e-04
p(-0.1)= 8.9541055396e-01 error=-4.24e-04
p(-0.1)= 9.0445897615e-01 error=-3.78e-04
p(-0.1)= 9.1359425872e-01 error=-3.37e-04
p(-0.1)= 9.2281752745e-01 error=-2.99e-04
p(-0.1)= 9.3212990787e-01 error=-2.64e-04
p(-0.1)= 9.4153252568e-01 error=-2.32e-04
p(-0.0)= 9.5102650710e-01 error=-2.03e-04
p(-0.0)= 9.6061297924e-01 error=-1.76e-04
p(-0.0)= 9.7029307045e-01 error=-1.52e-04
p(-0.0)= 9.8006791069e-01 error=-1.31e-04
p(-0.0)= 9.8993863184e-01 error=-1.11e-04
p( 0.0)= 9.9990636810e-01 error=-9.36e-05
p( 0.0)= 1.0099722563e+00 error=-7.79e-05
p( 0.0)= 1.0201374362e+00 error=-6.39e-05
p( 0.0)= 1.0304030510e+00 error=-5.15e-05
p( 0.0)= 1.0407702474e+00 error=-4.05e-05
p( 0.1)= 1.0512401762e+00 error=-3.09e-05
p( 0.1)= 1.0618139922e+00 error=-2.26e-05
p( 0.1)= 1.0724928553e+00 error=-1.53e-05
p( 0.1)= 1.0832779298e+00 error=-9.14e-06
p( 0.1)= 1.0941703855e+00 error=-3.90e-06
p( 0.1)= 1.1051713979e+00 error= 4.80e-07
p( 0.1)= 1.1162821480e+00 error= 4.08e-06
p( 0.1)= 1.1275038233e+00 error= 6.97e-06
p( 0.1)= 1.1388376174e+00 error= 9.23e-06
p( 0.1)= 1.1502847309e+00 error= 1.09e-05
p( 0.2)= 1.1618463714e+00 error= 1.21e-05
p( 0.2)= 1.1735237538e+00 error= 1.29e-05
p( 0.2)= 1.1853181005e+00 error= 1.32e-05
p( 0.2)= 1.1972306420e+00 error= 1.33e-05
p( 0.2)= 1.2092626170e+00 error= 1.30e-05
p( 0.2)= 1.2214152723e+00 error= 1.25e-05
p( 0.2)= 1.2336898639e+00 error= 1.18e-05
p( 0.2)= 1.2460876566e+00 error= 1.09e-05
p( 0.2)= 1.2586099242e+00 error= 9.91e-06
p( 0.2)= 1.2712579506e+00 error= 8.80e-06
p( 0.2)= 1.2840330291e+00 error= 7.61e-06
p( 0.3)= 1.2969364631e+00 error= 6.38e-06
p( 0.3)= 1.3099695665e+00 error= 5.12e-06
p( 0.3)= 1.3231336636e+00 error= 3.85e-06
p( 0.3)= 1.3364300897e+00 error= 2.60e-06
p( 0.3)= 1.3498601911e+00 error= 1.38e-06
p( 0.3)= 1.3634253254e+00 error= 2.11e-07
p( 0.3)= 1.3771268620e+00 error=-9.02e-07
p( 0.3)= 1.3909661818e+00 error=-1.95e-06
p( 0.3)= 1.4049446781e+00 error=-2.91e-06
p( 0.3)= 1.4190637563e+00 error=-3.79e-06
p( 0.4)= 1.4333248345e+00 error=-4.58e-06
p( 0.4)= 1.4477293436e+00 error=-5.27e-06
p( 0.4)= 1.4622787274e+00 error=-5.86e-06
p( 0.4)= 1.4769744431e+00 error=-6.35e-06
p( 0.4)= 1.4918179613e+00 error=-6.74e-06
p( 0.4)= 1.5068107666e+00 error=-7.02e-06
p( 0.4)= 1.5219543573e+00 error=-7.20e-06
p( 0.4)= 1.5372502460e+00 error=-7.28e-06
p( 0.4)= 1.5526999596e+00 error=-7.26e-06
p( 0.5)= 1.5683050400e+00 error=-7.15e-06
p( 0.5)= 1.5840670435e+00 error=-6.94e-06
p( 0.5)= 1.5999875420e+00 error=-6.65e-06
p( 0.5)= 1.6160681222e+00 error=-6.28e-06
p( 0.5)= 1.6323103869e+00 error=-5.83e-06
p( 0.5)= 1.6487159543e+00 error=-5.32e-06

[Execution complete with exit code 0