【入門】フーリエ級数(周期2L)Python【数値計算】

【入門】フーリエ級数(周期2L)Python【数値計算】 数値計算
【入門】フーリエ級数(周期2L)Python【数値計算】

MATLAB、Python、Scilab、Julia比較ページはこちら
https://www.simulationroom999.com/blog/comparison-of-matlab-python-scilab/

はじめに

の、

MATLAB,Python,Scilab,Julia比較 第5章 その52【フーリエ級数(周期2L)⑤】

を書き直したもの。

フーリエ級数、フーリエ係数を周期2πを任意周期2Lに改造した。
これをPythonで実現する。

【再掲】フーリエ係数を求めるプログラムフロー

まずは、プログラムフローを再掲。

  • csvファイル読み込み
  • 各種変数初期化
  • フーリエ係数算出
  • n=10,50,200のパターンでフーリエ級数で波形を合成
  • グラフにプロット

今回はフーリエ係数を求めるプログラムをPythonで実現。

【再掲】任意周期フーリエ級数、フーリエ係数

実現する数式はこれ

\(
\begin{eqnarray}
\displaystyle a_n&=&\frac{1}{L}\int_{-L}^{L}f(x)\cos(\frac{n\pi x}{L})dx\\
\displaystyle b_n&=&\frac{1}{L}\int_{-L}^{L}f(x)\sin(\frac{n\pi x}{L})dx\\
\displaystyle a_0&=&\frac{1}{L}\int_{-L}^{L}f(x)dx\
\end{eqnarray}
\)

Pythonコード

まず、使用する波形を取り込んだcsvファイル

Pythonコードは以下になる。

import numpy as np
import matplotlib.pyplot as plt

N=1000  # 係数算出項数(同定元波形のplotよりも少なく)
L=10    # 周期/2
wave=np.loadtxt('wave.csv',delimiter=',') # 同定波形読み込み

points=len(wave) # 波形のplot数取得
fx=np.array(wave) # 波形を行ベクトルへ
dx=2*L/points # 1plotあたりのx軸幅
x=np.linspace(-L,L,points); # -L~+Lの範囲で波形plot数分の等差数列

a = np.zeros(N) # a係数群格納用
b = np.zeros(N) # b係数群格納用

for n in range(1,N+1):
    # 係数a_n算出
    # a_n = (1/L)∫f(fx)cos(nxπ/L)dx
    a[n-1] = fx@np.cos(n*np.pi*x/L).T*dx/L;
    
    # 係数b_n算出
    # a_n = (1/L)∫f(fx)cos(nxπ/L)dx
    b[n-1] = fx@np.sin(n*np.pi*x/L).T*dx/L;

# 係数a_0算出
a0=np.sum(fx)*dx/L

Ns = [10,50,200]
fig = plt.figure()

for i in range(0,len(Ns)):
    NN = Ns[i]; # 今回のa_n,b_n項数

    # f(x)=a_0/2+Σ(a_n cos(nxπ/L)+ b_n sin(nxπ/L))
    Fourier_series=np.ones(points)*a0/2
    for n in range(1,NN+1):
        Fourier_series = Fourier_series+(a[n-1]*np.cos(n*np.pi*x/L)+b[n-1]*np.sin(n*np.pi*x/L))
    
    ax = fig.add_subplot(len(Ns), 1, i+1)
    
    # 元波形とフーリエ級数波形の表示
    ax.plot(x, fx,'-b',lw=3)
    ax.plot(x, Fourier_series,'-r',lw=2)
    ax.set_title('n={:d}'.format(NN))
    ax.set_ylim(-0.1,1.1)
    ax.set_xlim(-L,L)
    ax.grid(linestyle='dotted')

plt.show()

処理結果

処理結果は以下。

フーリエ係数算出任意周期(Python)

\(-\pi\sim\pi\)だけでなく、任意の周期を設定できることがわかる。
今回の場合は\(-10\sim10\)の範囲にしているが、当然変更することも可能。

まとめ

  • フーリエ級数、フーリエ係数の任意周期版のプログラムをPythonで作成。
  • -π~πだけでなく、-10~10のような任意の周期に適応可能。

MATLAB、Python、Scilab、Julia比較ページはこちら

マンガでわかるフーリエ解析

https://amzn.to/4azyPDj

手を動かしてまなぶ フーリエ解析・ラプラス変換

https://amzn.to/4bDH3vu

物理数学 量子力学のためのフーリエ解析・特殊関数

https://amzn.to/4bsygfU

単位が取れるフーリエ解析ノート

https://amzn.to/3V83fIl

今日から使えるフーリエ変換 普及版 式の意味を理解し、使いこなす

https://amzn.to/3ysbfvf

コメント

タイトルとURLをコピーしました