3 変分法

変分法は Sasaki (1958) により、物理法則による拘束を課す解析手法として提案され、Talagrand and Courtier (1987) に制御理論に基づく理論的な裏付けがされている。 Lorenc (1986) は、数値天気予報における様々な解析手法を比較し、変分法がベイズ推定から導かれることを指摘している。ここでは、Tsuyuki and Miyoshi (2007) に基づいて、変分法データ同化について述べる。

3.1 3次元変分法

3次元変分法同化（3DVar: three-dimensional variational assimilation）は、第一推定値\(\mathbf{x}^\mathrm{f}\)及び観測\(\mathbf{y}^\mathrm{o}\)が与えらたときに、後験pdfを最大化するものを解析値\(\mathbf{x}^\mathrm{a}\)とする。

\[ \begin{aligned} \mathbf{x}^\mathrm{a} &= \underset{\mathbf{x}}{\text{arg max}}[p(\mathbf{x}|\mathbf{x}^\mathrm{f},\,\mathbf{y}^\mathrm{o})]\\ &= \underset{\mathbf{x}}{\text{arg max}}[p(\mathbf{x}^\mathrm{f}|\mathbf{x})p(\mathbf{y}^\mathrm{o}|\mathbf{x})p(\mathbf{x})] \end{aligned} \tag{3.1}\] ガウス型のpdfを仮定し、負の対数pdfでコスト函数\(J(\mathbf{x})\)を定義すると、 \[ J(\mathbf{x}) = \frac{1}{2}(\mathbf{x}-\mathbf{x}^\mathrm{f})^\mathrm{T}\mathbf{B}^{-1}(\mathbf{x}-\mathbf{x}^\mathrm{f}) + \frac{1}{2}(h(\mathbf{x})-\mathbf{y}^\mathrm{o})^\mathrm{T}\mathbf{R}^{-1}(h(\mathbf{x})-\mathbf{y}^\mathrm{o}) + J_\mathrm{c}(\mathbf{x}) \tag{3.2}\] となる。ここで、\(\mathbf{B},\,\mathbf{R}\)はそれぞれ背景及び観測誤差共分散行列、\(h(\mathbf{x})\)は状態変数を観測変数に変換する観測演算子である。式(Equation 3.2)の右辺第一、二項はそれぞれ第一推定値、観測値からのずれを表す。式(Equation 3.2)の右辺第三項は、式(Equation 3.1)の右辺第三項に対応し、制約項と呼ばれている。

3DVarでは、数値最適化により\(J(\mathbf{x})\)の最小値として解析値\(\mathbf{x}^\mathrm{a}\)を求める。 \(h(\mathbf{x})\)の接線型演算子を\(\mathbf{H}\)で表すと、制約項を無視した解析値は \[ \mathbf{x}^\mathrm{a} = \mathbf{x}^\mathrm{f} + (\mathbf{B}^{-1} + \mathbf{H}^\mathrm{T}\mathbf{R}^{-1}\mathbf{H})^\mathrm{-1}\mathbf{H}^\mathrm{T}\mathbf{R}^{-1}[\mathbf{y}^\mathrm{o} - h(\mathbf{x}^\mathrm{f})] \tag{3.3}\] となる。

3.2 4次元変分法

4次元変分法同化（4DVar: four-dimensional variational assimilation）では、同化窓と呼ばれる期間に対するMAP推定を行う。時刻レベル\(i-1\)から\(i\)までの数値モデルによる予報を \[ \mathbf{x}_i = M_{i-1,i}(\mathbf{x}_{i-1})\;i=1,\dots,n \tag{3.4}\] で表す。 4DVarのコスト函数は \[ J(\mathbf{x_0}) = \frac{1}{2}(\mathbf{x}_0-\mathbf{x}_0^\mathrm{f})^\mathrm{T}\mathbf{B}^{-1}(\mathbf{x}_0-\mathbf{x}_0^\mathrm{f}) + \frac{1}{2}\sum_{i=0}^n(h_i(\mathbf{x}_i)-\mathbf{y}_i^\mathrm{o})^\mathrm{T}\mathbf{R}_i^{-1}(h_i(\mathbf{x}_i)-\mathbf{y}_i^\mathrm{o}) + J_\mathrm{c}(\mathbf{x}_0,\dots,\mathbf{x}_n) \tag{3.5}\] で表される。予報誤差共分散は、同化窓内において\(M\)の接線型演算子\(\mathbf{M}\)で時間発展し、解析値は式(Equation 3.3)の\(\mathbf{B}\)を\(\mathbf{M}^\mathrm{T}\mathbf{BM}\)で置き換えたものとなる。

4DVarでは、解析値を効率的に求めるために、\(J(\mathbf{x_0})\)の勾配を用いた最適化を用いる(Talagrand and Courtier 1987)。初期値\(\mathbf{x}_0\)に対する勾配 \[ \mathbf{p}_0 = \nabla_{\mathbf{x}_0}J(\mathbf{x}_0) \tag{3.6}\] を随伴方程式を用いて求める。随伴変数\(\mathbf{p}_{n+1} = \mathbf{0}\)を初期値として、時間逆向きに \[ \mathbf{p}_i = \mathbf{M}_i^\mathrm{T}\mathbf{p}_{i+1} + \frac{\partial J}{\partial\mathbf{x}_i}\;i=n,\dots,0 \tag{3.7}\] 右辺第一項は、接線型モデル\(\mathbf{M}_i^\mathrm{T}\)で一つ前の時間レベルまで、随伴変数を積分すること表している。以下に随伴モデルの作り方を示す。

3.3 随伴法の理論

随伴法の理論(Talagrand and Courtier 1987; Talagrand 1991)について述べる。

