Pythonによる気象・気候データ解析II: スペクトル解析・EOFとSVD・統計検定と推定

の章末問題を一つずつ解いていきます。まずは第1章

A

0, 1, 3, 5, 10, 100で試した。100でほぼ矩形波となっている。

for ii in [0, 1, 3, 5, 10, 100]:
  [a, b] = fourier_square(ii)

B

解析的に求めると、

\begin{align}
a_n&=\frac{2}{T}\int_{-T}^T{f(t)\cos(\frac{2\pi nt}{T})dt} \\
b_n&=\frac{2}{T}\int_{-T}^T{f(t)\sin(\frac{2\pi nt}{T})dt} \\
今回はT&=2 \\
f(t)&=\left\{
\begin{array}{ll}
-1 & (2m-1 \leq t \leq 2m)\\
\ \ \ 1 & (2m \leq t \leq 2m+1)\\
\end{array}
\right. \quad ただし m \in \mathbf{Z}
\end{align}

これをふまえて

\begin{align}
a_n&=\int_{0}^2{f(t)\cos(n\pi t)dt}\\
&=\int_{0}^{1}{1\cos(n\pi t)dt}+\int_{1}^{2}{(-1)\cos(n \pi t)dt}\\
&=\frac{1}{n\pi}\{[\sin{(n\pi t)}]_{0}^{1}-[\sin{(n\pi t)}]_{1}^{2}\}\\
&=\frac{1}{n\pi}\{(\sin{(n\pi)}-\sin{0})-(\sin{(2n\pi)}-\sin{(\pi)})\}\\
&=0\\
b_n&=\int_{0}^2{f(t)\sin(n\pi t)dt}\\
&=\int_{0}^{1}{1\sin(n\pi t)dt}+\int_{1}^{2}{(-1)\sin(n\pi t)dt}\\
&=\frac{1}{n\pi}\{-[\cos(n\pi t)]_{0}^{1}+[\cos(n\pi t)]_{1}^{2}\}\\
&=\frac{1}{n\pi}\{-(\cos{(n\pi)}-\cos{0})+(\cos{(2n\pi)}+\cos{(n\pi)})\}\\
&=\frac{2}{n\pi}\{1-\cos{(n\pi)}\}\\
\end{align}

$a_n$は0。$b_n$の20までの解析解と数値解の差は、大きくなっているのが気になるが、0.002以下と小さな値。

[a, b] = fourier_square(21)
b_analytic = [0]
for i in range(1, 21):
  b_analytic.append(2 / (i * np.pi)*(1-np.cos(i*np.pi)))
  print(i, b_analytic[i] - b[i])
plt.plot(b_analytic[0:21] - b[0:21], 'k--')
plt.show()

1 0.00010472147772788709
2 -2.7049675121795177e-15
3 0.0003142057846043933
4 -5.358656019388808e-15
5 0.0005238142235249876
6 -8.140542039023667e-15
7 0.0007336298087972093
8 -1.0678852785049947e-14
9 0.0009437359452774963
10 -1.3383820297229466e-14
11 0.0011542165868625442
12 -1.6319543331775015e-14
13 0.0013651563972762737
14 -1.8768115819103397e-14
15 0.001576640913829211
16 -2.148590308252548e-14
17 0.001788756714859785
18 -2.4165231092219888e-14
19 0.0020015915915926746
20 -2.6767640594458966e-14

C

三角波は色々置けるが今回は図のような偶関数とした。

t = np.arange(-3, 3, 0.01)

triangle_wave = np.abs(4 * (t / 2 - np.floor(t / 2 + 0.5))) -1
plt.plot(t, triangle_wave, 'k--')
plt.hlines(y=0, xmin=-3, xmax=3, color="black", lw=0.5, linestyles="dashed")
plt.show()

0, 1, 2, 3, 5, 10（早く収束するので小さめの値）で試した。10でほぼ三角波となっている。

def fourier_triangle(n):
  a = np.zeros(n+1)
  b = np.zeros(n+1)

  for ii in range(1, n+1):
    cosine = np.cos(np.real(ii)*np.pi*t)
    sine = np.sin(np.real(ii)*np.pi*t)
    a[ii] = np.polyfit(cosine, triangle_wave, 1)[0]
    b[ii] = np.polyfit(sine, triangle_wave, 1)[0]

  triangle_approx = np.zeros_like(t)
  for ii in range(0, n+1):
    triangle_approx += a[ii]*np.cos(np.real(ii)*np.pi*t) + b[ii]*np.sin(np.real(ii)*np.pi*t)
  plt.plot(t, triangle_wave, 'k--')
  plt.plot(t, triangle_approx, 'b')
  plt.title("n = " + str(n))
  plt.show()

  return a, b

for ii in [0, 1, 2, 3, 5, 10]:
  [a, b] = fourier_triangle(ii)

D

Cで示した三角波を式で表すと、

\begin{align}
f(t)&=\left\{
\begin{array}{ll}
2t-2m-1 & (2m \leq t \leq 2m+1)\\
2m+2-2t+1 & (2m+1 \leq t \leq 2m+2)\\
\end{array}
\right. \quad ただし m \in \mathbf{Z}
\end{align}

これをふまえて

\begin{align}
a_n&=\int_{0}^2{f(t)\cos{(n\pi t)}dt}\\
&=\int_{0}^{1}{(2t-1)\cos{(n\pi t)}dt}+\int_{1}^{2}{(-2t+3)\cos{(n \pi t)}dt}\\
&=\frac{1}{n\pi}\{\int_{0}^{1}{(2t-1)\sin'{(n\pi t)} dt}+\int_{1}^{2}{(-2t+3)\sin'{(n\pi t)}dt}\} \\
&=\frac{1}{n\pi}\{[(2t-1)\sin{(n\pi t)}]_{0}^{1}-2\int_0^1{\sin{n\pi t}dt}+[(-2t+3)\sin{(n\pi t)}]_{1}^{2}+2\int_1^2{\sin{(n\pi t)}dt}\}\\
&=\frac{1}{n\pi}\{\frac{2}{n\pi}[\cos{(n\pi t)}]_0^1-\frac{2}{n\pi}[\cos{(n\pi t)}]_1^2\} \\
&=\frac{2}{n^2\pi^2}\{(\cos{n\pi}-1)-(1-\cos{n\pi})\} \\
&=\frac{4}{n^2\pi^2}(\cos{n\pi}-1) \\
b_n&=\int_{0}^2{f(t)\sin(n\pi t)dt}\\
&=\int_{0}^{1}{(2t-1)\sin{(n\pi t)}dt}+\int_{1}^{2}{(-2t+3)\sin{(n\pi t)}dt}\\
&=\frac{1}{n\pi}\{\int_0^1{(2t-1)\cos'{(n\pi t)}dt)}+\int_1^2{(-2t+3)\cos'{(n\pi t)}dt)}\} \\
&=\frac{1}{n\pi}\{[(2t-1)\cos{(n\pi t)}]_0^1-2\int_0^1{\cos{(n\pi t)}dt}+[(-2t+3)\cos{(n\pi t)}]_1^2+2\int_1^2{\cos{(n\pi t)}dt}\} \\
&=\frac{1}{n\pi}\{(\cos{n\pi}+1)-\frac{2}{n\pi}[\sin{(n\pi t)}]_0^1+(-1-\cos{n\pi})+\frac{2}{n\pi}[\sin{n\pi}]_1^2\} \\
&=0
\end{align}

$a_n$の解析結果と数値の誤差は1e-5以下とかなり小さな値。

a_analytic = [0]
for i in range(1, 21):
  a_analytic.append(4 / (i * np.pi)**2 * (np.cos(i*np.pi) - 1))
plt.plot(a_analytic[0:21] - a[0:21], 'k--')
plt.show()

b_analytic = [0] * 21
plt.plot(b_analytic[0:21] - b[0:21], 'k--')
plt.show()

Pythonによる気象・気候データ解析I 1章章末問題