コンテンツ

第９回

フーリエ級数

周期関数 $u(t)$ が「（第一種不連続点のみをもつ）区分的に連続な周期 $2T$ をもつ周期関数」であるとき、フーリエ級数は元の関数u(t)に各点収束します。
よって、 $u(t)$ を $2T$ 周期の関数とすると、

$\displaystyle u(t) &=& \frac{a_0}{2} + \sum_{n=1}^\infty \left(a_n \cos \frac{\pi nt}{T} + b_n \sin \frac{\pi nt}{T} \right),\vspace{0.5}$
$\displaystyle a_n &=& \frac{1}{T}\int_{-T}^{T} u(t) \cos \frac{\pi nt}{T} dt,\vspace{0.5}$
$\displaystyle b_n &=&\frac{1}{T}\int_{-T}^{T} u(t) \sin \frac{\pi nt}{T} dt.$

が成り立ちます。
詳細は計算数理Aの講義で説明があると思いますので、ざっくりと概要だけ述べると、どんな周期関数も三角関数の線型結合で表現できて、各三角関数の重みが $a_n$ や $b_n$ である、ということです。

無限級数をとったとき、フーリエ級数が元の関数に収束することは数学的に証明されています。
では、無限に足し合わせるのではなくて、有限個だけ足し合わせた場合は、どのような曲線が得られるでしょうか？
今回の演習では、足し合わせる項の数を増やして行ったときに、どのようにもとの関数 $u(t)$ に収束して行くのかをシミュレーションを通して見ていきましょう。

また、 $u(t)$ が不連続面を含む場合において、「滑らかな三角関数だけを足し合わせて、不連続な関数を近似する」という一見無茶なことをどのように実現しているのかを確認してみてください（課題3、4）。

課題1（ジグザグ波）

$u(t) = |t|$ on $[-1, 1]$

を周期 $2T~(T=1)$ の関数として拡張したものを考える。この関数のフーリエ級数展開は、

$u(t) = \frac{1}{2} -\frac{4}{\pi^2}\sum_{n:odd} \frac{1}{n^2} \cos \pi n t$

