C++ - フーリエ級数展開!
Updated:
今回は、「フーリエ級数展開」を C++ で実装してみました。
ちなみに、テイラー展開は以前紹介しています。
0. 前提条件Permalink
- Linux Mint 14 Nadia (64bit) での作業を想定。
- g++ (Ubuntu/Linaro 4.7.2-2ubuntu1) 4.7.2
1. フーリエ級数展開について(簡単に)Permalink
(数式が多いので、一部 TEX で記載)
フーリエ級数展開の基本概念は、19 世紀前半にフランスの数学者フーリエ(Fourier,1764-1830)が熱伝導問題の解析の過程で考え出したものであり、「任意の周期関数は三角関数の和で表される」というものである。
さらに、与えられた関数がフーリエ級数の部分和で近似されるとき、項数をいくら増やしていっても、不連続点の近傍で誤差が生じる。これを「ギップス現象」という。
以下は1つの簡単な例。
2. C++ ソース作成Permalink
t
の範囲は −π∼π に限定している。- 計算項数は
N
の値を変更して対応する。
File: fourier_series_expansion.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
/*********************************************
* フーリエ級数展開 *
* f(t) = -1 (-pi < t <= 0 ) *
* 1 ( 0 < t <= pi) *
*********************************************/
#include <iostream> // for cout
#include <math.h> // for sin()
#include <stdio.h> // for printf()
#define N 3 // 計算項数
using namespace std;
/*
* 計算クラス
*/
class Calc
{
public:
void expandFourierSeries(); // フーリエ級数展開
private:
double calcTerm(int n, double x); //各項計算
};
/*
* フーリエ級数展開
*/
void Calc::expandFourierSeries()
{
int i; // LOOPインデックス
double t, y = 0; // 横軸、縦軸
FILE *pf; // ファイルポインタ
// 出力ファイルOPEN
pf = fopen("FourierSeriesExpansion.csv", "w");
// ヘッダ出力
fprintf(pf, "t,f(t)\n");
// 1 / 1000 刻みで計算
for (t = -M_PI; t < M_PI; t += 0.001) {
for (i = 1; i <= N; i++) y += calcTerm(i, t);
fprintf(pf, "%lf,%lf\n", t, 4 / M_PI * y);
y=0;
}
// 出力ファイルCLOSE
fclose(pf);
}
/*
* 各項計算
*/
double Calc::calcTerm(int n, double t)
{
return sin((2 * n - 1) * t) / (2 * n - 1);
}
/*
* メイン処理
*/
int main()
{
try
{
// 計算クラスインスタンス化
Calc objCalc;
// フーリエ級数展開
objCalc.expandFourierSeries();
}
catch (...) {
cout << "例外発生!" << endl;
return -1;
}
// 正常終了
return 0;
}
3. C++ ソースコンパイルPermalink
(-Wall
は警告出力、-O2
最適化のオプション)
$ g++ -Wall -O2 -o fourier_series_expansion fourier_series_expansion.cpp
何も出力されなければ成功。
4. 実行Permalink
$ ./fourier_series_xxpansion
コンソールには特に何も表示しない。
アプリと同じディレクトリに FourierSeriesExpansion.csv
という CSV ファイルが作成される。
内容は以下のようになっているはず。
File: FourierSeriesExpansion.csv
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
t,f(t)
-3.141593,-0.000000
-3.140593,-0.003820
-3.139593,-0.007639
-3.138593,-0.011459
-3.137593,-0.015278
-3.136593,-0.019098
-3.135593,-0.022917
-3.134593,-0.026735
-3.133593,-0.030554
-3.132593,-0.034372
-3.131593,-0.038190
-3.130593,-0.042007
-3.129593,-0.045824
-3.128593,-0.049640
-3.127593,-0.053456
-3.126593,-0.057271
-3.125593,-0.061085
-3.124593,-0.064899
-3.123593,-0.068712
:
====< 途中省略 >====
:
3.122407,0.073230
3.123407,0.069418
3.124407,0.065605
3.125407,0.061792
3.126407,0.057978
3.127407,0.054163
3.128407,0.050347
3.129407,0.046531
3.130407,0.042714
3.131407,0.038897
3.132407,0.035080
3.133407,0.031261
3.134407,0.027443
3.135407,0.023624
3.136407,0.019805
3.137407,0.015986
3.138407,0.012167
3.139407,0.008347
3.140407,0.004528
3.141407,0.000708
5. グラフ化Permalink
数字だけを眺めてもよく分からないので、R でグラフ化(プロット)してみた。
【元の関数グラフ】
以下、計算項数を 1, 2, 3, 5, 10, 20, 50, 100, 200, 500, 1000, 10000, 100000 個として計算した結果をグラフ化したもの。
項数を増やすにつれて元の関数のグラフに近付いていくのがよく分かる。
電気工学、音響学、振動、光学等でよく使用する重要な概念です。応用範囲は広いので他にも利用できるかと思います。
以上。
Comments