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

[Pythonによる気象・気候データ解析II: スペクトル解析・EOFとSVD・統計検定と推定](https://amzn.to/3wmyUfX)
<a href="https://amzn.to/3wmyUfX" target="_blank"><img src="https://m.media-amazon.com/images/I/7180iilHkRL._SY522_.jpg" width=150/></a>
の章末問題を一つずつ解いていきます。まずは第1章

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

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

![ダウンロード.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/33eb53c5-28ca-27ef-6e10-f80cfaacaac6.png)
![ダウンロード (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/f63a336c-54cd-fa60-3c76-7d8f453f89c4.png)
![ダウンロード (2).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/1e14d31c-581d-6e8b-094a-158fd67d5e41.png)
![ダウンロード (3).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/6da3a2d4-9040-222d-c1cb-0272df9738f7.png)
![ダウンロード (4).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/a4a8d388-5415-2ccf-1103-9c37563396b7.png)
![ダウンロード (5).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/4214172a-4aa4-765f-0637-2abd0ded60f7.png)


# B
- 解析的に求めると、

```math
\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}
```

これをふまえて

```math
\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以下と小さな値。

```python
[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
```

![ダウンロード (6).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/eb9efb91-f53e-a5ee-c409-85db8ec8a9c7.png)

# C

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

```python
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()
```

![ダウンロード (7).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/88e12e34-1b0e-7ef8-5619-97553c23ceca.png)

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

```python
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
```

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

![ダウンロード (8).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/514f3112-8e42-beac-9fe5-c37ce7fc26ac.png)
![ダウンロード (9).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/684a7337-cc8e-ca73-0e87-e30d78b851c0.png)
![ダウンロード (10).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/28564e1f-e270-4fc8-3327-e2c10d05f04f.png)
![ダウンロード (11).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/aab92a2b-4616-2b5f-c1f4-eb7b64a034e9.png)
![ダウンロード (12).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/06e14661-b9fb-7244-d741-9216e1167644.png)
![ダウンロード (13).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/94259446-2ee9-5fc8-49a2-63cc522c91a2.png)

# D

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

```math
\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}
```

これをふまえて

```math
\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以下とかなり小さな値。

```python
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()
```

![ダウンロード (14).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/2b7e1e61-a3f5-aae0-c3c1-4b17fdb20fa2.png)

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

![ダウンロード (15).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/1f386b26-7812-2568-6677-f98f1c6c3424.png)

---

- [Pythonによる気象・気候データ解析I 1章章末問題](https://qiita.com/aokikenichi/items/755925f2f16e4f9a022c)