next up previous
: 6 謝辞 : DCPAM3 第2部 離散化 : 4 参考文献


5 座標系・変換公式

9% latex2html id marker 7487
\setcounter{footnote}{9}\fnsymbol{footnote} 99 9% latex2html id marker 7488
\setcounter{footnote}{9}\fnsymbol{footnote} 99

A. 球面調和函数

ここでは連続系での球面調和函数を定義し, スペクトル計算の理解に必要な性質を挙げ, 証明する.

まず球面調和函数を定義し, 次いで球面調和函数が完全直交系をなすことを主張する. このことにより, 球面上に分布するあらゆる連続関数が 球面調和函数の重ね合わせで一意的に表されることになる.

球面調和函数は2次元ラプラシアンに関する固有関数であり, このために全波数という概念が生まれる. 参考までにこのことも記しておく.

さらに, 球面調和函数を空間微分した結果も書いておく.

  1. 定義と性質(球面調和函数, Legendre函数, Legendre陪函数)
  2. 空間微分
  3. 全波数の概念

また, イメージをつかむために, ルジャンドル(陪)関数のグラフを示す.

A..1 定義と性質

ここでは, 岩波公式集10の Legendre函数・陪函数 $\tilde{P}_n^m$ , 2 で規格化したLegendre函数・陪函数 $P_n^m$, $4 \pi$ で規格化した球面調和函数 $Y_n^m$ の順に定義する. さらにそれらの性質として, 従う微分方程式, 漸下式, 完全規格直交性について述べる.

A..1.1 岩波公式集の Legendre函数・陪函数$\tilde{P}_n^m$

A..1.2 2 で規格化した Legendre函数・陪函数$P_n^m$

A..1.3 球面調和函数$Y_n^m$

A..2 球面調和函数の空間微分

ここでは, 球面調和函数 $Y_n^m(\phi,\lambda)$

の計算をする.

A..2.1 x微分


$\displaystyle \frac{1}{r \cos \phi} \DP{Y_n^m}{\lambda}
= \frac{1}{r \cos \phi}...
...im \lambda) \right)
= \frac{im}{r \cos \phi}
P_n^m (\sin \phi) \exp(im \lambda)$     (97)

A..2.2 y微分


$\displaystyle \frac{1}{r} \DP{Y_n^m}{\phi}
= \frac{1}{r}
\DP{}{\phi}
\left( P_n...
...) \right)
= \frac{\sqrt{1-\mu^2} }{r}
\DD{}{\mu} P_n^{m} (\mu) \exp(im \lambda)$     (98)

A..2.3 2次元ラプラシアン


$\displaystyle \nabla_H^2
Y_n^m$ $\textstyle \equiv$ $\displaystyle \frac{1}{r^2}
\left[ \DP{}{\mu} \left( (1-\mu^2) \DP{}{\mu} \right)
+ \frac{1}{1-\mu^2} \DP[2]{}{\lambda}
\right] Y_n^m$ (99)
  $\textstyle =$ $\displaystyle \frac{1}{r^2}
\left[ \frac{1}{\cos \phi}
\DP{}{\phi} \left( \cos \phi \DP{}{\phi} \right)
+ \frac{1}{\cos^2 \phi} \DP[2]{}{\lambda}
\right] Y_n^m$ (100)
  $\textstyle =$ $\displaystyle - \frac{n(n+1)}{r^2} Y_n^m$ (101)

A..3 コメント -- 全波数について

球面調和函数 $Y^m_n(\lambda, \phi)$ において $n$ のことを全波数と呼ぶ.

全波数には, 座標系の回転に関して不変である, という特徴がある. すなわち, 任意の $Y_n^m(\lambda,\phi)$ は 回転して得られる座標系 $(\lambda', \phi')$ における 全波数 $n$ の球面調和函数 $\{ Y_n^m(\lambda', \phi') \vert m=-n,-n+1, \cdots, n \} $ の和で表現できる :

$\displaystyle Y^m_n(\lambda, \phi)
= \sum_{m'=-n}^{n} A_n^{m'} Y^{m'*}_n(\lambda',\phi')$     (102)

のである12. この特徴は, 球面調和函数が2次元ラプラシアンの 固有値であることによっている13.

A..4 グラフ

$P_n^m(\mu)$ の概形をつかむために, 2で規格化した $P_n,P_n^1,P_n^2$ 14 のグラフを示す.

\Depsf[10cm]{zahyou/lg1.ps}

岩波公式集の Legendre函数 $\tilde{P}_n$ のグラフ (森口, 宇田川, 一松, 1960)

\Depsf[][7cm]{zahyou/lg2.ps}

Legendre函数 $ \overline{P_n^1} = P_n^1/\sqrt{2},
\overline{P_n^2} = P_n^2/\sqrt{2}$ のグラフ (森口, 宇田川, 一松, 1960)
14 1414

B. 微分公式, GCMの変数の微分関係式

15

ここでは, スカラー量, ベクトルの微分を計算する. さらにそれらを元に, 発散$D$, 渦度 $\zeta$, 速度ポテンシャル$\chi$, 流線関数$\psi$$(U,V)$ との関係を付ける.

B..1 スカラー量の微分

スカラー量 $f(\lambda,\phi)$$x$ 微分は ${\displaystyle \frac{1}{r \cos \phi} \DP{f}{\lambda} }$ で与えられる.

$f$$y$ 微分は ${\displaystyle \frac{1}{r} \DP{f}{\phi}
\left( = \frac{\cos \phi}{r} \DP{f}{\mu} \right) }$ で与えられる.

$f$ の2次元ラプラシアンは

$\displaystyle \nabla^2_H f$ $\textstyle \equiv$ $\displaystyle \frac{1}{r^2}
\left[ \frac{1}{\cos \phi}
\DP{}{\phi} \left( \cos \phi \DP{}{\phi} \right)
+ \frac{1}{\cos^2 \phi} \DP[2]{}{\lambda}
\right] f$ (103)
  $\textstyle =$ $\displaystyle \frac{1}{r^2}
\left[
\DP{}{\mu}
\left\{ (1-\mu^2) \DP{}{\mu} \right\}
+ \frac{1}{1-\mu^2} \DP[2]{}{\lambda}
\right] f$ (104)

で与えられる.

B..2 ベクトル量の微分

2次元ベクトル場 $\Dvect{v}=(v_1,v_2)$ の水平発散は


$\displaystyle \mbox{div}_H \Dvect{v}$ $\textstyle \equiv$ $\displaystyle \frac{1}{r \cos \phi} \DP{v_1}{\lambda}
+ \frac{1}{r \cos \phi} \DP{}{\phi} (v_2 \cos \phi)$ (105)
  $\textstyle =$ $\displaystyle \frac{1}{r \sqrt{1-\mu^2} } \DP{v_1}{\lambda}
+ \frac{1}{r} \DP{}{\mu}( \sqrt{1-\mu^2} v_2 )$ (106)

で与えられる.



$\Dvect{v}$ の回転の $r$ 成分は,

$\displaystyle (\mbox{rot} \Dvect{v})_r$ $\textstyle \equiv$ $\displaystyle \frac{1}{r \cos \phi} \DP{v_2}{\lambda}
- \frac{1}{r \cos \phi} \DP{}{\phi}(v_1 \cos \phi)$ (107)
  $\textstyle =$ $\displaystyle \frac{1}{r \sqrt{1-\mu^2} } \DP{v_2}{\lambda}
- \frac{1}{r} \DP{}{\mu}(\sqrt{1-\mu^2}v_1)$ (108)

で与えられる.




以上で得られた微分公式を元に, 以下に実際にGCMで使用する便利な微分の公式を並べておく.

B..3 発散

水平分布する速度場

$\displaystyle (u,v) \equiv
\left( \frac{U}{\cos \phi}, \frac{V}{\cos \phi} \right)$     (109)

の水平発散 $D$ を, $U,V$ を用いて表す.


$\displaystyle D$ $\textstyle =$ $\displaystyle \frac{1}{r \cos \phi} \DP{u}{\lambda}
+ \frac{1}{r \cos \phi} \DP{}{\phi} (v \cos \phi)$ (110)
  $\textstyle =$ $\displaystyle \frac{1}{r \cos^2 \phi} \DP{U}{\lambda}
+ \frac{1}{r \cos \phi} \DP{V}{\phi}$ (111)
  $\textstyle =$ $\displaystyle \frac{1}{r (1-\mu^2)} \DP{U}{\lambda}
+ \frac{1}{r} \DP{V}{\mu}$ (112)

B..4 渦度

水平分布する速度場

$\displaystyle (u,v) = \left( \frac{U}{\cos \phi}, \frac{V}{\cos \phi} \right)$     (113)

の渦度 $\zeta$ を, $U,V$ を用いて表す.


$\displaystyle \zeta$ $\textstyle =$ $\displaystyle \frac{1}{r \cos \phi} \DP{v}{\lambda}
- \frac{1}{r \cos \phi} \DP{}{\phi}(u \cos \phi)$ (114)
  $\textstyle =$ $\displaystyle \frac{1}{r \cos^2 \phi} \DP{V}{\lambda}
- \frac{1}{r \cos \phi} \DP{U}{\phi}$ (115)
  $\textstyle =$ $\displaystyle \frac{1}{r (1-\mu^2)} \DP{V}{\lambda}
- \frac{1}{r} \DP{U}{\mu}$ (116)

B..5 速度ポテンシャル, 流線関数と $(U,V)$

速度ポテンシャル $\chi$, 流線関数 $\psi$

$\displaystyle D$ $\textstyle \equiv$ $\displaystyle \nabla_H^2 \chi$ (117)
$\displaystyle \zeta$ $\textstyle \equiv$ $\displaystyle \nabla_H^2 \psi$ (118)

で定義される. $(U,V)$$\chi ,\psi $ で表す.


$\displaystyle U$ $\textstyle =$ $\displaystyle - \frac{1-\mu^2}{r} \DP{\psi}{\mu}
+ \frac{1}{r} \DP{\chi}{\lambda}$ (119)
$\displaystyle V$ $\textstyle =$ $\displaystyle \frac{1}{r} \DP{\psi}{\lambda}
+ \frac{1-\mu^2}{r} \DP{\chi}{\mu}$ (120)

となる.

15 1515 15 1515

C. Legendre函数 $P_n$ の性質

ここでは Legendre函数 $P_n$の性質である

  1. $ n-1$ 次以下の多項式との積を $-1 \leq \mu \leq 1$ まで積分すると 零になること
  2. $ P_n(\mu)$$-1 < \mu < 1$$n$ 個の零点を持つこと,

を記す. 1 より Gauss 格子を定義することが保証される. また, 1, 2 は共に Gauss-Legendre の公式の証明に用いられる.

C..1 多項式とLegendre函数の積の積分

$ P_n(\mu)$ は, $\mu$$n$ 次多項式である. $ n-1$ 次以下の任意の多項式は $P_0 \sim P_{n-1}$ の和で表されること, $P_n$ の直交性から明らかに, $ n-1$ 次以下の任意の多項式 $f(\mu)$ との積を積分すると

$\displaystyle \int_{-1}^1 f(\mu) P_n(\mu) d \mu = 0$     (121)

が成り立つ ことがわかる.

C..2 Legendre函数の零点

$P_n$$-1 < \mu < 1$$n$ 個の互いに異なる零点を持っている. このことについて, 以下に証明しておく. (寺沢, 1983 の10.7 節より)

  1. $f(x)=(x-1)^n (x+1)^n$ を導入する.
  2. $f=0$ の解は, $x=-1,1$ である. ゆえに, Rolle の定理により, $f'$ は ある $\alpha$$-1<\alpha<1$ ) で $f'(\alpha)=0$ となる.
    $f'=2nx(x^2-1)^{n-1}$ より, $f'=0$ の解は $x=-1,\alpha,1$ のみである.
  3. 同様に, $f''=0$ の解は $x=-1,\beta_1,\beta_2,1$ $-1 < \beta_1 < \beta_2 < 1$)のみ.
  4. 以上を繰り返すと, $f^{(n)}=0$ の解は $-1$$1$ の間で互いに異なる $n$ 個の 解を持つ. ($x=-1,1$ は解でないことに注意せよ. )
  5. したがって, ${\displaystyle
P_n= \frac{1}{2^n n!}
\DD[n]{}{\mu} (\mu^2-1)^n }$$-1$$1$ の間で互いに異なる $n$ 個の解を持つ. (証明終り)