ヒルベルト空間\(\mathscr{F}\)での内積を\((,)\)で表す。 \(\mathbf{v}\rightarrow\mathscr{J}(\mathbf{v})\)は、\(\mathscr{F}\)上の微分可能な函数とする。 \(\mathscr{F}\)上の任意の点で\(\mathscr{F}\)の微分\(\delta\mathscr{F}\)は、内積により\(\mathbf{v}\)の微分\(\delta\mathbf{v}\)で表すことができる。

\[ \delta\mathscr{J} = (\nabla_\mathbf{v}\mathscr{F},\delta\mathbf{v}) \tag{3.8}\]

ここで\(\nabla_\mathbf{v}\mathscr{J}\)は、\(\mathbf{v}\)についての\(\mathscr{F}\)の勾配で唯一に定まる。 \(\mathscr{F}\)が有限次元で正規直交系\(v_i\)で表されるとき、\(\nabla_\mathbf{v}\mathscr{F}\)の要素は\(\partial\mathscr{J}/\partial v_i\)となるが、勾配の概念はより一般的なものを表している。

\(\mathscr{F}\)とは別のヒルベルト空間\(\mathscr{E}\)での内積を\(\langle,\rangle\)で表し、\(\mathbf{L}\)を\(\mathscr{E}\)から\(\mathscr{F}\)への連続な線型演算子とする。任意の\(\mathbf{u}\in\mathscr{E}\)と任意の\(\mathbf{v}\in\mathscr{F}\)に対して内積が等値である、\(\mathscr{F}\)から\(\mathscr{E}\)への連続な線型演算子の\(\mathbf{L}^*\)が唯一存在する。

\[ (\mathbf{v},\mathbf{Lu}) = \langle\mathbf{L}^*\mathbf{v},\mathbf{u}\rangle \tag{3.9}\]

\(\mathbf{L}^*\)は\(\mathbf{L}\)の随伴演算子（adjoint operator）と呼ばれる。 \(\mathscr{E}\)と\(\mathscr{F}\)が有限次元で正規直交系座標で表されるとき、\(\mathbf{L}^*\)は\(\mathbf{L}\)を表す行列の転置（要素が複素数なら共軛転置）である。

必ずしも線型とは限らない、\(\mathscr{E}\)から\(\mathscr{F}\)への可能な函数\(\mathbf{u}\rightarrow\mathbf{v} = \mathbf{G}(\mathbf{u})\)を考えると、\(\mathscr{J}(\mathbf{v}) = \mathscr{J}[\mathbf{G}(\mathbf{u})]\)は\(\mathbf{u}\)の合成函数となる。 \(\mathbf{u}\)の微分はEquation 3.8、\(\delta\mathbf{v}\)の微分は

\[ \delta\mathbf{v} = \mathbf{G}'\delta\mathbf{u} \tag{3.10}\]

で与えられる。ここで\(\mathbf{G}'\)は\(\mathbf{G}\)を微分した\(\mathscr{E}\)から\(\mathscr{F}\)への線型演算子である。 Equation 3.10をEquation 3.8に代入し、\(\mathbf{G}\)の随伴\(\mathbf{G}^{'*}\)を導入すると

\[ \delta\mathscr{J} = (\nabla_\mathbf{v}\mathscr{J},\mathbf{G}'\delta\mathbf{u}) = \langle\mathbf{G}^{'*}\nabla_\mathbf{v}\mathscr{J}, \delta\mathbf{u}\rangle \]

が得られる。勾配の定義により、\(\mathscr{J}\)の\(\mathbf{u}\)についての勾配\(\nabla_\mathbf{u}\mathscr{J}\)は

