ベッセル関数による多項式補間法のPython実装

参考ページ

準備

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

ソースコード

sample.py

# -*- coding: utf-8 -*-
import numpy as np
import math
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 bessel_j0(x, num_terms=50):
    """ゼロ次のベッセル関数 J0(x) を計算する"""
    x = np.asarray(x)  # 入力をNumPy配列に変換
    result = np.zeros_like(x)  # 結果を保存する配列
    for k in range(num_terms):
        term = ((-1)**k / (math.factorial(k) * math.factorial(k))) * (x/2)**(2*k)
        result += term
    return result

def main():
    m = 7
    Pi = np.pi
    dt = Pi / (m + 1)
    xi = bessel_j0(np.linspace(0, 10, m+1))  # ベッセル関数を使用
    b = [0] * (m+1)
    NewtonCoef(xi, m, b)
    print("degree={}".format(m))
    npoints = 10  # グラフ描画のために点数を増やす
    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 - f(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.0653055877e-01 error=-1.01e-07
p(-0.4)= 6.7032004068e-01 error=-5.35e-09
p(-0.3)= 7.4081822136e-01 error= 6.73e-10
p(-0.2)= 8.1873075312e-01 error= 4.34e-11
p(-0.1)= 9.0483741870e-01 error= 6.67e-10
p( 0.0)= 1.0000000001e+00 error= 1.25e-10
p( 0.1)= 1.1051709184e+00 error= 3.12e-10
p( 0.2)= 1.2214027601e+00 error= 1.96e-09
p( 0.3)= 1.3498588074e+00 error=-1.42e-10
p( 0.4)= 1.4918246872e+00 error=-1.04e-08
p( 0.5)= 1.6487212560e+00 error=-1.47e-08

[Execution complete with exit code 0]