Julia(SymPy)でフーリエ級数を計算する方法

どうも、木村(@kimu3_slime)です。

今回は、Julia(SymPy)でフーリエ級数を計算する方法を紹介します。

 

準備

SymPy, Plotsを使うので、持っていなければインストールしておきましょう。

準備として、以下のコードを実行しておきます。

 

フーリエ級数を計算する方法

はじめに、記号として使う変数を用意しておきます。

 

区間\([-\pi, \pi]\)における関数\(f\)のフーリエ級数展開は、

\[ \begin{aligned}f(x)=a_0 + \sum_{n=1}^\infty (a_n \cos nx +b_n \sin nx)\end{aligned} \]

\[ \begin{aligned}a_0=\frac{1}{2\pi} \int_{-\pi}^{\pi} f(x)dx\end{aligned} \]

\[ \begin{aligned}a_n=\frac{1}{\pi} \int_{-\pi}^{\pi} f(x)\cos nxdx\end{aligned} \]

\[ \begin{aligned}b_n=\frac{1}{\pi} \int_{-\pi}^{\pi} f(x)\sin nx dx\end{aligned} \]

定義されます

「sympy.fourier_series(関数,(x,範囲))」で特定の範囲における関数のフーリエ級数展開が求められます。

\[ \begin{aligned}2 \sin{\left(x \right)} – \sin{\left(2 x \right)} + \frac{2 \sin{\left(3 x \right)}}{3} + \ldots\end{aligned} \]

結果は無限級数であり、部分的に表示されています。

最初の0でない項を切り出して求めるには、「数式.truncate(n= 項数)」を使いましょう。

\[ \begin{aligned}2 \sin{\left(x \right)} – \sin{\left(2 x \right)} + \frac{2 \sin{\left(3 x \right)}}{3}\\ – \frac{\sin{\left(4 x \right)}}{2} + \frac{2 \sin{\left(5 x \right)}}{5}\end{aligned} \]

「plot(関数)」を使えば、この結果のグラフを描くことができます。

 

矩形波(ヘビサイド関数)について、同じことをやってみましょう。

\[ \begin{aligned}\frac{2 \sin{\left(x \right)}}{\pi} + \frac{2 \sin{\left(3 x \right)}}{3 \pi} \\+ \frac{2 \sin{\left(5 x \right)}}{5 \pi} + \frac{2 \sin{\left(7 x \right)}}{7 \pi} + \frac{1}{2}\end{aligned} \]

ただし、その最初の第1項(定数関数)をプロットしようとすると、エラーが表示されていまいます。

\[ \begin{aligned}\frac{1}{2}\end{aligned} \]

Symの定数は、変数\(x\)を含んでいないので、plotの対象になってくれないようです。

定数関数を作り出すために、ヘビサイド関数の足し合わせ「Heaviside(x)+Heaviside(-x)」を利用します。最適な方法ではない気がしますが、これを使えば定数関数として扱ってくれます。

 

以上の結果をまとめて、与えられた関数のフーリエ級数展開のグラフを重ねて描く関数を作ってみましょう。

最初に「p」でプロットの土台を作り、その後「plot!」でグラフを重ねて書いていきます。最初は関数そのもの、続いてフーリエ級数展開1項、3項、5項、100項の順です。1項のみのケースが少し汚い式ですが、さきほど見た問題への対処をしています。

 

のこぎり波のフーリエ級数展開は次の通り。

 

矩形波のフーリエ級数展開は次の通り。

 

絶対サイン波のフーリエ級数展開は次の通り。

\[ \begin{aligned}- \frac{4 \cos{\left(2 x \right)}}{3 \pi} – \frac{4 \cos{\left(4 x \right)}}{15 \pi} \\- \frac{4 \cos{\left(6 x \right)}}{35 \pi} – \frac{4 \cos{\left(8 x \right)}}{63 \pi} + \frac{2}{\pi}\end{aligned} \]

 

\(x^2\)のフーリエ級数展開は次の通り。

\[ \begin{aligned}- 4 \cos{\left(x \right)} + \cos{\left(2 x \right)} – \frac{4 \cos{\left(3 x \right)}}{9}\\ + \frac{\cos{\left(4 x \right)}}{4} + \frac{\pi^{2}}{3}\end{aligned} \]

 

以上、Julia(SymPy)でフーリエ級数を計算する方法を紹介してきました。

以前Juliaで級数をプロットする方法で、フーリエ級数のプロットも紹介しましたが、そちらはフーリエ係数を手計算で導いた上で与えています。

SymPyを使うと、フーリエ係数の計算も込みで解いてくれるので、簡単ですね。

木村すらいむ(@kimu3_slime)でした。ではでは。

 

1から始める Juliaプログラミング
進藤 裕之(著), 佐藤 建太(著)
コロナ社 (2020-03-26T00:00:01Z)
5つ星のうち4.5
¥7,353 (コレクター商品)

 

こちらもおすすめ

フーリエ係数に関するオイラーの公式の導出

矩形波とは:フーリエ級数展開の求め方、ギブス現象、ライプニッツの級数

のこぎり波とは:フーリエ級数展開の求め方

絶対サイン波|sin x|とは:フーリエ級数展開の求め方

単位ステップ関数とは、ラプラス変換と微分方程式への応用

Juliaで級数を計算し、グラフを描く方法

Julia(SymPy)で数列の和、無限級数、べき級数を求める方法

Julia(SymPy)によるテイラー級数展開の求め方(指数対数、三角、双曲線)