\[ \nabla_\mathbf{u} = \mathbf{G}^{'*}\nabla_\mathbf{v}\mathscr{J} \tag{3.11}\]

Equation 3.11は制御理論における随伴方程式の基本であり、数値的に\(\nabla_\mathbf{u}\mathscr{J}\)を計算する非常に効率的な方法を与える。ここでは、\(\mathbf{u}\rightarrow\mathbf{v} = \mathbf{G}(\mathbf{u})\)を想定しているので、\(\mathbf{v}\)は\(\mathbf{u}\)の明示的ではあるが複雑な関数になっているため、計算可能な\(\nabla_\mathbf{u}\mathscr{J}\)の解析的な式を得ることは困難である。 \(\nabla_\mathbf{u}\mathscr{J}\)の数値的に見積もるために、\(\mathbf{u}\)の各要素に摂動を加えて、それぞれの摂動に対して\(\mathbf{v} = \mathbf{G}(\mathbf{u})\)を計算するのは、大次元では現実的ではない。 Equation 3.11が示しているのは、与えられた\(\mathbf{w}\)に対して\(\mathbf{G}^{'*}\mathbf{w}\)を計算するプログラムがあれば、\(\nabla_\mathbf{u}\)は\(\nabla_\mathbf{v}\)から容易に計算できるということである。 \(\nabla_\mathbf{v}\)自体は\(\mathscr{J}\)が\(\mathbf{v}\)の簡単な函数であれば容易に求めることができる。て\(\mathbf{G}^{'*}\mathbf{w}\)の計算は、\(\mathbf{G}(\mathbf{u})\)を計算と同程度の複雑さであり、\(\mathbf{v} = \mathbf{G}(\mathbf{u})\)を計算して\(\nabla_\mathbf{v}\mathscr{J}\)を求め、Equation 3.11により最終的に\(\nabla_\mathbf{u}\)を求める計算は\(\mathbf{G}(\mathbf{u})\)の2〜3倍程度であり、直接明示的に\(\mathbf{u}\)の各要素に摂動を与えるよりも大幅に計算量を削減できる。

3.4 随伴モデル

次の形の常微分方程式を考える。

\[ \frac{\mathrm{d}\mathbf{w}}{\mathrm{d}t} = N(\mathbf{w}) + \mathbf{Lw} + \mathbf{f} \tag{3.12}\]

状態変数を\(\mathbf{w}\)、非線型項を\(N(\mathbf{w})\)、線型項を\(\mathbf{Lw}\)、強制を\(\mathbf{f}\)で表している。状態変数に摂動を与えた\(\mathbf{w}+\delta\mathbf{w}\)を(Equation 3.12)に代入し、摂動の時間発展\(\mathrm{d}\delta\mathbf{w}/\mathrm{d}t=\mathrm{d}(\mathbf{w}+\delta\mathbf{w})/\mathrm{d}t-\mathrm{d}\mathbf{w}/\mathrm{d}t\)を求めると接線型モデル（TLM: tangent linear model）が得られる。

\[ \frac{\mathrm{d}\delta\mathbf{w}}{\mathrm{d}t} = \left(\left.\frac{\partial N}{\partial\mathbf{w}}\right|_\mathbf{w} + \mathbf{L}\right)\delta\mathbf{w} \tag{3.13}\] ここで、摂動\(\delta\mathbf{w}\)が微小であると仮定し、その二次以上の項を無視した。強制項\(\mathbf{f}\)はTLMにはない。

\(\partial N/\partial\mathbf{w}|_\mathbf{w} + \mathbf{L}\)を改めて\(\mathbf{L}\)と置き、入力\(\delta\mathbf{w}\)及び出力\(\dot{\delta\mathbf{w}}=\mathrm{d}\delta\mathbf{w}/\mathrm{d}t\)を縦に並べたベクトルを用いて行列形式で表す。

\[ \begin{bmatrix} \delta\mathbf{w}\\ \dot{\delta\mathbf{w}} \end{bmatrix} = \begin{bmatrix} \mathbf{I} & \mathbf{0}\\ \mathbf{L} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \delta\mathbf{w}\\ \dot{\delta\mathbf{w}} \end{bmatrix} \tag{3.14}\] つまり \[ \begin{aligned} \delta\mathbf{w} &= \delta\mathbf{w}\\ \dot{\delta\mathbf{w}} &= \mathbf{L}\delta\mathbf{w} \end{aligned} \] である。随伴変数を\(\mathbf{w}^\mathrm{a}\)とすると、(Equation 3.14)の随伴モデル（ADM, adjoint model）は \[ \begin{bmatrix} \mathbf{w}^\mathrm{a}\\ \dot{\mathbf{w}}^\mathrm{a} \end{bmatrix} = \begin{bmatrix} \mathbf{I} & \mathbf{L}^\mathrm{T}\\ \mathbf{0} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{w}^\mathrm{a}\\ \dot{\mathbf{w}}^\mathrm{a} \end{bmatrix} \tag{3.15}\] と表すことができる。すなわち、 \[ \begin{aligned} \mathbf{w}^\mathrm{a} &= \mathbf{w}^\mathrm{a} + \mathbf{L}^\mathrm{T}\dot{\mathbf{w}}^\mathrm{a}\\ \dot{\mathbf{w}}^\mathrm{a} &= \mathbf{0} \end{aligned} \] である。 \(\dot{\mathbf{w}}^\mathrm{a} = \mathbf{0}\)もADMの一部であり、忘れるとバグの原因となりうる。

3.5 随伴モデルの作成

上述のように、TLMが行列\(\mathbf{L}\)で表せる場合は、ADMはその転置\(\mathbf{L}^\mathrm{T}\)を取ればよい。明示的に行列で表さなくても、TLMのソースコードを行毎に逆順にたどりながらADMを作成することもできる。ここでは、後で用いる、いくつかの簡単な操作に対する随伴を求める。

代入

代入\(A=B\)の接線型は \[ \begin{bmatrix} \delta B \\ \delta A \end{bmatrix}= \begin{bmatrix} 1 & 0 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} \delta B \\ \delta A \end{bmatrix} \tag{3.16}\] で、その随伴は \[ \begin{bmatrix} B^\mathrm{a} \\ A^\mathrm{a} \end{bmatrix}= \begin{bmatrix} 1 & 1 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} B^\mathrm{a} \\ A^\mathrm{a} \end{bmatrix} \tag{3.17}\] つまり \[ \begin{aligned} B^\mathrm{a} &= A^\mathrm{a} + B^\mathrm{a}\\ A^\mathrm{a} &=0 \end{aligned} \] である。

係数を掛けて和を取る操作

代入\(C = \alpha A + \beta B + \gamma D\)の接線型は \[ \begin{bmatrix} \delta A\\ \delta B \\ \delta D \\ \delta C \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ \alpha & \beta & \gamma & 0 \end{bmatrix} \begin{bmatrix} \delta A\\ \delta B \\ \delta D \\ \delta C \end{bmatrix} \tag{3.18}\] で、その随伴は \[ \begin{bmatrix} A^\mathrm{a} \\ B^\mathrm{a} \\ D^\mathrm{a} \\ C^\mathrm{a} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & \alpha\\ 0 & 1 & 0 & \beta\\ 0 & 0 & 1 & \gamma\\ 0 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} A^\mathrm{a} \\ B^\mathrm{a} \\ D^\mathrm{a} \\ C^\mathrm{a} \end{bmatrix} \tag{3.19}\] つまり \[ \begin{aligned} A^\mathrm{a} &= A^\mathrm{a} + \alpha C^\mathrm{a} \\ B^\mathrm{a} &= B^\mathrm{a} + \beta C^\mathrm{a} \\ D^\mathrm{a} &= D^\mathrm{a} + \gamma C^\mathrm{a} \\ C^\mathrm{a} &=0 \end{aligned} \] である。

