「線形代数の半歩先」をPythonで写経 ～ １章 エルミート多項式、ルジャンドル多項式、フーリエ級数展開

第２部「ならべた数に解釈を」

書籍の著者　大久保 潤 先生

この記事は、書籍「線形代数の半歩先」の 第２部「ならべた数に解釈を」に掲載の「関数の基底展開」に関する Python写経活動 のドキュメンタリーです。

第８話、第12話の エルミート多項式、ルジャンドル多項式、フーリエ級数展開 の可視化に取り組みます！
きれいで不思議な曲線をお楽しみください！

この３つの数式は私にとって初物です。
数学素人なので、どうぞお手柔らかにお願いいたします。

では書籍を開いて線形代数の旅に出発です🚀

はじめに

このブログシリーズは、書籍「線形代数の半歩先 データサイエンス・機械学習に挑む前の30話」（講談社サイエンティフィク、「テキスト」と呼びます）の Python 写経の実践を通じて得た個人的な知見を書きます。

書籍の紹介と引用表記はリンク先の記事に掲載しています。

第２部 ならべた数に解釈を

テキストの第８話「やっぱり基底は一つではない」では、$${N}$$次の多項式で作られる空間 $${\R[x]_N}$$を検討します。
この空間は定数項を加えて $${N+1}$$ 次元なので、生成元の個数が $${N+1}$$ であることと、一次独立性であることを確認できれば、多項式が基底であるとわかる、とのことです。
そして、便利で性質がよくわかっている基底の例として、エルミート多項式とルジャンドル多項式の紹介へと進みます。
また第10話に、エルミート多項式とルジャンドル多項式の直交性に関する説明があります。

では、この記事のコードで利用するライブラリのインポートからスタートしましょう。

### インポート

# 数値計算
import numpy as np

# 数学演算用の個別関数
from scipy.special import hermite   # エルミート多項式
from scipy.special import legendre  # ルジャンドル多項式
from scipy.integrate import quad    # 定積分（数値積分）

# sympy
import sympy as sp
import spb                          # sympyの可視化ライブラリ

# 可視化
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

１．エルミート多項式

① エルミート多項式の概要
テキストは「物理学者」のエルミート多項式を紹介しています。
テキストより、次数 $${n}$$ のエルミート多項式 $${H_n(x)}$$ を生み出す式を引用いたします。

次数 $${n=0,1,2,3,4}$$ の場合のエルミート多項式がテキスト 68ページで紹介されています。

python 実装がしやすかった 三項間漸化式 も引用いたします。

② エルミート多項式の実装
テキスト 69 ページの図 8.2 「エルミート多項式のいくつか」を３つの方法で実装します。

■ エルミート多項式の三項間漸化式を用いる方法
テキストの数式（8. 8）を手がかりにして、次のWebサイトでエルミート多項式の実装方法を学びました。
ありがとうございます！

コードの「エルミート多項式関数の定義」の箇所で三項漸化式を実装しています。

### p.69 図8.2 エルミート多項式のいくつか 関数作成

## 設定と準備

# エルミート多項式関数の定義　引数 n:次数、x:変数
# 参考サイト：https://qiita.com/hanaata/items/52a1a813385380124131
def H(n, x):
    match n:
        case 0:  # H_0(x) = 1
            return np.ones(len(x)) 
        case 1:  # H_1(x) = 2x
            return 2*x
        case _:  # H_2(x)以降: 式（8.8）三項間漸化式 H_n+1(x)=2xHn(x)-2nH_n-1(x)
            return 2*x*H(n-1, x) - 2*(n-1)* H(n-2, x)

# 変数の定義
xs = np.linspace(-3, 3, 1001)  # x軸の値
ns = list(range(5))            # 次数のリスト

## 描画処理
# 描画領域の設定
fig, ax = plt.subplots()
# 次数ごとにエルミート多項式の算出とチャート描画を繰り返し処理
for n in ns:
    # エルミート多項式の算出と曲線の描画
    ax.plot(xs, H(n, xs), label=f'$H_{n}(x)$')
# 修飾
ax.set(xlim=(-3, 3), ylim=(-30, 20), title='図 8.2 エルミート多項式のいくつか')
ax.set_xlabel('$x$', fontsize=12)
ax.legend();

【実行結果】

■ numpy を利用する方法
続いて numpy の多項式を扱う numpy.polynomial の Hermite() を使って実装します。
コードの「エルミート多項式インスタンスの生成」の箇所で、エルミート多項式を生み出す Hn を設定しています。
変数 xs , ns は先のコードで作成したものをそのまま使います。