となる。ただし、n:oddとは $n$ が奇数であること(つまり $n=1, 3, 5, \cdots$ について和をとること）を表している。ここで、足し合わせる $n$ の最大値を $N$ として、

$u_N(t) = \frac{1}{2} -\frac{4}{\pi^2}\sum_{n:odd}^N \frac{1}{n^2} \cos \pi nt$

とおき、 $N=1, 3, 5, \cdots$ の場合における、 $u_N(t)$ のグラフを図示せよ。
その際、 $u(t)$ のグラフも重ねて描画すること。ただし、描画する $t$ の範囲は $-2$ から $2$ とせよ。（ $N$ が大きいときについて： $\cos \pi Nt$ がそこそこまともに描画できる程度に、 $t$ の刻み幅を小さくすることを意識してください）

小ネタ：関数の引数の参照渡しと値渡し

（少し難易度高めです。詳しいことは述べません。興味がある人は調べてみてください。）
関数にリスト（配列）を渡してまとめてできたらなぁ、と思う人がいるかもしれません。思わないし、別に興味ない、と言う人はこの小ネタは飛ばしてください。Pythonでは、関数の引数にリストを指定することも可能です。例えば、課題1において、以下の関数を定義すると便利かもしれません。

def a_n(n):
    return -4/(np.pi*np.pi*n*n)

こうすれば、関数a_nにnが1, 3, 5, … , Nというリストの形で入ってきても、a_n(1), a_n(3), …, a_n(N)というリストの形で返ってきます（※ 「リスト * リスト」はリストの要素ごとの掛け算。内積ではない）。

では、 $u(t)$ についても同様に関数を定義しようと思います。

def u(tau):
    for i in range(len(tau)):
        while tau[i] > T:
            tau[i] = tau[i] - 2*T
        while tmp[i] < -T:
            tau[i] = tau[i] + 2*T
    return np.abs(tau)

t = np.linspace(-2,2,M+1)
u_ans = u(t)

こうすれば、周期関数もうまく値が入ってくれそうです。しかしこの状態で横軸をリストt、縦軸をu_ansでプロットすると、tがおかしなことになっていることがわかるとおもいます（下図）。

これは、Pythonでは、引数の受け渡しの実態は、値そのものではなく、値が格納されている場所（変数が参照している場所）を渡しているからです。これを参照渡しと呼びます。（対して、値そのものを受け渡していることを値渡しと呼びます。）そのため、関数の中で変数tau（仮引数）を変更すると、同時に関数の外の変数t（実引数）も変更していることになってしまっており、グラフ描画の際にtが元々の値から変わってしまったわけです。

ちなみに、Pythonでは引数を参照渡しで渡していますが、実変数を変更される型（mutable型、リスト・辞書型など）と変更されない型（immutable型、数値（スカラー）・文字列など）があります。immutable型の場合は、きちんと代わりの変数の場所が用意され、実引数は影響を受けません（そのため見かけ上は値渡しと同じように動作します）。

では、実際にリストを引数とするときに、実引数に影響を与えないようにするにはどうすれば良いのか、ですが、以下のようにcopyモジュールのdeepcopyを使います。

import copy

def u(tau):
    tmp = copy.deepcopy(tau)
    for i in range(len(tmp)):
        while tmp[i] > T:
            tmp[i] = tmp[i] - 2*T
        while tmp[i] < -T:
            tmp[i] = tmp[i] + 2*T
    return np.abs(tmp)

こうすることで、実引数の値そのものを受け渡すことができます。単にtmp = tauとしてしまうと、また参照渡しになってしまうので、無意味になってしまいますのでご注意を。

この辺りの話は、C言語でも配列を引数にする方法があるにはありますが、そのために「ポインタとアドレス」という概念を学ぶ必要があります（計算数学演習では扱いません）。このような概念を知っていると、他の言語への応用がきき、両方の言語のより深い理解に繋がります。

課題1′

課題1 では、すでに数値積分を解析的に解いた式をプロットしたが、すでに諸君は数値積分で近似的に求める方法を知っているはずである。そこで、関数u(t)さえ定義してしまえば、自動的に $a_n$ および $b_n$ が求まるような仕様に変更せよ。さらに、横軸を $n$ とし、縦軸に $a_n$ および $b_n$ としたグラフも表示するプログラムを作成せよ。（課題1の場合、 $b_n$ は $n$ によらずほぼ0になるはずである。）

課題2（ジグザグ波の続き）

$d(t, N) = |u_N(t) - u(t) |$

とする。このとき、 $-1 \le t \le 1$ における $d(t, N)$ の最大値を $D(N)$ とおき、
$D(N)$ を画面に表示できるように、課題1のプログラムを改良せよ。
$N$ を動かして、 $D(N)$ が $0$ に収束するかどうか検証せよ。

（ $d$ は元の関数と $u_N(t)$ とののズレを表している。 $D$ はもっともずれの大きい場所の $d$ の値に相当する。）

課題3　ギブス現象

次の式（矩形波）

$u(t) = \left\{ \begin{array}{cc} 1 & (0 < t < \pi) \\ -1 & (-\pi < t < 0) \end{array} \right.$

を周期 $2\pi$ の関数として拡張したものを考え、不連続面( $t=0$ )に注目する。

$u(t)$ を有限の $N$ に対してフーリエ級数展開した関数を $u_N(t)$ とし、 $N$ の値を変化させた場合において、不連続面近くのグラフを描画せよ。

ただし、描画する範囲を、 $-0.5<t<0.5$ として、
キーボードから $N$ の入力を促し、 $u_N(t)$ と $u(t)$ を重ねて描画せよ。

参考： $N=21$ のときのグラフ

課題4

課題３の式について、不連続点以外の $t$ について、 $u_N(t)$ が $u(t)$ に収束するかどうかを確認するため、

$d(t, N) = |u_N(t) - u(t)|$

とし、 $d(0.4, N)$ が $N=1, \cdots ,100$ でどのように変化するかグラフで示せ。

チャレンジ問題

ここから”exp_data.csv”をダウンロードし、次の課題をせよ。

ダウンロードしたファイルを読み込み、1行目のデータを横軸、2行目のデータを縦軸でプロットせよ。

なお、CSVデータの読み込みは、import pandas as pdでモジュールを読み込んだのち、
```
df = pd.read_csv('exp_data.csv')
t = df.t
u = df.u
```
で行うことができる。このとき、リストtに1行目の、リストuに2行目のデータが入るが、型はこれまでのリストやArrayではなく、Seriesという型になる。とはいえ、この授業の範囲では、基本的にリストと同じように扱って良い。（データ解析を行うようになるとPandasのお世話になる。）
※ いずれのリストもサイズは10001である。
得られたデータの組 $(t, u(t))$ に対して、10次までのフーリエ級数で表し、1.のグラフに重ね描きせよ。
ここで、本来であれば $t$ は連続ではあるが、t[j+1]-t[j]が十分小さいと仮定し、近似的に定積分を行うこと。なお、周期 $T=2$ とすること。

うまくできると、以下のような結果が得られるはずである。

Masashi FUJII

計算数理A演習第9回

第９回

フーリエ級数

課題1（ジグザグ波）

小ネタ：関数の引数の参照渡しと値渡し

課題1′

課題2（ジグザグ波の続き）

課題3　ギブス現象

課題4

チャレンジ問題

Department of Mathematical and Life Sciences, Graduate School of Integrated Sciences for Life, Hiroshima University

第９回

フーリエ級数

課題1（ジグザグ波）

小ネタ：関数の引数の参照渡しと値渡し

課題1′

課題2（ジグザグ波の続き）

課題3 ギブス現象

課題4

チャレンジ問題

Department of Mathematical and Life Sciences, Graduate School of Integrated Sciences for Life, Hiroshima University

課題3　ギブス現象