ループ

TLMの最後の行から逆順にADMを作る。ループは添字を逆に回す。 \(i = 1, 2, \cdots, m,\,1\le j \le m\)に対して

\[ \delta x_j = c_1\delta x_1 + c_2\delta x_2 + \cdots + c_ix_i + \cdots + c_m\delta x_m \tag{3.20}\]

の随伴は

\[ \begin{aligned} x_m^\mathrm{a} &= x_m^\mathrm{a} + c_m x_j^\mathrm{a}\\ & \vdots \\ x_i^\mathrm{a} &= x_i^\mathrm{a} + c_i x_j^\mathrm{a}\\ & \vdots \\ x_2^\mathrm{a} &= x_2^\mathrm{a} + c_2 x_j^\mathrm{a}\\ x_1^\mathrm{a} &= x_1^\mathrm{a} + c_1 x_j^\mathrm{a}\\ x_j^\mathrm{a} &= c_j x_j^\mathrm{a} \end{aligned} \tag{3.21}\]

TLMの行の右辺に現れる摂動変数\(\delta x_i\)に対応して、ADMの行では随伴変数\(x_i^\mathrm{a}\)が左辺に現れる。 \(x_i^\mathrm{a}\)には\(x_i^\mathrm{a}\)自身に係数\(c_i\)とTLMの左辺に現れる摂動変数の随伴\(x_j^\mathrm{a}\)との積を加えたものになる。左辺の変数と右辺の変数が場所が入れ替わっている。 \(x_j^\mathrm{a}\)は最後に現れる。 \(c_j=0\)の場合のように、TLMの右辺に\(x_j^\mathrm{a}\)がない場合でも、\(x_j^\mathrm{a}\)に0を割り当てないと、他の部分でエラーを引き起こしうる。

ADMにおける変分演算は、微分演算とは異なり、ある時刻における状態ベクトル全体に対する摂動なので、状態ベクトルの一部が変化しないからといって、それを省略したものはADMとしては正しくない。

3.6 Lorenz-63モデル

Lorenz (1963) のモデル（以下Lorenz-63）は熱対流を理想化したモデルで，パラメタ次第でカオスにふるまう。Lorenz-63は次の3変数の常微分方程式で表される。

\[ \begin{aligned} \dot{X} &= -\sigma X + \sigma Y \\ \dot{Y} &= -XZ + rX -Y \\ \dot{Z} &= XY - \beta Z \end{aligned} \]

パラメタ\(\sigma\), \(r\), \(\beta\)はそれぞれPrandtl数，変形されたRayleigh数，アスペクト比を表す。予報変数\(X, Y, Z\)はそれぞれ無次元化された対流の強さ，最大温度差，対流に伴う成層の変化を表す。左辺の変数の上の\(\cdot\)は時間微分\(\mathrm{d}/\mathrm{d}t\)を表す記号である。

\(\mathbf{x}=(1,\,3,\,5)\)（青）及び\(\mathbf{x}=(1.1,\,3.3,\,5.5)\)（橙）から無次元時刻10まで積分した結果の3次元可視化結果を示す。回転や拡大縮小（右クリック）ができる。

Code

l63 <- function(t, w, p, r, b){
  c(
            -p * w[1] + p * w[2],
    (r - w[3]) * w[1] -     w[2],
                     w[1] * w[2] - b * w[3]    
  )
}

rk4 <- function(f, t, y, h, ...) {
  k1 <- f(t, y, ...)
  k2 <- f(t + 0.5 * h, y + 0.5 * h * k1, ...)
  k3 <- f(t + 0.5 * h, y + 0.5 * h * k2, ...)
  k4 <- f(t + h, y + h * k3, ...)
  y + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6
}

step.fom <- function(f, x, nstep, dt, ...) {
  t <- 0
  for (i in 1:nstep) {
    x <- rk4(f, t, x, dt, ...)
    t <- t + dt
  }
  x
}

Code

library(rgl)

options(rgl.useNULL = TRUE)
rgl::setupKnitr()
palette("Tableau 10")

ns <- 3
x0 <- c(1, 3, 5)
x0.p <- c(1.1, 3.3, 5.5)
p <- 10
r <- 32
b <- 8 / 3
dt <- 0.01
nt <- 1000

x.rk4 <- matrix(rep(0, ns * nt), nrow=ns)
x.rk4[, 1] <- x0
x.rk4p <- matrix(rep(0, ns * nt), nrow=ns)
x.rk4p[, 1] <- x0.p
for (j in 2:nt) {
  t <- (j - 1) * dt
  x.rk4[, j] <- step.fom(l63, x.rk4[, j-1], 1, dt, p, r, b)
  x.rk4p[, j] <- step.fom(l63, x.rk4p[, j-1], 1, dt, p, r, b) 
}

plot3d(x.rk4[1,], x.rk4[2,], x.rk4[3,],
  axes=TRUE, xlab="X", ylab="Y", zlab="Z", type='l',
  lwd=2, col=1)
plot3d(x.rk4p[1,], x.rk4p[2,], x.rk4p[3,], type="l",
  lwd=3, col=2, add=TRUE)
plot3d(x.rk4[1,1], x.rk4[2,1], x.rk4[3,1],
  type="s", size=1, col=1, add=TRUE)
plot3d(x.rk4p[1,1], x.rk4p[2,1], x.rk4[3,1],
  type="s", size=1, col=2, add=TRUE)
palette("R4")