### p.69 図8.2 エルミート多項式のいくつか numpy利用

## 描画処理
# 描画領域の設定
fig, ax = plt.subplots()
# 次数ごとにエルミート多項式の算出とチャート描画を繰り返し処理
for n in ns:
    # エルミート多項式インスタンスの生成 ※numpy.polynomial.Hermite
    Hn = np.polynomial.Hermite([0] * n + [1])
    # エルミート多項式の算出と曲線の描画
    ax.plot(xs, Hn(xs), label=f'$H_{n}(x)$')
# 修飾
ax.set(xlim=(-3, 3), ylim=(-30, 20), title='図 8.2 エルミート多項式のいくつか')
ax.set_xlabel('$x$', fontsize=12)
ax.legend();

【実行結果】

Hermite() の引数には次数の係数リストを与えています。
例えば次数が 2 の場合は、[0, 0, 1] を渡しています。
0次、1次、2次の係数を順に設定しています。

■ scipy を利用する方法
最後に scipy の特殊関数を扱う scipy.special の hermite() を使って実装します。
コードの「エルミート多項式インスタンスの生成」の箇所で、エルミート多項式を生み出す Hn を設定しています。
変数 xs , ns は先のコードで作成したものをそのまま使います。

実装にあたっては次のWebサイトを参考にしました。
ありがとうございます！

### p.69 図8.2 エルミート多項式のいくつか scipy利用
# 参考サイト: https://python.atelierkobato.com/hermite/

## 描画処理
# 描画領域の設定
fig, ax = plt.subplots()
# 次数ごとにエルミート多項式の算出とチャート描画を繰り返し処理
for n in ns:
    # エルミート多項式インスタンスの生成 ※scipy.special.hermite
    Hn = hermite(n=n)
    # エルミート多項式の算出と曲線の描画
    ax.plot(xs, Hn(xs), label=f'$H_{n}(x)$')
# 修飾
ax.set(xlim=(-3, 3), ylim=(-30, 20), title='図 8.2 エルミート多項式のいくつか')
ax.set_xlabel('$x$', fontsize=12)
ax.legend();

【実行結果】

Hermite() の引数には次数 n を与えています。
numpy よりも直感的ですね！

次数のリスト ns の内容を変更して、さまざまな次数の曲線をお楽しみくださいませ。

エルミート多項式は 第11話、第23話あたりで活用されます！

２．ルジャンドル多項式

① ルジャンドル多項式の概要
テキストは次数 $${n=0, 1, 2, 3, 4}$$ のルジャンドル多項式 $${P_n(x)}$$を紹介しています。

こちらの Web サイトにて、ロドリゲスの公式と呼ばれるルジャンドル多項式の公式を教えていただきました。
ありがとうございます！

【ロドリゲスの公式】

② ルジャンドル多項式の実装
テキスト 71 ページの図 8.3 「ルジャンドル多項式のいくつか」を３つの方法で実装します。

■ ロドリゲスの公式を用いる方法
ロドリゲスの公式には微分が含まれています。
微分の得意な sympy ライブラリで実装します。
コードの「ルジャンドル多項式の算出」の箇所にロドリゲスの公式を書きました。
チャート描画は sympy と相性の良い spb (Sympy Plotting Backends) を利用します。

### p.71 図8.3 ルジャンドル多項式のいくつか sympy利用
# 参考サイト 高校数学の美しい物語： https://manabitimes.jp/math/1368

## 設定と準備
# sympy変数xの定義
x = sp.symbols('x')
# 次数のリスト
ns = list(range(5))

## ルジャンドル多項式の算出
# ルジャンドル多項式の結果を格納するリスト
Pn = []
# 次数ごとにルジャンドル多項式の算出を繰り返し処理
for n in ns:
    # ルジャンドル多項式の算出 ※参考サイトの公式参照（ロドリゲスの公式）
    leg_poly = 1/(2**n * sp.factorial(n)) * sp.diff((x**2 - 1)**n, x, n) 
    # ルジャンドル多項式の結果と凡例ラベルをリストに格納
    Pn.append((leg_poly, f'$P_{n}(x)$'))

## 描画処理
spb.plot(*Pn, (x, -1, 1), title='図 8.3 ルジャンドル多項式のいくつか',
         xlim=(-1, 1), ylim=(-1.5, 1.5), ylabel='$P_n(x)$');

【実行結果】

■ numpy を利用する方法
続いて numpy の多項式を扱う numpy.polynomial の Legendre() を使って実装します。
コードの「ルジャンドル多項式インスタンスの生成」の箇所で、ルジャンドル多項式を生み出す Pn を設定しています。
変数 ns は先のコードで作成したものをそのまま使います。

