Pythonで関数をもっと精密に！Chebyshev補間：`polynomial.chebyshev.chebinterpolate()`で精度の高い近似を実現

numpy.polynomial.chebyshev.chebinterpolate() 関数は、Chebyshev 第一種点と呼ばれる特殊な点における関数値に基づいて、Chebyshev 多項式による関数補間を行います。Chebyshev 多項式は、区間 [−1,1] 内で数値計算に非常に効率的な性質を持つ多項式です。

関数詳細

numpy.polynomial.chebyshev.chebinterpolate(func, deg, args=())

引数

args: 関数 func に渡される追加引数のタプル。
deg: 補間多項式の次数。
func: 補間対象の関数。x を引数とし、float 型の値を返す関数である必要があります。

戻り値

coef: Chebyshev 多項式の係数を含む (deg + 1,) 形状の ndarray。係数は低次から高次へと並んでいます。

詳細解説

chebinterpolate() 関数は、Chebyshev 第一種点と呼ばれる等間隔で配置された点における関数値に基づいて、Chebyshev 多項式による関数補間を行います。Chebyshev 第一種点は、以下の式で計算されます。

x = cos(np.pi * np.arange(deg + 1) / (deg + 1))

これらの点における関数値 func(x) を用いて、Chebyshev 多項式の係数を最小二乗法によって求めます。

例

以下の例では、x^2 関数を Chebyshev 第一種点で補間し、その結果をプロットします。

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return x**2

deg = 10

x = np.cos(np.pi * np.arange(deg + 1) / (deg + 1))
y_interp = np.polynomial.chebyshev.chebinterpolate(f, deg)

x_plot = np.linspace(-1, 1, 100)
y_plot = y_interp(x_plot)

plt.plot(x_plot, f(x_plot), label='Original function')
plt.plot(x_plot, y_plot, label='Chebyshev interpolation')
plt.legend()
plt.show()

この例では、Chebyshev 多項式による補間が x^2 関数を十分に精度良く再現していることが確認できます。

Chebyshev 多項式の積分・微分は、numpy.polynomial.chebyshev.chebder() および numpy.polynomial.chebintegral() 関数で行えます。
Chebyshev 多項式による関数評価は、numpy.polynomial.chebyshev.chebeval() 関数で行えます。
Chebyshev 第一種点における関数値は、numpy.polynomial.chebyshev.chebpts() 関数で取得できます。

多項式近似
データフィッティング
微分方程式の解法
数値積分

例 1: x^3 関数の補間

import numpy as np
import matplotlib.pyplot as plt

def f(x):
  return x**3

deg = 15

x = np.cos(np.pi * np.arange(deg + 1) / (deg + 1))
y_interp = np.polynomial.chebyshev.chebinterpolate(f, deg)

x_plot = np.linspace(-1, 1, 100)
y_plot = y_interp(x_plot)

plt.plot(x_plot, f(x_plot), label='Original function')
plt.plot(x_plot, y_plot, label='Chebyshev interpolation')
plt.legend()
plt.show()

例 2: データフィッティング

この例では、ランダムなノイズを含むデータに対して Chebyshev 補間を行い、元のデータと比較します。

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(123)

x = np.linspace(-1, 1, 20)
y = x**3 + np.random.rand(20) * 0.1

deg = 9

y_interp = np.polynomial.chebyshev.chebinterpolate(y, deg)

x_plot = np.linspace(-1, 1, 100)
y_plot = y_interp(x_plot)

plt.plot(x, y, 'o', label='Data with noise')
plt.plot(x_plot, y_plot, label='Chebyshev interpolation')
plt.legend()
plt.show()

ラグランジュ補間

利点

コードが簡潔で分かりやすい
任意の点における関数値に基づいて補間が可能

欠点

数値誤差の影響を受けやすい
Chebyshev 補間よりも精度が劣る場合がある

例

import numpy as np
from scipy.interpolate import lagrange

def f(x):
  return x**3

x = np.array([-1, 0, 0.5, 1])
y = f(x)

interp_poly = lagrange(x, y)

x_plot = np.linspace(-1, 1, 100)
y_plot = interp_poly(x_plot)

import matplotlib.pyplot as plt
plt.plot(x_plot, f(x_plot), label='Original function')
plt.plot(x_plot, y_plot, label='Lagrange interpolation')
plt.legend()
plt.show()

スプライン補間

利点

境界条件を指定できる
Chebyshev 補間よりも滑らかな曲線を生成できる

欠点

計算コストがかかる
コードが複雑になる

例

import numpy as np
from scipy.interpolate import CubicSpline

def f(x):
  return x**3

x = np.array([-1, 0, 0.5, 1])
y = f(x)

cs = CubicSpline(x, y)

x_plot = np.linspace(-1, 1, 100)
y_plot = cs(x_plot)

import matplotlib.pyplot as plt
plt.plot(x_plot, f(x_plot), label='Original function')
plt.plot(x_plot, y_plot, label='Spline interpolation')
plt.legend()
plt.show()

ニューラルネットワーク

利点

データ量が増えるほど精度が向上する
複雑な非線形関係を捉えることができる

欠点

過学習しやすい
学習に時間がかかる

例

import tensorflow as tf
from tensorflow.keras import layers

model = tf.keras.Sequential([
  layers.Dense(128, activation='relu', input_shape=(1,)),
  layers.Dense(1)
])

model.compile(loss='mse', optimizer='adam')

x = np.linspace(-1, 1, 100)
y = x**3

model.fit(x[:, np.newaxis], y, epochs=1000, batch_size=32)

y_pred = model.predict(x[:, np.newaxis])

import matplotlib.pyplot as plt
plt.plot(x, y, label='Original function')
plt.plot(x, y_pred, label='Neural network prediction')
plt.legend()
plt.show()

ガウス過程回帰

利点

外れ値に強い
不確実性を推定できる

欠点

コードが複雑になる
計算コストがかかる

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor

def f(x):
  return x**3

x = np.linspace(-1, 1, 100)
y = f(x) + np.random.rand(100) * 0.1

gpr = GaussianProcessRegressor()
gpr.fit(x[:, np.newaxis], y)

y_pred, y_std = gpr.predict(x[:, np.newaxis])

import matplotlib.pyplot as plt
plt.plot(x, f(x), label='Original function')
plt.plot(x, y_pred, label='Gaussian process prediction')
plt.fill_between(x, y_pred - y_std, y_pred + y_std, alpha=0.2, label='Uncertainty')
plt.legend()
plt.show()