Lorenz-63の随伴モデルの作成方法は Huang and Yang (1996) に詳述されている。この報告書やtenomoto/l63vdaでは時間積分にEuler法が用いられているが、4次のRunge-Kutta法（RK4）を用いると、より安定するので、パラメタ\(r\)を大きめにする。青い実線及び破線は、時間積分にRK4を用いて\(\mathbf{x}=(1, 3, 5)\)及び\((1.1, 3.3, 5.5)\)から積分したもので、橙色の実線及び破線はEuler法で\(\mathbf{x}=(1, 3, 5)\)及び\((1.1, 3.3, 5.5)\)から積分したものである。点線は報告書に掲載されている同じ設定とよく一致している。 RK4を用いると、解は安定していて初期値の差が現れるのが遅い。

Code

ns <- 3
x0 <- c(1, 3, 5)
x0.p <- c(1.1, 3.3, 5.5)
p <- 10
r <- 32
b <- 8 / 3
dt <- 0.01
nt <- 800

x.rk4 <- matrix(rep(0, ns * nt), nrow=ns)
x.rk4[, 1] <- x0
x.rk4p <- matrix(rep(0, ns * nt), nrow=ns)
x.rk4p[, 1] <- x0.p
x.eul <- matrix(rep(0, ns * nt), nrow=ns)
x.eul[, 1] <- x0
x.eulp <- matrix(rep(0, ns * nt), nrow=ns)
x.eulp[, 1] <- x0.p

for (j in 2:nt) {
  t <- (j - 1) * dt
  x.rk4[, j] <- step.fom(l63, x.rk4[, j-1], 1, dt, p, r, b)
  x.rk4p[, j] <- step.fom(l63, x.rk4p[, j-1], 1, dt, p, r, b) 
  x.eul[, j] <- x.eul[, j-1] + l63(t, x.eul[, j-1], p, r, b) * dt
  x.eulp[, j] <- x.eulp[, j-1] + l63(t, x.eulp[, j-1], p, r, b) * dt
}


title <- paste("L63 free run")
off <- c(0, 40, 60)
t <- (1:nt) * dt
x.rk4 <- x.rk4 + off
x.rk4p <- x.rk4p + off
x.eul <- x.eul + off
x.eulp <- x.eulp + off
plot(1, type = "n", xlab="time", ylab="state", main=title,
     xlim=c(min(t), max(t)), ylim=c(-20, 150))
for (i in 1:ns) {
  lines(t, x.rk4[i,], lwd=2, lty=1, col=1)
  lines(t, x.rk4p[i,], lwd=2, lty=2, col=1)
  lines(t, x.eul[i,], lwd=2, lty=1, col=2)
  lines(t, x.eulp[i,], lwd=2, lty=2, col=2)
  abline(h=off[i], lty=3, col=10)
}
legend("topright", lty=c(1, 2, 1, 2), ncol=2,
       legend=c("RK4", "RK4p", "Euler", "Eulerp"),
       col=c(1, 1, 2, 2))

Lorez-63の接線型モデルは次のように書ける。

\[ \begin{aligned} \dot{\delta X} &= -\sigma \delta X + \sigma \delta Y \\ \dot{\delta Y} &= -\delta XZ -X\delta Z + r\delta X - \delta Y \\ \dot{\delta Z} &= \delta XY + X\delta Y - b\delta Z \end{aligned} \] \[ \begin{bmatrix}\delta\dot{X}\\\delta\dot{Y}\\\delta\dot{Z}\end{bmatrix} = \begin{bmatrix} -\sigma & \sigma & 0\\ r-Z & -1 & -X\\ Y & X & -\beta \end{bmatrix} \begin{bmatrix}\delta X\\\delta Y\\\delta Z \end{bmatrix} \] 随伴モデルは行列を転置して次のように書ける。

\[ \begin{bmatrix}X^\mathrm{a}\\Y^\mathrm{a}\\Z^\mathrm{a}\end{bmatrix} = \begin{bmatrix}X^\mathrm{a}\\Y^\mathrm{a}\\Z^\mathrm{a}\end{bmatrix} +\begin{bmatrix} -\sigma & r-Z & Y\\ \sigma & -1 & X\\ 0 & -X & -\beta \end{bmatrix} \begin{bmatrix}\dot{X}^\mathrm{a}\\\dot{Y}^\mathrm{a}\\\dot{Z}^\mathrm{a} \end{bmatrix} \] ここで、\((X^\mathrm{a}, Y^\mathrm{a}, Z^\mathrm{a})^\mathrm{T}\)は随伴変数を表す。つまり \[ \begin{aligned} X^{\mathrm{a}} &= X^\mathrm{a} - \sigma\dot{\delta X}^\mathrm{a} + (r-Z)\dot{Y}^\mathrm{a} + Y\dot{Z}^\mathrm{a} \\ Y^{\mathrm{a}} &= Y^{\mathrm{a}} + \sigma\dot{X}^\mathrm{a} - \dot{ Y}^\mathrm{a} + X\dot{Z}^\mathrm{a} \\ Z^\mathrm{a} &= Z^\mathrm{a} - X\dot{Y}^\mathrm{a} - b\dot{Z}^\mathrm{a} \end{aligned} \] となる。

3.7 Lorenz-96モデル