### p.71 図8.3 ルジャンドル多項式のいくつか numpy利用

## 設定と準備
# 変数の定義
xs = np.linspace(-1, 1, 1001)  # x軸の値

## 描画処理
# 描画領域の設定
fig, ax = plt.subplots()
# 次数ごとにルジャンドル多項式の算出とチャート描画を繰り返し処理
for n in ns:
    # ルジャンドル多項式インスタンスの生成 ※numpy.polynomial.Legendre
    Pn = np.polynomial.Legendre([0] * n + [1])
    # ルジャンドル多項式の算出と曲線の描画
    ax.plot(xs, Pn(xs), label=f'$P_{n}(x)$')
# 修飾
ax.set(xlim=(-1, 1), ylim=(-1.5, 1.5), xticks=np.arange(-1, 1.1, 0.5),
       title='図 8.3 ルジャンドル多項式のいくつか')
ax.set_xlabel('$x$', fontsize=12)
ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left');

【実行結果】

Legendre() の引数には、Hermite() と同様に、次数の係数リストを与えています。

■ scipy を利用する方法
最後に scipy の特殊関数を扱う scipy.special の legendre() を使って実装します。
コードの「ルジャンドル多項式インスタンスの生成」の箇所で、ルジャンドル多項式を生み出す Pn を設定しています。
変数 xs , ns は先のコードで作成したものをそのまま使います。

### p.71 図8.3 ルジャンドル多項式のいくつか scipy利用

## 描画処理
# 描画領域の設定
fig, ax = plt.subplots()
# 次数ごとにルジャンドル多項式の算出とチャート描画を繰り返し処理
for n in ns:
    # ルジャンドル多項式インスタンスの生成 ※scipy.special.legendre
    Pn = legendre(n=n)
    # ルジャンドル多項式の算出と曲線の描画
    ax.plot(xs, Pn(xs), label=f'$P_{n}(x)$')
# 修飾
ax.set(xlim=(-1, 1), ylim=(-1.5, 1.5), xticks=np.arange(-1, 1.1, 0.5),
       title='図 8.3 ルジャンドル多項式のいくつか')
ax.set_xlabel('$x$', fontsize=12)
ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left');

【実行結果】

次数のリスト ns の内容を変更して、さまざまな次数の曲線をお楽しみくださいませ。

３．フーリエ級数展開

① フーリエ級数展開の概要
第12話「波の分解と再構築」では、第２部のまとめに加えて「波の関数」が紹介されています。
特に、周期関数を $${\cos}$$ 関数や $${\sin}$$ 関数を基底として表現する フーリエ級数展開 に焦点を当てています。

フーリエ級数展開の数式をテキストよりお借りします。

係数 $${c_0, c_n, d_n}$$ の算出に関しては、次のWebサイトにお世話になりました。
ありがとうございます！

係数算出の計算式をお借りします。

フーリエ級数展開はさまざまな波を重ねることで、任意の関数を表現できるのです。
ということで、任意の関数をフーリエ級数展開で近似させるコードを作っていきましょう！

（参考情報）
次の Web サイトでは「どんな形でも三角関数のたし合わせで表すことができる」ことを丁寧に説明されています。
合わせてお読みください！

② フーリエ級数展開の実装
この記事では、公式を利用して、フーリエ級数展開でいくつかの関数 $${f(t)}$$ を近似してみようと思います。

■ フーリエ級数展開で関数 $${f(t)}$$ の近似値を求める関数の定義
テキストの式（12. 2）と参考Webサイトの公式を利用して、関数化しました。
積分計算には scipy.integrate の quad （数値積分）を利用しています。

## フーリエ級数展開で関数fの近似値を求める関数の定義
# 参考サイト：高校数学の美しい物語 https://manabitimes.jp/math/1156

# 引数 f:関数、t:時間、T:周期、N:次数
def fourier_series_expansion(f, t, T, N):

    ## cn, dnの数式内で積分する関数の定義 ※参考サイトの係数の導出
    c_func = lambda t, n, T: f(t) * np.cos(2*n*np.pi*t / T)
    d_func = lambda t, n, T: f(t) * np.sin(2*n*np.pi*t / T)

    ## c0, cn, dnの算出 ※参考サイトの公式an, bn。scipy.integrate.quadで数値積分
    c0 = 2/T * quad(func=f, a=0, b=T)[0]
    cn = [2/T * quad(c_func, a=0, b=T, args=(n, T))[0] for n in range(1, N+1)]
    dn = [2/T * quad(d_func, a=0, b=T, args=(n, T))[0] for n in range(1, N+1)]
    
    # フーリエ級数展開の計算 式（12.2）
    f_t = 1/2 * c0 + sum([
        cn[i] * np.cos(2*n*np.pi*t / T) + dn[i] * np.sin(2*n*np.pi*t / T)
        for i, n in enumerate(range(1, N+1))])
    
    # 戻り値：フーリエ級数展開によるf(t)の近似値
    return f_t