この零点の求め方としては, ${\displaystyle x_j=\cos \frac{j-1/2}{n}\pi }$ を近似解として Newton 法を用いるという方法がある.

15 1515 15 1515

D. 積分評価

D..1 Gauss の台形公式

ここでは Gauss の台形公式を示す.

波数 $M$ 以下の三角函数で表現される $g(\lambda)$ $0 \le \lambda < 2 \pi$

$\displaystyle g(\lambda) = \sum_{m=-M}^{m=M} g_m \exp( i m \lambda )$     (122)

について $M < I$ を満たすように $I$ をとると,
    $\displaystyle \frac{1}{2\pi} \int_0^{2\pi} g(\lambda) d \lambda
= \frac{1}{I} \sum_{n=1}^{I} g(\lambda_n)$ (123)
    $\displaystyle \lambda_n = \frac{2\pi (n-1)}{I} \ (n=1,2,\cdots,I)$ (124)

が成り立つ. これを Gauss の台形公式という.

より実用的な公式は,

    $\displaystyle \sum_{n=1}^{I} \exp(i m \lambda_n)
= \left\{
\begin{array}{ll}
I & \ \ (m=0) \\
0 & \ \ (0 < \vert m\vert < I) \\
\end{array}\right.$ (125)
    $\displaystyle \lambda_n = \frac{2\pi (n-1)}{I} \ (n=1,2,\cdots,I)$ (126)

である. この証明は, $I > M$ ($\vert m\vert$ の最大値) より $m \neq 0$ の時には ${\displaystyle \exp(i m \lambda_n)
= \exp \left( \frac{2\pi i m (n-1)}{I} \right) }$ において, 全ての $n$ について $m(n-1)$$I$ の整数倍になることがないことを考慮すると明らかである ( $m$, $ n-1$ はともに $I$ よりも小さい整数なので, $m(n-1)$$I$ の整数倍にならない) 16.

以下に Gauss の台形公式の証明を記す.
まず, 左辺を計算すると,

$\displaystyle \frac{1}{2\pi} \int_0^{2\pi} g(\lambda) d \lambda$ $\textstyle =$ $\displaystyle \sum_{m=-M}^{M} \frac{1}{2\pi} g_m
\int_0^{2\pi} \exp(i m \lambda) d \lambda
=g_0$ (128)

である. ここで, $\int_0^{2\pi} \exp(i m \lambda) d \lambda $$m=0$ の項しか残らないことを使った. 一方右辺は
$\displaystyle \frac{1}{I} \sum_{n=1}^{I} g(\lambda_n)$ $\textstyle =$ $\displaystyle \frac{1}{I}
\sum_{n=1}^{I} \sum_{m=-M}^{M}
g_m \exp(i m \lambda_n)$ (129)
  $\textstyle =$ $\displaystyle g_0
+ \sum_{m=-M, m\neq0}^{M} \frac{g_m}{I}
\sum_{n=1}^{I}
\left( \exp(\frac{2 \pi i m}{I}) \right)^{n-1}$ (130)

ここで, 上に示した「より実用的な公式」により
$\displaystyle \sum_{n=1}^{I}
\left( \exp(\frac{2 \pi i m}{I}) \right)^{n-1} =0
\ \ \ \ (m \neq 0)$     (131)

が成り立つ. したがって,
$\displaystyle \frac{1}{2\pi} \int_0^{2\pi} g(\lambda) d \lambda
= \frac{1}{I} \sum_{n=1}^{I} g(\lambda_n)$     (132)

となる.

D..2 Gauss-Legendreの公式

$f(\mu)$$2J-1$ 次以下の多項式とする. $P_n$ を2で規格化した n 次の Legendre函数とする. このとき, ${\displaystyle \int_{-1}^1 f d \mu }$$P_J$ の零点である Gauss 格子 $\mu_j(j=1,2,\cdots,J)$ における $f$ の値 $f(\mu_j)$ のみを用いて, 次式にもとづいて正確に評価することができる.

    $\displaystyle \int^{1}_{-1} f(\mu) d \mu = 2 \sum_{j=1}^{J} f(\mu_j) w_j$ (133)
    $\displaystyle w_j = \frac{1}{2} \int_{-1}^{1}
\frac{P_J(\mu)}{(\mu-\mu_j) P^{'}_J (\mu_i)}d \mu
= \frac{(2J-1)(1-\mu_j^2)}{(J P_{J-1}(\mu_j))^2 } .$ (134)

ここで, $w_j$ は Gauss 荷重と呼ばれる.


以下では上の式を証明する. ただし, Legendre函数としては, 最初は岩波公式集のLegendre函数 $\tilde{P_n}$ を用い, 最後に2で規格化したLegendre函数 $P_n$ に 直すことにする17.

STEP 1 Lagrange 補間の導入

$f(\mu)$$K$ 次多項式( $0 \le K \le 2J-1$)とする. $\tilde{P}_n$ を岩波公式集のLegendre函数(Rodriguesの公式) とする.


$\displaystyle \int_{-1}^1 \tilde{P}_n(\mu) \tilde{P}_{n'}(\mu) d \mu
= \frac{2}{2n+1} \delta_{nn'}$     (135)

$L(\mu)$ を, $f(\mu_j)$ を Lagrange 補間公式にしたがって補間した多項式として 定義する.


$\displaystyle L(\mu) \equiv
\sum_{j=1}^{J} f(\mu_j)
\prod_{k=1,k \neq j}^{J}
\frac{ \mu-\mu_k}{\mu_j-\mu_k}$     (136)

このとき, 各$j$ について $ L(\mu_j) = f(\mu_j) $ である. ここで $L$ は, $0 \le K \le J-1$ の時($f$$J-1$ 次以下の多項式)のときは 厳密に $L=f$ になる18 ことに注意せよ.

したがって, 関数 $f(\mu)-L(\mu)$

$f(\mu)-L(\mu)$$\mu$ について $-1$ から 1 まで積分する. $J \le K \le 2J-1$ の時については Legendre函数の直交性より, $\tilde{P}_J(\mu) S(\mu)$ の積分は零である. したがって,

$\displaystyle \int_{-1}^{1} f(\mu) d \mu$ $\textstyle =$ $\displaystyle \int_{-1}^1 L(\mu) d \mu$ (138)
  $\textstyle =$ $\displaystyle \sum_{j=1}^{J} f(\mu_j)
\int^1_{-1}
\frac{ {\displaystyle
\prod_{...
...
}
{ (\mu-\mu_j)
{\displaystyle
\prod_{k=1,k \neq j}^{J}(\mu_j-\mu_k) }
}
d \mu$ (139)
  $\textstyle =$ $\displaystyle \sum_{j=1}^{J} f(\mu_j)
\int_{-1}^1
\frac{\tilde{P}_J(\mu)}
{(\mu-\mu_j)\tilde{P}^{'}_J(\mu_j)} d \mu$ (140)
  $\textstyle =$ $\displaystyle 2 \sum_{j=1}^{J} f(\mu_j) w_j$ (141)

ここで, 証明すべき式の $P_J$ は規格化されていて, 上の式の $\tilde{P}_J$ は規格化されていないのにもかかわらず 同じ $w_j$ が使われているが, $\tilde{P}_J$$P_J$ の規格化定数は同じなので consistent である.

STEP 2 ${\displaystyle w_j
= \frac{1}{2} \int_{-1}^1
\frac{\tilde{P}_J(\mu)}
{(\mu-\mu_j)\tilde{P}^{'}_J(\mu_j)} d \mu
}$ の漸化式を用いた変形

漸化式 (岩波の Lgendre 関数・陪関数の従う漸化式) において $m=0$ とした式

$\displaystyle (n+1) \tilde{P}_{n+1}(\mu)
= (2n+1) \mu \tilde{P}_n(\mu) - n \tilde{P}_{n-1}(\mu)
\ \ \ \ (n=0,1,2,\cdots)$     (142)

より,
$\displaystyle (n+1) \left\vert
\begin{array}{ll}
\tilde{P}_{n+1}(x) & \tilde{P}_n(x) \\
\tilde{P}_{n+1}(y) & \tilde{P}_n(y)
\end{array}\right\vert$ $\textstyle =$ $\displaystyle \left\vert
\begin{array}{ll}
(2n+1)x \tilde{P}_n(x)-n\tilde{P}_{n...
...+1)y \tilde{P}_n(y)-n\tilde{P}_{n-1}(y)
& \tilde{P}_n(y)
\end{array}\right\vert$ (143)
  $\textstyle =$ $\displaystyle (2n+1)(x-y)\tilde{P}_n(x)\tilde{P}_n(y)$ (144)
    $\displaystyle +n (- \tilde{P}_{n-1}(x) \tilde{P}_n(y)
+ \tilde{P}_{n-1}(y) \tilde{P}_n(x) )$ (145)
  $\textstyle =$ $\displaystyle (2n+1)(x-y)\tilde{P}_n(x)\tilde{P}_n(y)
+
n \left\vert
\begin{arr...
...de{P}_{n-1}(x) \\
\tilde{P}_{n}(y) & \tilde{P}_{n-1}(y)
\end{array}\right\vert$ (146)

となる. この式を $n=0,1,\cdots,n-1$ について加えると,
$\displaystyle n \left\vert
\begin{array}{ll}
\tilde{P}_n(x) & \tilde{P}_{n-1}(x) \\
\tilde{P}_n(y) & \tilde{P}_{n-1}(y)
\end{array}\right\vert$ $\textstyle =$ $\displaystyle \sum_{k=0}^{n-1} (2k+1)(x-y)\tilde{P}_k(x)\tilde{P}_k(y)$ (147)

が成り立つ. ここで $n=J,x=\mu,y=\mu_j$ とすると $\tilde{P}_J(\mu_j)=0$ より,
$\displaystyle J \tilde{P}_J (\mu) \tilde{P}_{J-1} (\mu_j)$ $\textstyle =$ $\displaystyle \sum_{k=0}^{J-1} (2k+1) (\mu-\mu_j)
\tilde{P}_k(\mu) \tilde{P}_k(\mu_j).$ (148)

よって,
$\displaystyle \frac{\tilde{P}_J (\mu) }{\mu-\mu_j}$ $\textstyle =$ $\displaystyle \frac{
{\displaystyle
\sum_{k=0}^{J-1} (2k+1)
\tilde{P}_k(\mu) \tilde{P}_k(\mu_j) }
}
{J \tilde{P}_{J-1}(\mu_j)}$ (149)

である. したがって,
$\displaystyle w_j$ $\textstyle =$ $\displaystyle \frac{1}{2} \int_{-1}^1
\frac{\tilde{P}_J(\mu)}
{(\mu-\mu_j) \tilde{P}^{'}_J(\mu_j)}
d \mu$ (150)
  $\textstyle =$ $\displaystyle \frac{1}{2 J \tilde{P}_{J-1}(\mu_j)
\tilde{P}^{'}_{J} (\mu_j)}
\sum_{k=0}^{J-1}
(2k+1) \tilde{P}_k(\mu_j)
\int_{-1}^1 \tilde{P}_k(\mu) d \mu$ (151)
  $\textstyle =$ $\displaystyle \frac{1}{J \tilde{P}_{J-1}(\mu_j)
\tilde{P}^{'}_{J} (\mu_j)}$ (152)

である. ただし, (A.80) における積分は, k=0 の時のみ0でない値を 持つこと, および $\tilde{P}_0 = 1$ を使った. さらに, 漸化式
$\displaystyle (1-\mu^2) \DP{\tilde{P}_n}{\mu}
= n \tilde{P}_{n-1}(\mu) - n \mu \tilde{P}_n(\mu)$     (153)

$n=J, \mu=\mu_j$ とする. $\tilde{P}_J(\mu_j)=0$ より,
$\displaystyle w_j$ $\textstyle =$ $\displaystyle \frac{1-\mu_j^2}{(J \tilde{P}_{J-1}(\mu_j))^2 }$ (154)

となる.

STEP3 $\tilde{P}_n$ の規格化

$P_n$

$\displaystyle \int_{-1}^1 P_n(\mu) P_n'(\mu) d \mu = 2$     (155)

になるように規格化する. ${\displaystyle
\tilde{P}_{J-1}
=\sqrt{ \frac{1}{2(\mbox{J}-1)+1} } P_{J-1}
}$ より,
$\displaystyle w_j = \frac{1-\mu_j^2}
{(J \sqrt{ \frac{1}{2J-1} } P_{J-1}(\mu_j))^2 }
= \frac{(2J-1)(1-\mu_j^2)}
{(J P_{J-1}(\mu_j))^2 }$     (156)

となる.

まとめ

以上より

    $\displaystyle \int^{1}_{-1} f(\mu) d \mu = 2 \sum_{j=1}^{J} f(\mu_j) w_j$ (157)
    $\displaystyle w_j = \frac{(2J-1)(1-\mu_j^2)}{(J P_{J-1}(\mu_j))^2 }$ (158)

18 1818 18 1818 18 1818

E. 球面調和函数の離散的直交関係

ここでは 球面直交関数の離散的直交関係である選点直交性を示す.


$\displaystyle \sum_{j=1}^{J}
\sum_{i=1}^{I}
P_n^m (\mu_j) P_{n'}^{m'} (\mu_j)
\exp(i m \lambda_i) \exp(-i m' \lambda_i) w_j
= I \delta_{nn'} \delta_{mm'}$     (159)

ここで, $ i,j,m,m',n,n',I,J,M,N(m)$ は整数で, $1 \le i \le I , 1 \le j \le J, $ $0 \le \vert m\vert,\vert m'\vert \le M, \vert m\vert \le n \le N, \vert m'\vert \le n' \le N$ であり, ${\displaystyle M \le \left[ \frac{I}{2} \right] ,
N(m) \le J-1 }$ を満たす. また, $w_j$ は Gauss 荷重, ${\displaystyle \lambda_i=\frac{2\pi(i-1)}{I} }$, $\mu_j$$P_J (\mu)$ の零点である. $[ \ ]$ は それを越えない最大の整数を表す. これは, 有限な直交多項式系において成り立つ 選点直交性と呼ばれる性質である19.

この式を証明する. Legendre函数・陪函数の定義・(連続系での)直交性, Gauss の台形公式, Legendre函数の零点を用いた多項式の積分評価 を既知とすると,

    $\displaystyle \sum_{j=1}^{J}
\sum_{i=1}^{I}
P_n^m (\mu_j) P_{n'}^{m'} (\mu_j)
\exp(i m \lambda_i) \exp(-i m' \lambda_{i})
w_j$ (160)
  $\textstyle =$ $\displaystyle I \sum_{j=1}^{J}
P_n^m (\mu_j) P_{n'}^{m'} (\mu_j) w_j \delta_{mm'}$ (161)

ここで Gauss の台形公式を使った. 更に変形すると
    $\displaystyle \sum_{j=1}^{J}
\sum_{i=1}^{I}
P_n^m (\mu_j) P_{n'}^{m'} (\mu_j)
\exp(i m \lambda_i) \exp(-i m' \lambda_{i})
w_j$ (162)
  $\textstyle =$ $\displaystyle I \sum_{j=1}^{J}
P_n^m (\mu_j) P_{n'}^{m} (\mu_j) w_j$ (163)
  $\textstyle =$ $\displaystyle \frac{I}{2} \int_{-1}^1
P_n^m (\mu) P_{n'}^{m} (\mu) d\mu$ (164)

ここで, Gauss-Legendreの公式を使った. 更に, 連続系の Legendre函数・陪函数の直交性より
    $\displaystyle \sum_{j=1}^{J}
\sum_{i=1}^{I}
P_n^m (\mu_j) P_{n'}^{m'} (\mu_j)
\exp(i m \lambda_i) \exp(-i m' \lambda_{i})
w_j$ (165)
  $\textstyle =$ $\displaystyle I \delta_{nn'} \delta_{mm'}$ (166)

が得られる. 以上により, 離散化した球面調和関数の選点直交性が示された.





余談ではあるが, 直交多項式系においては 離散的な直交関係としては選点直交性のほかに 次のような直交関係も知られている20. $ \{ f_k(\mu) \}(k=0,1,2,\cdots) $$\ [ a ,b ] \ $ で定義された 重み $w(\mu)$, 規格化定数 $\lambda_k$ の直交多項式 ${\displaystyle
\left( \int_a^b f_k (\mu) f_{k'} (\mu) w(\mu) d \mu
= \lambda_k \delta_{kk'} \right) }$ とする. $\mu_j, \mu_{j'}(1 \le j,j' \le J)$$f_J(\mu)$ の零点, $w_j=w(\mu_j)$ とすれば, 選点直交性

$\displaystyle \sum_{j=0}^{J-1}
f_k (\mu_j) f_{k'} (\mu_{j}) w_j
= \lambda_k \delta_{kk'}$     (167)

のほかに,
$\displaystyle \sum_{k=0}^{J-1}
\frac{f_k (\mu_j) f_k (\mu_{j'}) }{\lambda_k}
= \frac{1}{w_j} \delta_{jj'}$     (168)

が成り立つ.

実際, Legendre函数 $\{ P_n \}(n=0,1,2,\cdots,J-1) $ については この関係が成り立つ. すなわち, $w_j$ を GCM で用いている Gauss 荷重として,

$\displaystyle \sum_{n=0}^{J-1} P_n (\mu_j) P_n (\mu_{j'})
= \frac{1}{w_j} \delta_{jj'}$     (169)

である. しかし, GCM では Legendre函数 $P_{J}$ の零点でのみ 値を計算することと, 波数切断の関係とから, Legendre陪函数 $\{ P_n^m \} (n=\vert m\vert,\vert m\vert+1,\vert m\vert+2,\cdots,N)$ の 離散的直交関係は意味がない21. Legendre函数の直交関係についても, 波数切断により $P_n$ $n=0,1,2,\cdots,N < J-1$ しか扱わないので22 実際には意味がない.

三角関数についても同様な離散的直交関係がある. 選点直交性

$\displaystyle \sum_{i=0}^{I-1}
\exp(i m \lambda_i) \exp(-i m' \lambda_i)
= I \delta_{mm'}$     (170)

のほかに,
$\displaystyle \sum_{m=-\frac{I}{2}+1}^{\frac{I}{2}}
\exp(i m \lambda_i) \exp(-i m \lambda_{i'})
= I \delta_{ii'}$     (171)

も成り立つ. (ただし, $I$ は偶数で $I=2M$. $I$ が奇数の場合には, $I=2M+1$ として, $m$についての和は ${\displaystyle -\frac{I-1}{2} \sim \frac{I-1}{2} }$ でとる. ) しかし GCM では, 波数切断により $\vert m\vert$ の最大値 $M$ ${\displaystyle \frac{I}{3} }$ 以下の値なので やはり意味がない23. 23 2323

F. スペクトルの係数と格子点値とのやり取り

ここではスペクトルの係数と格子点値との変換法について述べる. 実際の GCM 計算において必要になるのは

である.

F..1 スペクトルの係数と格子点値との値のやり取り

スカラー関数 $A(\lambda,\phi)$ の 格子点値とスペクトルの係数とのやり取りは 以下のとおりである. ただし, 格子点値は $A_{ij} \; (i=1,2,\cdots,I, \; j=1,2,\cdots,J)$ , スペクトルの係数は $\tilde{A}_n^m \;
(m=-M,-M+1, \cdots,M, \; n=\vert m\vert,\vert m\vert+1,\cdots,N(m))$ とする.


$\displaystyle A_{ij}$ $\textstyle \equiv$ $\displaystyle \sum_{m=-M}^{M} \sum_{n=\vert m\vert}^{N}
\tilde{A}_n^m
Y_n^m (\lambda_i,\phi_j)$ (172)
$\displaystyle \tilde{A}_n^m$ $\textstyle =$ $\displaystyle \frac{1}{I}
\sum_{i=1}^{I} \sum_{j=1}^{J}
A_{ij} Y_n^{m*} (\lambda_i,\phi_j) w_j$ (173)
$\displaystyle w_j$ $\textstyle =$ $\displaystyle \frac{(2J-1)(1-\sin^2 \phi_j)}
{(J P_{J-1}(\sin \phi_j))^2 }$ (174)

以後この文書では簡単のために, ${\displaystyle \sum_{m=-M}^{M} \sum_{n=\vert m\vert}^{N} }$ ${\displaystyle \sum_{m,n} }$ と, ${\displaystyle \sum_{i=1}^{I} \sum_{j=1}^{J} }$ ${\displaystyle \sum_{i,j} }$ と表記する.

F..2 スペクトルの係数と格子点値との値のやり取り〜東西微分編

まず,

\begin{eqnarray*}
g&\equiv& \DP{f}{\lambda}
\end{eqnarray*}

を考える.

東西微分($\lambda$ 微分)は次式で評価する.

$\displaystyle g_{ij}$ $\textstyle \equiv$ $\displaystyle \left[
\DP{}{\lambda} \left(
\sum_{m,n} \tilde{f}_n^m Y_n^m (\lambda, \phi)
\right)
\right]_{ij}$ (175)

すなわち,
$\displaystyle g_{ij}$ $\textstyle =$ $\displaystyle \sum_{m,n} im \tilde{f}_n^m
Y_n^m (\lambda_i,\phi_j)$ (176)

である. 変換公式 (A.103)で $A$$g$ とみなしたものと (A.106) とを比較すれば明らかに24,


$\displaystyle \tilde{g}_n^m$ $\textstyle =$ $\displaystyle im \tilde{f}_n^m$ (177)

よって,
$\displaystyle \tilde{g}_n^m$ $\textstyle =$ $\displaystyle \frac{1}{I}
\sum_{i,j} im f_{ij} Y_n^{m*} (\lambda_i, \phi_j) w_j$ (178)

である.

次に,

\begin{eqnarray*}
h \equiv \frac{g}{r \cos^2 \phi}
= \frac{1}{r \cos^2 \phi} \...
...left[ =
\DP{}{x} \left( \frac{f}{\cos \phi} \right)
\right]
\end{eqnarray*}

とする. $f$$h$ とのやり取りを考える. (A.105) より明らかに,

\begin{eqnarray*}
h_{ij} &=& \frac{1}{r \cos^2 \phi_i} g_{ij} \\
∴\ h_{ij} &...
...2 \phi_j}
\sum_{m,n} im \tilde{f}_n^m Y_n^m (\lambda_i, \phi_j)
\end{eqnarray*}

一方, (A.107) より

$\displaystyle \tilde{h}_n^m$ $\textstyle =$ $\displaystyle \widetilde{
\left[
\DP{}{\lambda} \left( \frac{f}{r \cos^2 \phi} ...
...ht)
\right]_n^m
}
= im \widetilde{ \left( \frac{f}{r \cos^2 \phi} \right)_n^m }$  
  $\textstyle =$ $\displaystyle \frac{1}{I}
\sum_{i,j}
im
\left( \frac{f}{r \cos^2 \phi} \right)_{ij}
Y_n^{m*} (\lambda_i, \phi_j) w_j$  
  $\textstyle =$ $\displaystyle \frac{1}{I}
\sum_{i,j} im f_{ij}
Y_n^{m*} (\lambda_i \phi_j)
\frac{w_j}{r \cos^2 \phi_j}$ (179)

F..3 スペクトルの係数と格子点値との値のやり取り〜南北微分編

まず,

\begin{eqnarray*}
p \equiv \DP{f}{\mu}
\end{eqnarray*}

を考える.

南北微分($\mu$微分)は次式で評価する.

$\displaystyle p_{ij} \equiv \left[ \DP{}{\mu}
\left( \sum_{m,n} \tilde{f}_n^m Y_n^m \right)
\right]_{ij}$     (180)

すなわち,
$\displaystyle p_{ij} = \sum_{m,n} \tilde{f}_n^m
\left. \DD{P_n^m}{\mu} \right\vert _j
\exp(im \lambda_i)$     (181)

である.

\begin{eqnarray*}
p_n^m &=& \frac{1}{I} \sum_{i,j} p_{ij} Y_n^{m*} w_j \\
&=&...
...eft. \DD{P_n^m}{\mu}\right\vert _j
\exp(-im \lambda_i) w_j \\
\end{eqnarray*}

ここで, 2行目から 3行目の等号では,
    $\displaystyle \sum_{i=1}^I \sum_{j=1}^J
f_{n'}^{m'}
P_{n}^{m}(\mu_j) \exp(im \lambda_i)
\left. \DD{P_{n'}^{m'}}{\mu}\right\vert _j
\exp(-im' \lambda_i)
w_j$  
  $\textstyle =$ $\displaystyle - \sum_{i=1}^I \sum_{j=1}^J
f_{n'}^{m'}
\left. \DD{P_{n}^{m}}{\mu}\right\vert _j
\exp(-im \lambda_i)
P_{n'}^{m'}(\mu_j) \exp(im' \lambda_i)
w_j$ (182)

を用いた25.

次に,

\begin{eqnarray*}
q \equiv (1-\mu^2) \DP{f}{\mu} = (1-\mu^2) p
\end{eqnarray*}

とする.

(A.110) より明らかに,

\begin{eqnarray*}
q_{ij} = (1-\mu_j^2) \sum_{m,n} \tilde{f}_n^m
\left. \DD{P_n^m}{\mu} \right\vert _j
\exp(im \lambda_i)
\end{eqnarray*}

である. 一方,

\begin{eqnarray*}
\tilde{q}_n^m &=& \frac{1}{I} \sum_{i,j} q_{ij} Y_n^{m*} w_j ...
... (1-\mu^2)P_n^m \right) \right\vert _j
\exp(-im \lambda_i) w_j
\end{eqnarray*}

が成り立つ. ここで, 2行目から3行目において,

\begin{eqnarray*}
& & \sum_{i=1}^I \sum_{j=1}^J
f_{n'}^{m'} (1-\mu_j^2)
P_{n}...
...xp(-im \lambda_i)
P_{n'}^{m'}(\mu_j) \exp(im' \lambda_i)
w_j
\end{eqnarray*}

を用いた26.

F..4 速度の格子点値から発散・渦度のスペクトルの係数への変換

速度場を

$\displaystyle (u,v) = \left( \frac{U}{\cos \phi}, \frac{V}{\cos \phi} \right)$     (183)

とする. ここでは, $(U_{ij},V_{ij})$ から $\tilde{D}_n^m, \; \tilde{\zeta}_n^m$ を求める27.



まず,

$\displaystyle D = \frac{1}{r(1-\mu^2)} \DP{U}{\lambda}
+ \frac{1}{r} \DP{V}{\mu}$     (184)

より,


$\displaystyle \tilde{D}_n^m$ $\textstyle =$ $\displaystyle \frac{1}{I}
\sum_{ij}
\left(
\frac{im}{r \cos^2 \phi_j} U_{ij} P_...
...} V_{ij} \left. \DD{P_n^m}{\mu} \right\vert _j
\right)
\exp(- im \lambda_i) w_j$  
  $\textstyle =$ $\displaystyle \frac{1}{I}
\sum_{ij}
\left(
im U_{ij} P_n^m
- V_{ij}
\cos^2 \phi...
...m}{\mu} \right\vert _j
\right)
\exp(- im \lambda_i)
\frac{w_j}{r \cos^2 \phi_j}$ (185)

同様に,

$\displaystyle \zeta = \frac{1}{r(1-\mu^2)} \DP{V}{\lambda}
- \frac{1}{r} \DP{U}{\mu}$     (186)

より,


$\displaystyle \tilde{\zeta}_n^m$ $\textstyle =$ $\displaystyle \frac{1}{I}
\sum_{ij}
\left(
\frac{im}{\cos^2 \phi} V_{ij} P_n^m
...
...} U_{ij} \left. \DD{P_n^m}{\mu} \right\vert _j
\right)
\exp(- im \lambda_i) w_j$  
  $\textstyle =$ $\displaystyle \frac{1}{I}
\sum_{ij}
\left(
im V_{ij} P_n^m
+ U_{ij}
\cos^2 \phi...
...m}{\mu} \right\vert _j
\right)
\exp(- im \lambda_i)
\frac{w_j}{r \cos^2 \phi_j}$ (187)

F..5 $\chi ,\psi $ のスペクトルの係数から速度の格子点値への変換

ここでは $\chi_n^m,\psi_n^m$ から $U_{ij},V_{ij}$を求める方法を記す.

まず,

$\displaystyle U = - \frac{1-\mu^2}{r} \DP{\psi}{\mu}
+ \frac{1}{r} \DP{\chi}{\lambda}$     (188)

より,
$\displaystyle U_{ij}$ $\textstyle =$ $\displaystyle \sum_{m,n}
\left(
- \frac{\cos^2 \phi_j}{r} \tilde{\psi}_n^m
\lef...
...+ \frac{1}{r} im \tilde{\chi}_n^m P_n^m(\sin \phi_j)
\right)
\exp(im \lambda_i)$  
  $\textstyle =$ $\displaystyle \sum_{m,n}
\left(
- \frac{1}{r} \tilde{\psi}_n^m
\cos^2 \phi_j \l...
...+ \frac{1}{r} im \tilde{\chi}_n^m P_n^m(\sin \phi_j)
\right)
\exp(im \lambda_i)$ (189)

である. 同様に,
$\displaystyle V = \frac{1}{r} \DP{\psi}{\lambda}
+ \frac{1-\mu^2}{r} \DP{\chi}{\mu}$     (190)

より,
$\displaystyle V_{ij}$ $\textstyle =$ $\displaystyle \sum_{m,n}
\left(
\frac{1}{r} im \tilde{\psi}_n^m P_n^m (\sin \ph...
...de{\chi}_n^m
\left. \DD{P_n^m}{\mu} \right\vert _{j}
\right)
\exp(im \lambda_i)$  
  $\textstyle =$ $\displaystyle \sum_{m,n}
\left(
\frac{1}{r} im \tilde{\psi}_n^m P_n^m (\sin \ph...
...cos^2 \phi_j \left. \DD{P_n^m}{\mu} \right\vert _{j}
\right)
\exp(im \lambda_i)$ (191)

である. 27 2727

G. スペクトルの係数同士の関係

ここではスペクトルの係数同士の便利な公式を挙げておく. $g = \DP{f}{\lambda}$ の時

$\displaystyle \tilde{g}_n^m = im \tilde{f}_n^m$     (192)

$h = \nabla_H^2 f$ の時
$\displaystyle \tilde{h}_n^m = -\frac{n(n+1)}{r^2} \tilde{f}_n^m$     (193)




(A.122) については 「スペクトルの係数と格子点値とのやり取り」に証明を示した. ここでは, (A.123) について証明しておく.

微分評価の定義より,

\begin{eqnarray*}
h_{ij}
= \left. \left(\nabla_H^2 \sum_{m,n}
\tilde{f}_n^m...
...{n(n+1)}{r^2}
\left. \tilde{f}_n^m Y_n^m \right\vert _{ij} \\
\end{eqnarray*}

である. ところで,

\begin{eqnarray*}
h_{ij} = \sum_{m,n}
\left. \tilde{h}_n^m Y_n^m \right\vert _{ij} \\
\end{eqnarray*}

である. この2つの式の右辺に左から ${\displaystyle \sum_{i,j} Y_{n'}^{m'*}\vert _{ij} }$ を演算して 比較すると,

\begin{eqnarray*}
\tilde{f}_{n'}^{m'}
= -\frac{n(n+1)}{r^2} \tilde{h}_{n'}^{m'}
\end{eqnarray*}

を得る.

27 2727

H. 波数切断

GCM では, 物理量を 球面調和函数 $P_n^m(\sin \phi) \exp(im \lambda)$ で展開したり 波数空間で計算するときに, 計算資源の都合上, ある一定波数以下の波数のみを考慮して計算する. そのことを波数切断するという28. 以下ではまず, 切断の基礎知識として切断の仕方・流儀を述べ, ついで, 切断における事情を述べた上で切断波数の決め方を記す.

H..1 波数切断の仕方

波数切断の仕方については, 東西波数($m$), 南北波数($n-m$)の それぞれの切断の方法にいくつかの流儀がある. 一般によく用いられるものは 三角形切断(Triangle), 平行四辺形切断(Rhomboidal : 偏菱形)と 呼ばれるものである. 三角形切断の場合について計算する波数領域を波数平面上に書くと (A.1)のようになる. 平方四辺形切断の場合は, (A.2)である.

図 A.1: 三角形切断の場合の波数領域
\begin{figure}\begin{center}
\begin{picture}(300,100)(10,-10)
\put(20 , 10){\v...
...00,-20) \{ shortstack\{三角形切断\}\}
\par
\end{picture}\end{center}\end{figure}

図 A.2: 平方四辺形切断の場合の波数領域
\begin{figure}\begin{center}
\begin{picture}(300,170)(10,-20)
\put(20 , 10){\v...
...20) \{ shortstack\{平方四辺形切断\}\}
\par
\end{picture}\end{center}\end{figure}

三角形切断, 平行四辺形切断, という名称は 波数平面上($(n,m)$平面)での形状による29.

より一般的な切断方法は五角形切断 ((A.3)) である.

図 A.3: 五角形切断の場合の波数領域
\begin{figure}\begin{center}
\begin{picture}(300,170)(10,-20)
\put(20 , 10){\v...
...00,-20) \{ shortstack\{五角形切断\}\}
\par
\end{picture}\end{center}\end{figure}

三角形切断, 平行四辺形切断はそれぞれ, 五角形切断において

であるような特別な場合である30.

三角形切断と平行四辺形切断の違いについて, 世の中では次のように言われている31.

H..2 切断波数の決め方

ここでは切断波数と南北格子点数の決め方について記す. これらは 切断の仕方を決めた後に, 使用する計算資源がネックになって決まる. その際, FFT の仕様, aliasing の回避, という2つの数値的な事情を 考慮した上で決める必要がある.

FFT の仕様の事情というのは, 話は簡単で, 東西方向に 「格子 $\Leftrightarrow$ スペクトル」 変換する ために用いる FFT が 効率よく動くための格子点数・波数がある33ことである.

一方, aliasing に関する事情は複雑である. ここで扱っているスペクトルモデルでは, 格子点でのみ値を計算している. いわゆるスペクトルを使うのは, 単に格子点上での水平微分項の評価をする時のみである. その意味で, 「微分の評価にのみスペクトルを用いるグリッドモデル」 と言ってもよい. そのように受け止めると, 格子点値を"正しく"計算することを目指し, また, 考慮する波数は 厳密にスペクトルの係数と格子との変換を行なうことのできる波数, すなわち変換において情報の落ちないだけの波数を とらねばならないように思える. ところが実際には, スペクトルモデル的な配慮 -- ある波数以下についてのみ 正しく計算し, それ以上の波数については計算しない -- により 切断波数・格子点数が決められている. また, 後述する理由により 情報は(非線形 aliasing のことを考えずとも) 必ず落ちてしまうのである34.

さて, 以下では aliasing に関する事情を具体的に述べながら, 切断波数に対する格子点数の決め方を記そう. 球面上に連続分布している物理量を球面調和函数で展開する. ある波数 $M,N(m)$ 以下(例えば, T42 ならば $M=42, N=42$ ) については 線形項・非線形項の両方について 厳密に計算できるように $I,J$ を決めることを目指す.

$M,N$ を仮に固定したとして, まずは線形項について 切断波数以下のスペクトルの係数のわかっている物理量 $A$ を 格子点値に変換しさらにスペクトルの係数に正しくもどすことを考える. $A$ $ -M \le m \le M,\ \vert m\vert \le n \le N(m)$$m,n$ については $\tilde{A}_n^m$ が わかっているとする. 格子点値は, $1 \le i \le I,1 \le j \le J$ について

$\displaystyle A_{ij}$ $\textstyle \equiv$ $\displaystyle \sum_{m=-M}^M \sum_{n=\vert m\vert}^N
\tilde{A}_n^m P_n^{m}(\sin \phi_j)
\exp(im \lambda_i)$ (194)

で与えられる. これらの格子点値から逆に $\tilde{A}_n^m ( -M \le m \le M, \ \vert m\vert \le n \le N)$ を計算する. 離散化した系での積分を Gauss の公式, Gauss-Legendre の公式で評価すれば,
$\displaystyle \tilde{A}_n^m$ $\textstyle =$ $\displaystyle \frac{1}{I}
\sum_{i=1}^I \sum_{j=1}^J
A_{ij} P_n^{m}(\sin \phi_j)
\exp(-im \lambda_i) w_j$ (195)

である. ここで, $w_j$$\phi_j$ における重みである. $A_{ij}$ の定義を代入すれば,
$\displaystyle \tilde{A}_n^m$ $\textstyle =$ $\displaystyle \frac{1}{I} \sum_{i=1}^I \sum_{j=1}^J
\left(
\sum_{m'=-M}^{M} \su...
... \phi_j)
\exp(im' \lambda)
\right)
P_n^{m}(\sin \phi_j)
\exp(-im \lambda_i) w_j$  
  $\textstyle =$ $\displaystyle \frac{1}{I} \sum_{m'=-M}^{M} \sum_{n'=\vert m'\vert}^N
\tilde{A}_...
...i(m'-m) \lambda)
\sum_{j=1}^J
P_n^{m}(\sin \phi_j)
P_{n'}^{m'}(\sin \phi_j) w_j$ (196)

となる. この計算が $\tilde{A}_n^m$ を 正しく評価している(すなわち元にもどる)ための $I,J$ の条件は, $ -M \le m \le M,\ \vert m\vert \le n \le N$ を満たす $m,n$ について
    $\displaystyle \sum_{i=1}^I \exp(i(m'-m) \lambda) = I \delta_{mm'}$ (197)
    $\displaystyle \sum_{j=1}^J P_n^{m}(\sin \phi_j)
P_{n'}^{m}(\sin \phi_j) w_j
= \delta_{nn'}$ (198)

が成り立つことである. 三角関数の和による評価が正しいための条件は, ここに登場する波数 $\vert m'-m\vert$ が最大で $2M$ の値をとるので, Gauss の公式の適用条件より, 格子点数 $I$$I \ge 2M+1$ を満たすことである. Legendre函数の積の和による評価が正しいための条件は, ここに登場する計算が $n+n'$ 次の多項式35の 評価であることから, Gauss - Legendre の公式の適用条件より, 格子点数 $J$ $2J-1 \ge max[n+n'] = 2 max[N]$ を 満たすことである. ここで, $max[n+n']$$n+n'$ の最大値を, $max[N]$$N$ の最大値を表す.

ちなみに, 格子点値からスペクトルの係数に変換し格子点値にもどす という立場からすれば, この Gauss-Legendre の公式の適用条件というのが 情報を落とさずには済まない理由である36. このことを以下に述べる. 情報を落とさずに 格子点値をスペクトルの係数に変換し格子点値にもどすには, あらゆる東西波数について 南北方向の格子点数 $J$ と同じだけの個数のLegendre函数が必要である. 東西波数 $m$ の場合, 登場するLegendre陪函数の $n$ $n=\vert m\vert, \vert m\vert+1, \cdots, \vert m\vert+J-1$ である. $ P_n^m P_{n'}^{m} $ の次数は $n+n'$ であるから, 最大で $2J+2\vert m\vert-2$ である. これが $2J-1$ 以下になるのは $m=0$ の時のみである. $m \neq 0$ の場合は高次のLegendre函数は計算してはならない. つまり情報を落とさざるをえない37.

改めて $M,N$ を固定するという立場にもどって, 切断波数以下のスペクトルの係数のわかっている物理量 $B,C$ の積から それらの格子点値を用いて $B$$C$ との積(非線形項) $A$ のスペクトルの係数を正しく求めるための $I,J$ の条件を考える.

    $\displaystyle A= BC$ (199)
    $\displaystyle B= \sum_{m=-M}^{M} \sum_{n=\vert m\vert}^{N}
\left( \tilde{B}_n^m \exp(im \lambda) \right)
P_n^m(\sin \phi)$ (200)
    $\displaystyle C= \sum_{m=-M}^{M} \sum_{n=\vert m\vert}^{N}
\left( \tilde{C}_n^m \exp(im \lambda) \right)
P_n^m(\sin \phi)$ (201)

なる物理量 $A,B,C$ があるとする38. $B$,$C$ $ -M \le m \le M,\ \vert m\vert \le n \le N$ における スペクトルの係数 $\tilde{B}_n^m,\tilde{C}_n^m$ を用いて $A$ の スペクトルの係数 $\tilde{A}_n^m$ $0 \le m \le M,\ \vert m\vert \le n \le N$ については 正しく計算することを考える.


$\displaystyle \tilde{A}_n^m$ $\textstyle \equiv$ $\displaystyle \widetilde{(BC)_n^m}$ (202)
  $\textstyle =$ $\displaystyle \frac{1}{I} \sum_{i=1}^I \sum_{j=1}^J
B_{ij} C_{ij}
P_n^{m}(\sin \phi_j) \exp(-im \lambda_i) w_j$ (203)
  $\textstyle =$ $\displaystyle \frac{1}{I} \sum_{i=1}^I \sum_{j=1}^J
\left( \sum_{m'=-M}^{M} \su...
...ert}^N
\tilde{B}_{n'}^{m'}
\exp(im' \lambda_i)
P_{n'}^{m'}(\sin \phi_j)
\right)$ (204)
    $\displaystyle \vspace{4cm}
\times
\left( \sum_{m''=-M}^{M} \sum_{n''=\vert m''\...
...P_{n''}^{m''}(\sin \phi_j)
\right)
P_n^{m}(\sin \phi_j) \exp(-im \lambda_i) w_j$ (205)
  $\textstyle =$ $\displaystyle \frac{1}{I}
\sum_{m'=-M}^{M} \sum_{n'=\vert m'\vert}^N
\sum_{m''=-M}^{M} \sum_{n''=\vert m''\vert}^N
\tilde{B}_{n'}^{m'} \tilde{C}_{n''}^{m''}$ (206)
    $\displaystyle \vspace{4cm}
\times
\sum_{i=1}^I
\exp(i(m'+m''-m) \lambda_i)
\sum...
...^J
P_{n'}^{m'}(\sin \phi_j)
P_{n''}^{m''}(\sin \phi_j)
P_n^{m}(\sin \phi_j) w_j$ (207)

この計算が $\tilde{A}_n^m$ $0 \le m \le M,\ \vert m\vert \le n \le N$ について 正しく評価しているための, $I,J$ の条件を 線形項の場合と同様に考えると, 格子点数 $I$$I \ge 3M+1$ を, 格子点数 $J$ $2J-1 \ge max[n+n'+n''] = 3 max[N]$ を 満たすことである. ここで, $max[n+n'+n'']$$n+n'+n''$ の最大値を, $max[N]$$N$ の最大値を表す.

再び 格子点値からスペクトルの係数に変換し格子点値にもどす という立場からすれば, これらの $I,J$ に関する条件から, 南北成分のみならず, 東西成分についても 変換によって情報が落ちてしまうことがわかる.

これまでに述べた $M,N$ を固定したときに 格子点数 $I,J$ がとらねばならない個数について, 線形項・非線形項の2つの場合のうち条件が厳しいのは, 明らかに非線形項の場合である. この条件以下の格子点数しかとらない場合には, aliasing をおこすことになる.

以上, FFT, aliasing という2つの事情を考えて 格子点数と切断波数とは同時に決められる. 具体的手順は以下のとおりである.

  1. 波数切断の仕方を決める.
  2. FFT のかけやすい数を選ぶ. それを東西格子点数 $I$ とする.
  3. 東西方向の波数の最大値 $M$ ${\displaystyle M= \left[ \frac{I-1}{3} \right] }$ にする. ただし $[ \ ]$ はそれを越えない最大の整数を表す記号である.
  4. 最大全波数 $N_{max}$ を決める. 三角形切断ならば $N_{max}=M$, 平行四辺形切断ならば $N_{max}=2M$ である.
  5. 南北方向の格子点数 $J$ $J \ge \frac{3N_{max}+1}{2}$ を 満たす数に選ぶ. (今の GCM では偶数でなくてはならない. )
例えば, T42 の場合には $M=42, N=42$, 東西格子点数 $I$ が 128, 南北格子点数 $J$ が 64 である. R21 の場合には $M=21, N=42$, 東西格子点数 $I$ が 64, 南北格子点数 $J$ が 64 である.



参考までに, 線形モデルの場合について決め方を示しておく.

  1. 波数切断の仕方を決める.
  2. FFT のかけやすい数を選ぶ. それを東西格子点数 $I$ とする.
  3. 東西方向の波数の最大値 $M$ ${\displaystyle M= \left[ \frac{I}{2} \right] }$ にする. ただし $[ \ ]$ は それを越えない最大の整数を表す記号である39.
  4. 最大全波数 $N_{max}$ を決める. 三角形切断ならば $N_{max}=M$, 平行四辺形切断ならば $N_{max}=2M$ である.
  5. 南北方向の格子点数 $J$ $J \ge \frac{2N_{max}+1}{2}$ を 満たす数に選ぶ.
例えば, 三角形切断の場合には, $I=128$ とすると, $M=64$, $N=64$, $J=65$ となる. つまり T64 では $I=128, J=65$ である. 平方四辺形切断の場合には, $I=64$ とすると, $M=32$, $N=64$, $J \ge 65$ となる. つまり R32 では $I=64, J=65$ でよい40.

40 4040

I. スペクトルモデルと差分モデル

世の中の多くの GCM の離散化の方法としては, 鉛直方向については必ず レベルと称する差分による離散化を行なうが, 水平方向については, 差分する方法(この方法を用いるモデルをグリッドモデルという)と 球面調和函数で展開してその係数の時間変化を計算する方法 (力学過程において41この方法を用いるモデルを スペクトルモデルという)と が用いられる. その二つの方法については一長一短がある. ここでは双方の特徴について列挙しておく42.

ちなみに, dcpam はスペクトルモデルに分類される.



... 岩波公式集10
... 三角関数11
$\exp(i m \lambda)$ $\int_0^{2 \pi} \exp(im\lambda) \exp(-im'\lambda)
d \lambda = 2 \pi \delta_{mm'}$ を満たす. ただし, $m,m'$ は整数である.
... のである12
この特徴を 言い替えれば, 全波数 $n$ の球面調和函数の重ね合わせで表現できる分布関数は 座標系を回転させた系においても 全波数 $n$ の球面調和函数の重ね合わせで表現できることになる.
... 固有値であることによっている13
$\nabla_H^2 \equiv
\frac{1}{r^2}
\left[ \DP{}{\phi} \left( \cos \phi \DP{}{\phi} \right)
+ \frac{1}{\cos^2 \phi} \DP[2]{}{\lambda}
\right]
$ の, 固有値を $- \frac{n(n+1)}{r^2}$ とする固有関数であることと, スカラー演算子 $\nabla_H^2$ が 座標系の回転に関して不変な演算子であることとに起因する.

すなわち, $\nabla_H^2 Y^m_n(\lambda, \phi)
=- \frac{n(n+1)}{r^2} Y^m_n(\lambda,\phi)$ より, 球面調和函数 $Y_n^m \exp(im\lambda)$ は 固有値を $- \frac{n(n+1)}{r^2}$ とする $\nabla_H^2$ の固有関数である. $\{Y_n^m \vert n=0,1,2,\cdots,\ m= -n, -n+1,\cdots,n \}$ の完全直交性より, $\{Y_n^m \vert m= -n, -n+1,\cdots,n \}$ $\nabla^2_H f = -\frac{n(n+1)}{r^2} f$ の 解空間を張っている基底である.

座標系を回転させて, 新たな座標系での球面調和函数 $Y_n^m(\lambda',\phi')$ の和の形で 前の座標系での球面調和函数 $Y_n^m(\lambda,\phi)$ を 表現することを考えよう.

絶対系で見て同じ位置の値を比べると, 2次元ラプラシアンを演算した値は不変なので, 前の座標系での球面調和函数 $Y_n^m(\lambda',\phi')$ は 新たな座標系においても $\nabla^{'2}_H Y_n^m = -\frac{n(n+1)}{r^2} Y_n^m $ の解である. 新たな座標系の球面調和函数の集合 $\{ Y_n^m(\lambda', \phi') \vert m=-n,-n+1, \cdots, n \} $ $\nabla^{'2}_H Y_n^m = -\frac{n(n+1)}{r^2} Y_n^m $ の 解空間の基底である. したがって, 前の座標系の球面調和函数は 新たな座標系の球面調和函数の和の形で書ける.

... 14
(2005/4/4 石渡) 関数形も書いておきたい. グラフは自分で描きたい.
...P15
(2004/5/27 石渡) ここにあるのは変かも.
... の整数倍にならない)16
等比級数の和を直接計算しても良い.
$\displaystyle \sum^{I}_{n=1} \exp \left\{ im \frac{2 \pi (n-1)}{I} \right\}
= \...
...im 2 \pi}{I} } }
= \frac{1 - e^{im 2 \pi} }
{1 - e^{ \frac{im 2 \pi}{I} } }
= 0$     (127)

... 直すことにする17
私が混乱しないためにこのような手続きを踏む. 実際, 公式集を含む他の文献には $\tilde{P}_n^m$ の公式が書かれていることが多いので, このように書く方が他と参照しやすいであろう.
... になる18
このことは $L-f$$J-1$ 次以下の多項式であること, $J$ 個の零点 $\mu_j$ を持つことから明らか.
... 選点直交性と呼ばれる性質である19
別の離散的直交関係については後で述べる.
... 次のような直交関係も知られている20
以下については, 森, 1984 「数値解析法」が詳しい.
... 離散的直交関係は意味がない21
そもそも, ここで述べている直交関係は $f_k(k=0,1,2,\cdots,K-1)$$k$次多項式であるような直交多項式系において 成り立つものである. Legendre陪函数は $m$ が奇数のときは多項式でないし, $m$ が偶数であっても $P_n^m$$n$次多項式であって, $n-m$ 次多項式ではない. その場合にも直交多項式の議論を拡張して ここで述べている直交関係を 使えるのか, については未調査である.
... しか扱わないので22
T42 ならば, $m=0$$J=63, N=42$, R21 ならば, $m=0$$J=63, N=21$, である.
... やはり意味がない23
T42 ならば $I=128$ に対して $M=42$ , R21 ならば $I=64$ に対して $M=21$ である.
... とを比較すれば明らかに24
より正確には, ${\displaystyle (g_{ij}=)
\sum_{m,n} im \tilde{f}_n^m Y_n^m
= \sum_{m,n} \tilde{g}_n^m Y_n^m }$ の両辺に左から ${\displaystyle \sum_{i,j} Y_n^{m*}(\lambda_i,\phi_j)w_j}$ を演算すれば, $im' \tilde{f}_{n'}^{m'}
= \tilde{g}_{n'}^{m'} $ として 得られる.
... を用いた25
... を用いた26
... を求める27
この項の計算については 後ろの補遺 「$\mu$ 微分についての補助計算」 を参照せよ.
... そのことを波数切断するという28
...平面)での形状による29
平方四辺形切断には, $n$ の最大値を $m$ の最大値の2倍に しないようなとり方もある. 詳しくは五角形切断に関する脚注参照.
... であるような特別な場合である30
... 世の中では次のように言われている31
... 経度方向のみならず緯度方向にも一定である32
分解能が緯度方向に変化することについては, 平行四辺形切断に限らず, 三角形切断以外のどれでも起こる.
... 効率よく動くための格子点数・波数がある33
コード依存性がある. 通常, 2のべき乗が好ましいとされる. コードによっては, 2,3,5 のべき乗の積でもよいものもある.
... 必ず落ちてしまうのである34
実際の GCM では 格子点値からスペクトルに変換する際に情報は落ちている. したがって, 格子 - スペクトル - 格子という変換を 行なうと元にはもどらない.
例えば T42 の場合, 自由度は $ 1+ (2 \times1 +1) + \cdots
+ (2 \times 42 +1) = 43^2 = 1849$ に対して 格子点数は $ 128 \times 64 = 8192 $ である. R21 の場合も, 自由度は $ (2 \times 21 +1) \times (21+1)= 946$ に対して, 格子点数は $ 64 \times 64 = 4096 $ である. すなわち, $3/4$ 以上の情報は 格子点値からスペクトルに変換するときに落ちている.
工夫すれば 情報が落ちないうまい方法があるかも知れないが, 今のところ見つけていないし多分見つからない.
もちろん, スペクトル - 格子 - スペクトルという変換では 元にもどる(ように決めている).
... 次の多項式35
ここで, 三角関数の和が $I \delta_{mm'}$ となることを用いた. 一般には($m,m'$ の偶奇が一致しない場合には) $ P_n^m P_{n'}^{m'} $ は多項式にならない.
... 情報を落とさずには済まない理由である36
Gauss の公式の適用条件と情報欠落との関係について コメントしておく. 格子点数 $I$ が奇数の場合には, スペクトルで同じ情報量を持つためには 波数 $\frac{I-1}{2}$ までを考慮すればよいので, 情報は欠落しないことは明らかである. 一方, $I$ が偶数の場合には, 情報は欠落させないためには 波数 $\frac{I}{2}$ が必要であるが, この波数は Gauss の公式の適用条件を満たさない. しかしこの場合にも, (私は根拠を調べていないが, 少なくとも) 経験的には FFT および 逆FFT によって 格子 - スペクトル - 格子変換によって 情報が落ちないことが知られている.
... つまり情報を落とさざるをえない37
この事情により, 非線形項の場合を考えて さらに著しく落とすことが必要になることが次節からわかる.
... があるとする38
$A,B,C$ とも 実数である. すなわち, $\tilde{B}_n^m = \tilde{B}_n^{m*}, {\it etc.}$ となっている.
... それを越えない最大の整数を表す記号である39
ここで, $I$ が偶数のときについては Gauss の公式の適用条件を越えて 最大波数 $\frac{I}{2}$ まで 計算できるという知識を用いた.
... でよい40
これらの場合でも, 南北方向の細かい情報は 格子 - スペクトル - 格子変換によって 落ちていることに注意せよ.
... (力学過程において41
adjustment 等の意味をなど考えると, 特に物理過程においては, 格子点で考える方が物理的に当然であるように思う. そのためであろうか, スペクトルモデルである東大版GCM でも 物理過程を格子点で計算している. 他のスペクトルモデルについても そうであるかどうかは未調査.
... ここでは双方の特徴について列挙しておく42
出典は, スペクトル法による数値予報(その原理と実際) (1.6)
... 単純には扱えない43
問題点その1. グリッドモデルでは 緯度経度図で等間隔に格子点をとると, 極でも CFL を満たすようにするために, 時間差分を細かくしなければならない. 他は未調査.

next up previous
: 6 謝辞 : DCPAM3 第2部 離散化 : 4 参考文献
Morikawa Yasuhiro 平成17年12月6日