Lorenz-96モデル(Lorenz and Emanuel 1998) \[ \frac{\mathrm{d}X_i}{\mathrm{d}t} = (X_{i+1} - X_{i-2})X_{i-1} - X_i + F \tag{3.22}\] のTLMとADMを作る。非線型項を\(N(X)=(X_{i+1} - X_{i-2})X_{i-1}\)を線型化するには、掛け算の一方を順に摂動に置き換えるか、 \[ \begin{aligned} \frac{\partial N}{\partial X_{i-2}} &= -X_{i-1} \\ \frac{\partial N}{\partial X_{i-1}} &= X_{i+1}-X_{i-2}\\ \frac{\partial N}{\partial X_{i+1}} &= X_{i-1} \end{aligned} \tag{3.23}\] のように微分を計算すれば、次のように求められる。 \[ \frac{\mathrm{d}\delta X_i}{\mathrm{d}t} = -X_{i-1}\delta X_{i-2} + (X_{i+1} - X_{i-2})\delta X_{i-1} + X_{i-1}\delta X_{i+1} - \delta X_i \tag{3.24}\]

(Equation 3.15)に基づいて、随伴モデルを作る。 (Equation 3.24)の右辺が表す\(\mathbf{L}\)の随伴\(\mathbf{L}^\mathrm{T}\)は、\(\delta X\)の添字が\(i\)となる行における係数を考えればよい。右辺第1、2、3項の\(\delta X_{i-2},\,\delta X_{i-1},\,\delta X_{i+1}\)なので、それぞれ\(+2,\,+1,\,-1\)加えれば\(i\)行目になる。さらに、左辺の変数（時間微分）と右辺の変数を入れ替え\(\mathbf{L}^\mathrm{T}\dot{\mathbf{w}}^\mathrm{a}\)が得られる。

L96の随伴モデルを書き下すと \[ X_i^\mathrm{a} = X_i^\mathrm{a} - X_{i+1}\frac{\mathrm{d}X_{i+2}^\mathrm{a}}{\mathrm{d}t} + (X_{i+2} -X_{i-1})\frac{\mathrm{d}X_{i+1}^\mathrm{a}}{\mathrm{d}t} +X_{i-2}\frac{\mathrm{d}X_{i-1}^\mathrm{a}}{\mathrm{d}t} - \frac{\mathrm{d}X_{i}^\mathrm{a}}{\mathrm{d}t}\\ \tag{3.25}\] 及び \[ \frac{\mathrm{d}X_{i}^\mathrm{a}}{\mathrm{d}t} = 0 \] となる。

3.8 拡散方程式

1次元の拡散方程式

\[ \frac{\partial u}{\partial t} = \sigma\frac{\partial^2 u}{\partial x^2} \] を差分で離散化して、接線型モデルと随伴モデルを導出する(Kalnay 2003)。

時刻（\(t\)、添字\(k\)）方向に前方、空間（\(x\)、添字\(i\)）方向に中央差分すると

\[ u_i^{k+1} = u_i^k + \alpha(u_{i+1}^k - 2u_i^k + u_{i-1}^k) \tag{3.26}\]

ここで\(\alpha = \sigma\Delta t/\Delta x^2\)である。

拡散方程式は線型なので、\(u = u_\mathrm{b} + \delta u\)のように基本場と摂動に分けると、接線型モデルはEquation 3.26の変数を摂動に置き換えて次のように書ける。

\[ \begin{aligned} \delta u_i^{k+1} &= \delta u_i^k + \alpha(\delta u_{i+1}^k - 2\delta u_i^k + \delta u_{i-1}^k)\\ &= \begin{pmatrix}\alpha & 1 - 2\alpha & \alpha\end{pmatrix} \begin{pmatrix}\delta u_{i-1}^k \\ \delta u_i^k \\ \delta u_{i+1}^k\end{pmatrix} \end{aligned} \tag{3.27}\]

上式に現れる変数を列ベクトルで表し、行列ベクトル積で表すと次のようになる。

\[ \begin{pmatrix}\delta u_{i-1}^k \\ \delta u_i^k \\ \delta u_{i+1}^k \\\delta u_i^{k+1}\end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ \alpha & 1 - 2\alpha & \alpha & 0 \end{pmatrix} \begin{pmatrix}\delta u_{i-1}^k \\ \delta u_i^k \\ \delta u_{i+1}^k \\\delta u_i^{k+1}\end{pmatrix} \tag{3.28}\]

摂動をベクトル\(\delta \mathbf{u}^k\)で表すと

\[ \begin{pmatrix}\delta\mathbf{u}^k \\ \delta\mathbf{u}^{k+1}\end{pmatrix} = \begin{pmatrix}\mathbf{I} & 0 \\ \mathbf{L} & 0\end{pmatrix} \begin{pmatrix}\delta\mathbf{u}^k \\ \delta\mathbf{u}^{k+1}\end{pmatrix} \]

と表すことができる。ここで\(\mathbf{I}\)は単位行列で\(\mathbf{L}\)は、境界条件に合わせた（準）三重対角行列になる。例えば周期境界条件の場合 \[ \mathbf{L} = \begin{pmatrix} 1 - 2\alpha & \alpha & 0 & \dots & 0 & 0 & \alpha\\ \alpha & 1- 2\alpha & \alpha & 0 & \dots & 0 &0\\ 0 & \alpha & 1-2\alpha & \alpha & 0 & \dots & 0\\ \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots\\ 0 & \dots & 0 & \alpha &1-2\alpha & \alpha &0 \\ 0 & \dots & 0 & 0 & \alpha & 1-2\alpha &\alpha \\ \alpha & 0 & 0 & \dots & 0 & \alpha & 1-2\alpha \end{pmatrix} \]

Equation 3.28を転置すると、随伴モデル \[ \begin{pmatrix}\delta^* u_{i-1}^k \\ \delta^* u_i^k \\ \delta^* u_{i+1}^k \\ \delta^* u_i^{k+1}\end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 & \alpha\\ 0 & 1 & 0 & 1-2\alpha\\ 0 & 0 & 1 & \alpha\\ 0 & 0 & 0 & 0 \end{pmatrix} \begin{pmatrix}\delta^* u_{i-1}^k \\ \delta^* u_i^k \\ \delta^* u_{i+1}^k \\ \delta^* u_i^{k+1}\end{pmatrix} \tag{3.29}\] あるいは