あわせて、チャート描画関数を作っておきます。

## 描画関数の定義
# 引数 f:関数、cycle:周期の数、N:次数のリスト
def plot_fourier(f, cycle, N_list):
    
    ## 定義域の設定
    # 周期 T の設定 ※周期は2πで固定
    T = np.pi*2
    # 時間変数 t の設定
    t = np.arange(0, T*cycle, 0.01)
    # 周期Tの時間を繰り返す時間変数 tcycle の設定（関数fの描画用に利用）
    tcycle = np.tile(np.arange(0, T, 0.01), reps=cycle)[:len(t)]

    ## 描画処理
    # 関数 f の描画
    plt.plot(t, f(tcycle), color='black', ls='--', label='元の関数 $f\ (t)$')
    # リストの次数 N の値ごとにフーリエ級数展開の算出と描画を繰り返し処理
    for N in N_list:
        # フーリエ級数展開の計算
        fourier = fourier_series_expansion(f, t, T, N)
        # フーリエ級数展開の描画
        plt.plot(t, fourier, lw=1, label=f'フーリエ級数展開 $N={N}$')
    # 修飾
    plt.xlabel('時間 $t$', fontsize=12)
    plt.ylabel('$f\ (t)$', fontsize=12)
    plt.title('フーリエ級数展開による関数 $f\ (t)$ の近似')
    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
    plt.grid(lw=0.5)
    plt.show()

では、３つの関数を試してみましょう。
次数 $${N}$$ の値が大きくなるにつれて、フーリエ級数展開の結果は元の関数に近づいていきます。

■ $${\cos}$$ 関数の足し合わせ関数
実はこの関数はフーリエ級数展開と同じ計算をしているので、近似というよりも一致する結果を得られます。

N_list で描画する次数を、cycle で周期の数を設定します。
周期は $${2\pi}$$ 固定としています。

## フーリエ級数展開の設定と準備
# 対象関数 f の設定 （cosの和。フーリエ級数展開と同じ計算のはず）
f = lambda t: np.cos(t) + np.cos(3*t) + np.cos(5*t)
# 次数のリストの設定
N_list = [1, 3, 5]  # 69以上で警告が出る
# 周期の数の設定
cycle = 2

## 描画処理
plot_fourier(f, cycle, N_list)

【実行結果】
黒点線が元の関数 $${f(t)}$$ です。
次数 $${N=5}$$ のフーリエ級数展開（緑色）は $${f(t)}$$ とぴったり一致します。

■ 区間 $${[-\pi, \pi)}$$ の二次関数 $${f(t) = t^2}$$
区間 $${[-\pi, \pi)}$$ の二次関数 $${f(t) = t^2}$$ が３つ並びます。

## フーリエ級数展開の設定と準備
# 対象関数 f の設定 （二次関数 t² を区間[-π, π)で繰り返し算出する関数）
f = lambda t: (t - np.pi)**2
# 次数のリストの設定
N_list = [1, 2, 5, 70]  # 72以上で警告が出る
# 周期の数の設定
cycle = 3

## 描画処理
plot_fourier(f, cycle, N_list)

【実行結果】
次数 $${N}$$ が大きくなるにつれて、元の関数に近似していく様子が見て取れます。

■ $${0,1}$$ の信号のような関数
区間 $${[-\pi, \pi)}$$ の時間 $${t}$$ を $${\pi}$$ で割った商を利用して、周期的に $${0, 1}$$ を行き来する関数にしました。

## フーリエ級数展開の設定と準備
# 対象関数 f の設定 （0,1の信号出力っぽい関数。時間を円周率πで割り算するときの商）
f = lambda t: t // np.pi
# 次数のリストの設定
N_list = [1, 5, 10, 140]  # 147以上で警告が出る
# 周期の数の設定
cycle = 3

## 描画処理
plot_fourier(f, cycle, N_list)

【実行結果】
次数 $${N=100}$$（赤色）は直線的な元の関数をかなり近似できているのではないでしょうか。

コードの「対象関数 f の設定 」の部分にさまざまな関数を設定して、フーリエ級数展開による近似の様子をお楽しみくださいませ。

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。