\[ \begin{pmatrix}\delta^*\mathbf{u}^k \\ \delta^*\mathbf{u}^{k+1}\end{pmatrix} = \begin{pmatrix}\mathbf{I} & \mathbf{L} \\ 0 & 0\end{pmatrix} \begin{pmatrix}\delta\mathbf{u}^k \\ \delta\mathbf{u}^{k+1}\end{pmatrix} \]

となる。この例では\(\mathbf{L}^* = \mathbf{L}\)である。

Equation 3.29を書き下すと

\[ \begin{aligned} \delta^*u_{i-1}^k &= \delta^*u_{i-1}^k + \alpha \delta^*u_i^{k+1} \\ \delta^*u_i^k &= \delta^*u_i^k + (1 - 2\alpha)\delta^*u_i^{k+1}\\ \delta^*u_{i+1}^k &= \delta^*u_{i+1}^k + \alpha \delta^*u_i^{k+1}\\ \delta^*u_i^{k+1} &= 0 \end{aligned} \]

となる。 \(k+1\)における摂動の更新（0にする）は最後に行う。

\(u\)のフーリエ変換と逆変換を

\[ \begin{aligned} u_m &= \frac{1}{\sqrt{2\pi}}\int_0^{2\pi}u(x)\exp(-imx)\mathrm{d}x\\ u(x) &= \frac{1}{\sqrt{2\pi}}\sum_{-\infty}^\infty u_m\exp(imx) \end{aligned} \] とすると、2階微分は次のように書ける。

\[ \frac{\partial^2 u}{\partial x^2} = -\frac{m^2}{\sqrt{2\pi}}\sum_{-\infty}^\infty u_m\exp(imx) \] 拡散方程式は波数空間で次のように表される。

\[ \mathbf{u}^{k+1} = \alpha \mathbf{u}^k \] ここで、\(\alpha = 1 - \sigma m^2\Delta t\)である。

接線型モデルは \[ \delta\mathbf{u}^{k+1} = \alpha\delta\mathbf{u}^k \]

\[ \begin{pmatrix} \delta\mathbf{u}^k \\ \delta\mathbf{u}^{k+1} \end{pmatrix} = \begin{pmatrix} \mathbf{I} & 0 \\ \mathbf{L} & 0 \end{pmatrix} \begin{pmatrix} \delta\mathbf{u}^k \\ \delta\mathbf{u}^{k+1} \end{pmatrix} \]

となる。ここで\(\mathbf{L} = \alpha\mathbf{I}\)である。

随伴モデルは \[ \begin{pmatrix} \delta^*\mathbf{u}^k \\ \delta^*\mathbf{u}^{k+1} \end{pmatrix} = \begin{pmatrix} \mathbf{I} & \mathbf{L} \\ 0 & 0 \end{pmatrix} \begin{pmatrix} \delta^*\mathbf{u}^k \\ \delta^*\mathbf{u}^{k+1} \end{pmatrix} \] すなわち

\[ \begin{aligned} \delta^*\mathbf{u}^k &= (\mathbf{I} + \mathbf{L})\delta^*\mathbf{u}^{k+1} = (1 + \alpha)\mathbf{I}\delta^*\mathbf{u}^{k+1}\\ \delta^*\mathbf{u}^{k+1} &= 0 \end{aligned} \] となる。

3.9 4次のRunge-Kutta法

\[ \dot{y} = f(t, y) \tag{3.30}\]

を4次のRunge-Kutta（RK4）法で積分する。 \[ \begin{aligned} y_{n+1} &= y_{n} + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \\ k_1 &= f(t_n,\,y_n) \\ k_2 &= f\left(t_n + \frac{h}{2},\,y_n + \frac{h}{2}k_1\right) \\ k_3 &= f\left(t_n + \frac{h}{2},\,y_n + \frac{h}{2}k_2\right) \\ k_4 &= f(t_n + h,\,y_n + hk_3) \\ \end{aligned} \tag{3.31}\]

\(k\)の係数はButcherの表(Butcher and Wanner 1996) \[ \begin{array}{c|c} \mathbf{c} &\mathbf{A}\\ \hline & \mathbf{b}^\mathrm{T} \end{array} \tag{3.32}\] に整理できる。 \(s\)次の陽的Runge-Kutta法は \[ \begin{aligned} y_{n+1} &= y_n + h\sum_{i=1}^sb_ik_i\\ k_i &=f(t_n+c_ih, y_n+ h\sum_{j=1}^{i-1}a_{ij}k_j) \end{aligned} \] と書ける。 \(\mathbf{A}\)は対角成分が0である下三角行列になる。

RK4に対するButcherの表は次のようになる。 \[ \begin{array}{c|ccccc} 0 \\ 1/2 & 1/2\\ 1/2 & 0 & 1/2\\ 1 & 0 & 0 & 1\\ \hline & 1/6 & 1/3 & 1/3 & 1/6 \end{array} \tag{3.33}\]

RK4の接線型は次の通りである。 \[ \begin{aligned} \delta y_{n+1} &= \delta y_{n} + \frac{h}{6}(\delta k_1 + 2\delta k_2 + 2\delta k_3 + \delta k_4) \\ \delta k_1 &= \delta f(t_n,\,y_n,\,\delta y_n) \\ \delta k_2 &= \delta f\left(t_n + \frac{h}{2},\,y_n + \frac{h}{2}k_1,\,\delta y_n + \frac{h}{2}\delta k_1\right) \\ \delta k_3 &= \delta f\left(t_n + \frac{h}{2},\,y_n + \frac{h}{2}k_2,\,\delta y_n + \frac{h}{2}\delta k_2\right) \\ \delta k_4 &= \delta f(t_n + h,\,y_n + hk_3,\,\delta y_n + h\delta k_3) \end{aligned} \tag{3.34}\] ここで\(\delta f=\partial f/\partial y\)は接線型モデルである。 \(\delta k_1,\,\delta k_2,\,\delta k_3,\,\delta k_4\)は出力しないので、\(\delta y\)を更新した後 \[ \delta\mathbf{k} \equiv [\delta k_1,\,\delta k_2,\,\delta k_3,\,\delta k_4]^\mathrm{T}=0 \] と置くと、行列形式では \[ \begin{bmatrix} \delta\mathbf{k} \\ \delta y \end{bmatrix} = \begin{bmatrix} 0 & 0 \\ h\mathbf{b}^\mathrm{T} & 1 \end{bmatrix} \begin{bmatrix} \delta\mathbf{k} \\ \delta y \end{bmatrix} \tag{3.35}\] となる。

随伴は \[ \begin{bmatrix} \mathbf{k}^\mathrm{a} \\ y^\mathrm{a} \end{bmatrix} = \begin{bmatrix} 0 & h\mathbf{b}\\ 0 & 1 \end{bmatrix} \begin{bmatrix} \mathbf{k}^\mathrm{a} \\ y^\mathrm{a} \end{bmatrix} \tag{3.36}\] つまり、\(k_1^\mathrm{a} = hy^\mathrm{a}/6,\,k_2^\mathrm{a} = hy^\mathrm{a}/3,\,k_3^\mathrm{a} = hy^\mathrm{a}/3,\,k_4^\mathrm{a} = hy^\mathrm{a}/6\)と初期化される。

\(L\)を線型演算子、\(\alpha\)をスカラ係数として \(C=L(A+\alpha B)\)は \[ \begin{bmatrix} A \\ B \\ C \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ L & \alpha L & 0 \end{bmatrix} \begin{bmatrix} A \\ B \\ C \end{bmatrix} \] と書けるので、その随伴は \[ \begin{bmatrix} A^\mathrm{a} \\ B^\mathrm{a} \\ C^\mathrm{a} \end{bmatrix} = \begin{bmatrix} 1 & 0 & L^\mathrm{T}\\ 0 & 1 & \alpha L^\mathrm{T}\\ 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} A^\mathrm{a} \\ B^\mathrm{a} \\ C^\mathrm{a} \end{bmatrix} \] つまり\(A^\mathrm{a} = A^\mathrm{a} + LC^\mathrm{a},\,B^\mathrm{a} = B^\mathrm{a} + \alpha LC^\mathrm{a},\,C^\mathrm{a}=0\) である。これを用いるとRK4の随伴は次の順序で計算される。 \[ \begin{aligned} f_4^\mathrm{a} &\equiv f^\mathrm{a}(t_n + h,\,y_n + hk_3,\,k_4^\mathrm{a})\\ y^\mathrm{a} &= y^\mathrm{a} + f_4^\mathrm{a}\\ k_3^\mathrm{a} &= k_3^\mathrm{a} + hf_4^\mathrm{a} \\ k_4^\mathrm{a} &= 0 \\ f_3^\mathrm{a} &\equiv f^\mathrm{a}\left(t_n + \frac{h}{2},\,y_n + \frac{h}{2}k_2, \,k_3^\mathrm{a}\right)\\ y^\mathrm{a} &= y^\mathrm{a} + f_3^\mathrm{a}\\ k_2^\mathrm{a} &= k_2^\mathrm{a} + \frac{h}{2}f_3^\mathrm{a} \\ k_3^\mathrm{a} &= 0 \\ f_2^\mathrm{a} &\equiv \delta f^\mathrm{a}\left(t_n + \frac{h}{2},\,y_n + \frac{h}{2}k_1,\, k_2^\mathrm{a}\right)\\ y^\mathrm{a} &= y^\mathrm{a} + f_2^\mathrm{a}\\ k_1^\mathrm{a} &= k_1^\mathrm{a} + \frac{h}{2}f_2^\mathrm{a} \\ k_2^\mathrm{a} &= 0 \\ y^\mathrm{a} &= y^\mathrm{a} + f^\mathrm{a}(t_n,\,y_n,\,k_1^\mathrm{a})\\ k_1^\mathrm{a} &= 0 \end{aligned} \]

3.10 接線型及び随伴の動作確認

接線型は \[N(\mathbf{w}+\Delta\mathbf{w})-N(\mathbf{w})\approx \mathbf{L}(\mathbf{w})\Delta\mathbf{w} \tag{3.37}\] が近似的に成り立つことを確認する。ここで\(N\)は非線型モデル、\(\mathbf{L}\)は接線型モデル、\(\mathbf{w}\)は状態、\(\Delta\mathbf{w}\)は摂動を表す。精度は\(\mathbf{w}\)における\(N\)の非線型の強さや摂動\(\Delta\mathbf{w}\)の振幅に依存する。

随伴は \[(\mathbf{L}(\mathbf{w})\Delta\mathbf{w})^\mathrm{T}(\mathbf{L}(\mathbf{w})\Delta\mathbf{w})=\Delta\mathbf{w}^\mathrm{T}[\mathbf{L}(\mathbf{w})^\mathrm{T}(\mathbf{L}(\mathbf{w})\Delta\mathbf{w})] \tag{3.38}\] が成り立つことを確認する。厳密に成り立つか、機械精度以下の誤差（下1〜2桁の差異）でなければバグが存在していることを示す。