概要

　複雑な（またはカオス的な）偏微分方程式を解くと、初期値の設定によって予測結果が大きく異なる場合があります。昔から気象などの分野では、天気予報の精度を上げるため、観測データを使い、データ同化の手法（変分法）を用いて数値モデルの修正を行なっています。

　近年では、IoTやビックデータ活用などによりデータ同化が注目を集めています。データ同化は、大きく２つの手法があり、カルマンフィルタと変分法があります。

　変分法の特徴は、数値モデルが観測結果と一致するように、初期値や境界条件を修正する方法です。今回の記事では、4次元変分法について紹介します。

前回の記事では、3次元変分法について書きました。

４次元変分法と３次元変分法との違いは、方程式が時間発展するかしないかの違いです。

４次元変分法の評価関数の定義と勾配計算

　時間発展するモデルを以下のように定義する。

\begin{align}
x_t = M(x_{t-1})
\end{align}

$M$ は演算子であり、$x_t$ は時刻 $t$ ステップにおける状態変数を表す。また、$x_t$ はベクトルである。変分法では、基本的に初期値を求めたいので、初期時刻を陽に表した

\begin{align}
x_t = M^t(x_{0})
\end{align}

の形で書くこともできる。$M^t$ は演算子 $M$ を $t$ 回作用させたことを意味する。

　3次元変分法では、時間発展しないので評価関数は

\begin{align}
J(x_0)&= \frac{1}{2} \left(x_0-x^b \right)^T B^{-1} \left(x_0-x^b \right) + \frac{1}{2}  \left(H(x_0)-y_0 \right)^T R^{-1} \left(H(x_0)-y_0 \right)
\end{align}

であった。$y$は観測値であり、$H$は観測演算子、$B$は背景誤差共分散行列、$R$は観測誤差共分散行列である。$x^b$は、$x_0$の事前情報から推定される値である。

　４次元変分法では、これを拡張し評価関数を、初期値を陽に表して

\begin{align}
J(x_0)&= \frac{1}{2} \left(x_0-x^b \right)^T B^{-1} \left(x_0-x^b \right) + \frac{1}{2} \sum_{t=0}^T \left(H(x_t)-y_t \right)^T R_t^{-1} \left(H(x_t)-y_t \right)\\
&=\frac{1}{2} \left(x_0-x^b \right)^T B^{-1} \left(x_0-x^b \right) + \frac{1}{2} \sum_{t=0}^T \left(H(M^t(x_{0}))-y_t \right)^T R_t^{-1} \left(H(M^t(x_{0}))-y_t \right)\\
\end{align}

とする。下付き添字の$t$ は各時間ステップにおける値および演算子を示す。

　次に評価関数を初期値 $x_0$で微分する。まず、以下の量を定義しておく。

\begin{align}
\mathcal{H}_t \equiv \frac{\partial H_t}{\partial x_t} \ \ , \ \ \mathcal{M}_{t} \equiv  \frac{\partial  M(x_{t})}{\partial x_{t}} 
\end{align}

観測演算子 $H_t$の $x_0$ 微分は連鎖則を使い

\begin{align}
 \frac{\partial H_t}{\partial x_0} &= \frac{\partial H_t}{\partial x_t}\frac{\partial x_t}{\partial x_0} = \frac{\partial H_t}{\partial x_t}\frac{\partial  M(x_{t-1})}{\partial x_0} \\
&=\frac{\partial H_t}{\partial x_t}\frac{\partial  M(x_{t-1})}{\partial x_{t-1}} \cdots \frac{\partial  M(x_{0})}{\partial x_{0}}\\
& =\mathcal{H}_t\mathcal{M}_{t-1} \cdots \mathcal{M}_{0}
\end{align}

と計算できる。よって、評価関数の$x_0$の微分、つまり勾配は

\begin{align}
g \equiv \frac{\partial J}{\partial x_0}&=B^{-1} \left(x_0-x^b \right) + \sum_{t=0}^T \frac{\partial H_t^T}{\partial x_0} R_t^{-1} \left(H(M^t(x_{0}))-y_t \right)\\
&=B^{-1} \left(x_0-x^b \right) + \sum_{t=0}^T \mathcal{M}_{0}^T \cdots  \mathcal{M}_{t-1}^T \mathcal{H}_t^T R_t^{-1} \left(H(M^t(x_{0}))-y_t \right)\\
\end{align}

と計算される。時間を遡って勾配を求めて$x_0$ を修正する方法は、深層学習で使われる誤差逆伝播法と同様な計算法である。

　数値計算においては、上式で計算された勾配を用いて

\begin{align}
x_0 \leftarrow x_0 - \alpha g
\end{align}

と更新する。$\alpha$ はステップ幅である。

　実際の数値計算では、単純な勾配法ではなく、準ニュートン法を使う場合が散見される。またステップ幅は、直線探索法を使う場合が多い。

離散的なアジョイント方程式

　上記で導出した評価関数の勾配は、アジョイント方程式と呼ばれる、時間を逆に発展させた方程式を使い求めることができる。

　ラグランジュ未定乗数法を使い、条件付き最適化問題を解くことによりアジョイント方程式を導出する。

　時間発展するモデルは

\begin{align}
x_t = M(x_{t-1})
\end{align}

であり、評価関数を

\begin{align}
J(x_0,x_1,\ldots,x_T)&= \frac{1}{2} \left(x_0-x^b \right)^T B^{-1} \left(x_0-x^b \right) + \frac{1}{2} \sum_{t=0}^T \left(H(x_t)-y_t \right)^T R_t^{-1} \left(H(x_t)-y_t \right)
\end{align}

と書くことにする。

　ラグランジュ関数を、ラグランジュ未定乗数 $\lambda_t$ を使い

\begin{align}
\mathcal{L} = J(x_0,x_1,\ldots,x_T) + \sum_{t=1}^T \left\langle \lambda_t,M(x_{t-1})-x_t \right\rangle
\end{align}

と定義する。ここでブラケット$\langle \cdot, \cdot \rangle$は、内積を意味する。

　また、$x_t$ は演算子$M$ と初期値 $x_0$ で表せることから

\begin{align}
J(x_0) = J(x_0,x_1,\ldots,x_T)
\end{align}

と書けることを使うと、時間発展のモデルが成立するならば

\begin{align}
\mathcal{L}-J(x_0)=0
\end{align}

が成立する。

ラグランジュ関数を変分すると

\begin{align}
\delta\mathcal{L} &=  \sum_{t=0}^T \left\langle \frac{\partial J}{\partial x_t}\delta x_t \right\rangle + \sum_{t=1}^T \left\langle \delta\lambda_t,M(x_{t-1})-x_t \right\rangle + \sum_{t=1}^T \left\langle \lambda_t, \frac{\partial M(x_{t-1})}{\partial x_{t-1}} \delta x_{t-1} - \delta x_t \right\rangle \\
&= \sum_{t=0}^T \left\langle \frac{\partial J}{\partial x_t}\delta x_t \right\rangle + \sum_{t=1}^T \left\langle \delta\lambda_t,M(x_{t-1})-x_t \right\rangle + \sum_{t=1}^T \left\langle \mathcal{M}_{t-1}^{\ast} \lambda_t, \delta x_{t-1} \right\rangle - \sum_{t=1}^T \left\langle \lambda_t, \delta x_t \right\rangle \\
\end{align}

となる。行列$\mathcal{M}_t$ は、アジョイント行列と呼ばれ、ベクトル $a,b$ に対して

\begin{align}
\langle a,\mathcal{M}_t b \rangle = \langle \mathcal{M}_t^* a, b \rangle
\end{align}

が成立することを仮定している。（実際、成り立つ。）

　$\delta x_t$ を $\delta x_0,\delta x_T$ とそれ以外に比例する部分を分け、式を整理すると

\begin{align}
\delta\mathcal{L} &= \left\langle \frac{\partial J}{\partial x_0}+ \mathcal{M}_0^* \lambda_1 ,\delta x_0  \right\rangle 
+ \left\langle \frac{\partial J}{\partial x_T} -\lambda_T ,\delta x_T  \right\rangle \\
&+ \sum_{t=1}^{T-1} \left\langle \frac{\partial J}{\partial x_t}+ \mathcal{M}_t^* \lambda_{t+1} -\lambda_t ,\delta x_t  \right\rangle  
+ \sum_{t=1}^T \left\langle \delta\lambda_t,M(x_{t-1})-x_t \right\rangle 
\end{align}

となる。また、$\mathcal{L}-J(x_0)=0$ なので、 $\delta (\mathcal{L}-J(x_0))=0$ も成立する。つまり、

\begin{align}
\delta\mathcal{L}-\delta J(x_0)= \delta\mathcal{L}- \left\langle \frac{\partial J}{\partial x_0} ,\delta x_0  \right\rangle  = \delta\mathcal{L}- \left\langle g ,\delta x_0  \right\rangle  = 0
\end{align}

となる。これが成立するには、任意の $\delta x_t,\delta\lambda_t$ で成立しなければならないので

\begin{align}
&\frac{\partial J}{\partial x_0}+ \mathcal{M}_0^* \lambda_1 - g = 0 \\
&\frac{\partial J}{\partial x_t}+ \mathcal{M}_t^* \lambda_{t+1} -\lambda_t = 0 \ \ \ \ (t=1,2,\ldots,T-1) \\
&\frac{\partial J}{\partial x_T} -\lambda_T=0 \\
&M(x_{t-1})-x_t = 0 
\end{align}

となり方程式が得られる。最後の方程式は、時間発展モデルの方程式である。ここで、発展方程式の初期値 $\lambda_{T+1}=0$ とすれば、アジョイント方程式

\begin{align}
&\lambda_t = \mathcal{M}_t^* \lambda_{t+1} + \frac{\partial J}{\partial x_t}  \ \ \ \ (t=0,1,\ldots,T) \\
& \lambda_{T+1}=0 \ \ , \ \ g=\lambda_0 
\end{align}

が得られる。評価関数を代入すれば、

\begin{align}
\lambda_t &= \mathcal{M}_t^T \lambda_{t+1} + \mathcal{H}_t^T R_t^{-1} \left(H(x_{t})-y_t \right)
  \ \ \ \ (t=1,2,\ldots,T) \\
\lambda_0 &= \mathcal{M}_0^T \lambda_{1} +B^{-1} \left(x_0-x^b \right) +  \mathcal{H}_0^T R_0^{-1} \left(H(x_{0})-y_0 \right)  \\
g&=\lambda_0 
\end{align}

となる。行列は実行列としているので、複素共役は転置となる。
　
上記の式の $\lambda_0$ を出発として、 $\lambda_t$ を次々代入していけば、最初に登場した勾配

\begin{align}
g =B^{-1} \left(x_0-x^b \right) + \sum_{t=0}^T \mathcal{M}_{0}^T \cdots  \mathcal{M}_{t-1}^T \mathcal{H}_t^T R_t^{-1} \left(H(M^t(x_{0}))-y_t \right)\\
\end{align}

であることがわかる。時間発展するモデルを解いた後、アジョイント方程式を時間を遡り計算していくことで、勾配を求めることができる。

例：強制振動の方程式のアジョイントコード

　
　比較的簡単な強制振動の方程式のアジョイントコードについて説明する。強制振動の支配方程式は、

\begin{align}
\frac{\partial x}{\partial t}&=v \\
m\frac{\partial v}{\partial t}&=-kx -rv+\omega
\end{align}

であり、$m$は質量、$k$はバネ定数、$r$は減衰定数、$\omega$は外力である。離散化すれば、

\begin{align}
x_t &= x_{t-1}+\Delta t v_{t-1} \\
v_t &= -\frac{k\Delta t}{m} x_{t-1} + \left(1- \frac{r\Delta t}{m} \right)v_{t-1} +\frac{\Delta t}{m}\omega_{t-1}
\end{align}

と書ける。
　
　推定したい値は、初期値$x_0$と外力$\omega_t$とする。よって勾配は、$T+1$ 次元ベクトルになる。

　ラグランジュ関数は、

\begin{align}
\mathcal{L} =& J + \sum_{t=1}^{T}\lambda_t^{x}\left(x_{t-1}+\Delta t v_{t-1}- x_t\right) \\
&+\sum_{t=1}^{T}\lambda_t^{v}\left(-\frac{k\Delta t}{m} x_{t-1} + \left(1- \frac{r\Delta t}{m} \right)v_{t-1} +\frac{\Delta t}{m}\omega_{t-1} - v_t\right) \\
\end{align}

と定義する。　
　ラグランジュ関数の変分は、

\begin{align}
\delta\mathcal{L} =& \sum_{t=0}^{T}\frac{\partial J}{\partial x_t} \delta x_t + \sum_{t=0}^{T-1}\frac{\partial J}{\partial \omega_t} \delta \omega_t \\
&+\sum_{t=1}^{T}\delta\lambda_t^{x}\left(x_{t-1}+\Delta t v_{t-1}- x_t\right) \\
&+\sum_{t=1}^{T}\lambda_t^{x}\left(\delta x_{t-1}+\Delta t \delta v_{t-1}- \delta x_t\right) \\
&+\sum_{t=1}^{T}\delta\lambda_t^{v}\left(-\frac{k\Delta t}{m} x_{t-1} + \left(1- \frac{r\Delta t}{m} \right)v_{t-1} +\frac{\Delta t}{m}\omega_{t-1} - v_t\right) \\
&+\sum_{t=1}^{T}\lambda_t^{v}\left(-\frac{k\Delta t}{m} \delta x_{t-1} + \left(1- \frac{r\Delta t}{m} \right)\delta v_{t-1} +\frac{\Delta t}{m}\delta\omega_{t-1} - \delta v_t\right) \\
\end{align}

となる。また、

\begin{align}
\delta\mathcal{L}-\delta J(x_0,\omega_0,\omega_1,\ldots,\omega_{T-1})  = \delta\mathcal{L}- g_{x_0}\delta x_0 -\sum_{t=0}^{T-1} g_{\omega_t}\delta \omega_t  = 0
\end{align}

であり、任意の$\delta x_t ,\delta v_t,\delta \omega_t$ で成り立つので、アジョイント方程式

\begin{align}
\lambda_t^x &=  \lambda_{t+1}^x -\frac{k\Delta t}{m} \lambda_{t+1}^{v}+\frac{\partial J}{\partial x_t} \\
\lambda_t^v &= \Delta t\lambda_{t+1}^x + \left(1- \frac{r\Delta t}{m} \right)\lambda_{t+1}^{v}  \\
\lambda_t^{\omega} &= \frac{\Delta t}{m}\lambda_{t+1}^{v} + \frac{\partial J}{\partial \omega_t}\\
g_{\omega_t}&=\lambda_t^{\omega}\\
g_{x_0}&=\lambda_0^x\\

\end{align}

が得られる。（$\delta\lambda$からは支配方程式が得られる。）また、アジョイント方程式の初期値は、

\begin{align}
\lambda_{T+1}^x=0 \ \ ,\ \ \lambda_{T+1}^v=0 
\end{align}

である。行列を

\mathcal{M}=
\begin{pmatrix}
1                    & \Delta t & 0\\\ 
\frac{k\Delta t}{m}  & \left(1- \frac{r\Delta t}{m} \right)&  \frac{\Delta t}{m} \\\
0 & 0& 0
\end{pmatrix}

とおくと、支配方程式は、

\begin{pmatrix}
x_{t} \\\ 
v_{t}\\\
\omega_{t}
\end{pmatrix}=
\mathcal{M}
\begin{pmatrix}
x_{t-1} \\\ 
v_{t-1}\\\
\omega_{t-1}
\end{pmatrix}

と書ける。一方、アジョイント方程式は、評価関数の微分の部分は抜いて

\begin{pmatrix}
\lambda^x_{t}\\\ 
\lambda^v_{t}\\\
\lambda^{\omega}_{t}
\end{pmatrix}=
\mathcal{M}^{\ast}
\begin{pmatrix}
\lambda^x_{t+1} \\\ 
\lambda^v_{t+1}\\\
\lambda^{\omega}_{t+1}
\end{pmatrix}

と書ける。よって、行列$\mathcal{M}$ はアジョイント行列になっている。

　外力を

\begin{align}
\omega_t=\sin\left(\frac{2\pi t}{10} \right)
\end{align}

として、支配方程式の計算を以下のようなコードで書く。

! 真値
time = 0.0
do t =1,tt
  time = time + deltaT
  w(t-1) = sin(2.0*pi*time/10.0)
  x(t) = x(t-1) + deltaT*v(t-1) 
  v(t) = -(k*deltaT/m)*x(t-1) + (1.0 - (r*deltaT/m))*v(t-1) + (deltaT/m)*w(t-1)

enddo

! シミュレーションの解
time = 0.0

do t =1,tt
  time = time + deltaT
! 外力に乱数を加える。　　本当は、ガウス分布に従うようにする。
  call random_number(rnd)
  w_sim(t-1) = sin(2.0*pi*time/10.0) + rnd*0.1
  x_sim(t) = x_sim(t-1) + deltaT*v_sim(t-1) 
  v_sim(t) = -(k*deltaT/m)*x_sim(t-1) + (1.0 - (r*deltaT/m))*v_sim(t-1) + (deltaT/m)*w_sim(t-1)

enddo

観測値$y$ は偶数ステップ時に得られるとすると、評価関数は

\begin{align}
J = \frac{1}{2} \frac{(x_o-x^b)^2}{\sigma_b^2}+\frac{1}{2}\sum_{t=0}^T
 \frac{(\omega_t-\omega_t^b)^2}{\sigma_{\omega}^2}+\frac{1}{2}\sum_{m=1}^{T/2}
 \frac{(x_{2m}-y_{2m})^2}{\sigma_{o}^2}
\end{align}

と書ける。アジョイントコードは、以下のように書ける。


  ! アジョイントコード
  Ax = 0.0
  Av = 0.0

  do t = tt ,1 ,-1
    Ax_ = Ax - (k*deltaT/m)*Av
    if( mod(t,2) == 0 ) then
      Ax_ = Ax_ + (x(t) - x_sim(t))/(sigma_o**2)
    endif

    Av_ = deltaT*Ax + (1.0 - (r*deltaT/m))*Av
    g(t) = (deltaT/m)*Av + (w(t-1)-w_sim(t-1))/(sigma_w**2)

    Ax = Ax_
    Av = Av_
  enddo

  g(0) = Ax - (k*deltaT/m)*Av + (x(0) - x_sim(0))/(sigma_b**2)
enddo

連続的なアジョイント方程式

　時間発展する微分方程式を数値的に解く場合、（離散化されるので）上記で説明した離散的なアジョイント方程式を解き、最適化（ニュートン法）を行えばよいが、連続的なアジョイント方程式も導出しておく。

　初期値が決まれば、以下の微分方程式によって、$x$ の時間発展が計算できるとする。

\begin{align}
\frac{\partial x}{\partial t} = M(x)
\end{align}

ラグランジュ関数を、評価関数を用いて

\begin{align}
\mathcal{L} &= J + \int_0^T \mathrm{d}t \left\langle \lambda, M(x)-\frac{\partial x}{\partial t}\right\rangle \\
&J=  \int_0^T \mathrm{d}t \mathcal{J} 
\end{align}

と定義する。ラグランジュ関数の変分を行い、

\begin{align}
\delta\mathcal{L} = \int_0^T \mathrm{d}t &\left\{  \left\langle \frac{\partial \mathcal{J}}{\partial x},\delta x \right\rangle
+ \left\langle \delta\lambda, M(x)-\frac{\partial x}{\partial t}\right\rangle
\right. \\
& \left.
+ \left\langle \lambda,\frac{\partial M}{\partial x} \delta x \right\rangle 
- \left\langle \lambda, \frac{\partial }{\partial t}\delta x\right\rangle
\right\}
\end{align}

時間微分に関して、部分積分を行と

\begin{align}
\delta\mathcal{L} &= \int_0^T \mathrm{d}t  \left\langle \lambda, \frac{\partial }{\partial t}\delta x\right\rangle \\
&=\left\langle \lambda(T), \delta x(T) \right\rangle -\left\langle \lambda(0),\delta x(0)\right\rangle -\int_0^T \mathrm{d}t \left\langle \frac{\partial \lambda}{\partial t}, \delta x \right\rangle \\
\end{align}

ラグランジュ関数の変分は、アジョイント行列の性質を使い

\begin{align}
\delta\mathcal{L} &= \int_0^T \mathrm{d}t \left\langle \frac{\partial \mathcal{J}}{\partial x},\delta x\right\rangle
+\int_0^T \mathrm{d}t \left\langle \delta\lambda, M(x)-\frac{\partial x}{\partial t}\right\rangle　
+\int_0^T \mathrm{d}t \left\langle \frac{\partial M^{\ast}}{\partial x}\lambda, \delta x \right\rangle \\
&-\left\langle \lambda(T), \delta x(T) \right\rangle +\left\langle \lambda(0),\delta x(0)\right\rangle + \int_0^T \mathrm{d}t \left\langle \frac{\partial \lambda}{\partial t}, \delta x \right\rangle \\
\end{align}

と書ける。評価関数は、初期値$x(0)$ の関数なので

\begin{align}
&\delta\mathcal{L}-\delta J(x(0))= \delta\mathcal{L}- \left\langle \frac{\partial J}{\partial x(0)} ,\delta x(0) \right\rangle  = \delta\mathcal{L}- \left\langle g ,\delta x(0) \right\rangle  = 0\\
\end{align}

が成り立ち、以下のようにまとめられる。

\begin{align}
\delta\mathcal{L}-\delta J(x(0))&=\int_0^T \mathrm{d}t \left\langle \frac{\partial \mathcal{J}}{\partial x}+\frac{\partial M^{\ast}}{\partial x}\lambda+\frac{\partial \lambda}{\partial t},\delta x\right\rangle
+\int_0^T \mathrm{d}t \left\langle \delta\lambda, M(x)-\frac{\partial x}{\partial t}\right\rangle　\\
&-\left\langle \lambda(T), \delta x(T) \right\rangle +\left\langle \lambda(0)-g,\delta x(0)\right\rangle 
\end{align}

　任意の $\delta x,\delta \lambda$ で成立しなければならない。よって、連続的なアジョイント方程式と時間発展を記述する微分方程式

\begin{align}
& \lambda(T) = 0 \\
&\frac{\partial \lambda}{\partial t} + \frac{\partial M^{\ast}}{\partial x}\lambda + \frac{\partial \mathcal{J}}{\partial x} = 0 \\
& g =\lambda(0) \\
& \frac{\partial x}{\partial t}=M(x)
\end{align}

が得られる。

例：トレーサーの初期状態の推定

　トレーサー $C(t,x)$ の時空間分布は、１次元線形移流拡散方程式

\begin{align}
\frac{\partial C}{\partial t} = -u \frac{\partial C}{\partial x} +\nu \frac{\partial^2 C}{\partial x^2} 
\end{align}

で記述される。$u$ は移流速度、$\nu$ は拡散係数である。

　時刻 $t=T$ で位置 $x=0$ にトレーサーが存在するとして、 $t=0$において、トレーサーがどこに存在するか推定する。この推定は、勾配 $g(x)$ を求めることで得られる。初期値 $C(0,x)$ を勾配法で更新する際、$g(x)$ が大きくなる $x$ が、初期値 $C(0,x)$ を大幅に更新する。つまり、$g(x)$ が大きくなる $x$ が、$t=T$ で位置 $x=0$にいるトレーサーに最も影響を与え、その $x$ が推定したい値である。

　ラグランジュ関数と評価関数を以下のように定義する。

\begin{align}
\mathcal{L} &= J + \int_0^T \mathrm{d}t \int \mathrm{d}x  \lambda\left( -u \frac{\partial C}{\partial x} +\nu \frac{\partial^2 C}{\partial x^2} -\frac{\partial C}{\partial t} \right)
\\
&J= C(T,x=0)= \int \mathrm{d}x C(T,x)\delta(x)
\end{align}

変分すると、

\begin{align}
\delta\mathcal{L} &= \delta J + \int_0^T \mathrm{d}t \int \mathrm{d}x  \delta\lambda\left( -u \frac{\partial C}{\partial x} +\nu \frac{\partial^2 C}{\partial x^2} -\frac{\partial C}{\partial t} \right)\\
&+ \int_0^T \mathrm{d}t \int \mathrm{d}x \lambda\left( -u \frac{\partial }{\partial x}\delta C +\nu \frac{\partial^2 }{\partial x^2}\delta C -\frac{\partial }{\partial t}\delta C \right)
\end{align}

となる。部分積分を行い、各項を計算すると

\begin{align}
\int_0^T\mathrm{d}t\int\mathrm{d}x\lambda \frac{\partial}{\partial x}\delta C = \int_0^T\mathrm{d}t \left[\lambda \delta C \right]_{-\infty}^{\infty}-\int_0^T\mathrm{d}t\int\mathrm{d}x\frac{\partial \lambda}{\partial x} \delta C
\end{align}

\begin{align}
& \int_0^T \mathrm{d}t \int \mathrm{d}x \lambda\frac{\partial^2 }{\partial x^2}\delta C = \int_0^T\mathrm{d}t \left[\lambda \frac{\partial }{\partial x}\delta C  \right]_{-\infty}^{\infty} -  \int_0^T \mathrm{d}t \int \mathrm{d}x \frac{\partial \lambda}{\partial x}\frac{\partial }{\partial x}\delta C \\
&=\int_0^T\mathrm{d}t \left[\lambda \frac{\partial }{\partial x}\delta C  \right]_{-\infty}^{\infty} - \int_0^T\mathrm{d}t \left[ \frac{\partial \lambda}{\partial x}\delta C  \right]_{-\infty}^{\infty}
+\int_0^T \mathrm{d}t \int \mathrm{d}x \frac{\partial^2 \lambda}{\partial x^2}\delta C \\
\end{align}

\begin{align}
 \int_0^T \mathrm{d}t &\int \mathrm{d}x \lambda\frac{\partial }{\partial t}\delta C =\int \mathrm{d}x \left[\lambda\delta C  \right]_{0}^{T} - \int_0^T \mathrm{d}t \int \mathrm{d}x \frac{\partial \lambda}{\partial t}\delta C \\
&=\int \mathrm{d}x \left[\lambda(T,x)\delta C(T,x) -\lambda(0,x)\delta C(0,x)  \right] - \int_0^T \mathrm{d}t \int \mathrm{d}x \frac{\partial \lambda}{\partial t}\delta C \\
\end{align}

\begin{align}
&\delta J = \int \mathrm{d}x \delta(x) \delta C(T,x)
\end{align}

となる。無限遠では物理量はゼロとすると、$\delta C(\pm \infty)=0, \lambda(\pm \infty)=0$ である。また、

\begin{align}
\delta\mathcal{L}-\delta J(C(0,x))&=\delta\mathcal{L}-\left\langle g,\delta C(0,x)\right\rangle　\\
&=\delta\mathcal{L}-\int\mathrm{d}x g \delta C(0,x)\\
&=0\\
\end{align}

であり、任意の$\delta C$ で成立するので、アジョイント方程式

\begin{align}
&\frac{\partial \lambda}{\partial t} + u\frac{\partial \lambda}{\partial x}+\nu\frac{\partial^2 \lambda}{\partial x^2}=0 \\
&\lambda(T,x)= \delta(x) \\ 
&\lambda(0,x) = g
\end{align}

が導かれる。

次に、具体的にアジョイント方程式を解く。$\lambda(t,x)$ を

\begin{align}
\lambda(t,x) = \int \mathrm{d}k \hat{\lambda}(k) e^{i\left(kx-\omega t  \right)}
\end{align}

とする。アジョイント方程式に代入し、$\omega$ を求めると、

\begin{align}
 &\int \mathrm{d}k \hat{\lambda}(k)e^{i\left(kx-\omega t  \right)
}\left\{i\omega - iuk+k^2\nu \right\}=0\\
& \therefore \omega = uk+i \nu k^2
\end{align}

となる。また、上式から$\lambda(T,x)$を求めると、デルタ関数の定義から、

\begin{align}
\lambda(T,x)&= \delta(x) = \frac{1}{2\pi} \int \mathrm{d}k e^{ikx}\\
&=\int \mathrm{d}k \hat{\lambda}(k) e^{i\left(kx-\omega T \right)}\\
\end{align}

　$\hat{\lambda}(k)$ が

\begin{align}
\therefore \ \ &\hat{\lambda}(k) e^{-i\omega T} = \frac{1}{2\pi}\\
& \hat{\lambda}(k) = \frac{1}{2\pi} \exp\left\{iukT-\nu k^2T \right\}
\end{align}

求まる。最終的に勾配は、複素ガウス積分を使い

\begin{align}
g(x) = \lambda(0,x)&=\frac{1}{2\pi}\int \mathrm{d}k \exp\left\{ik(x+uT)-\nu k^2T \right\}\\
&=\frac{1}{2\pi}\int \mathrm{d}k \exp\left\{-\nu T\left(k-\frac{i(x+uT)}{2\nu T} \right)^2 - \frac{(x+uT)^2}{4\nu T} \right\}\\
&=\frac{1}{2\sqrt{\pi \nu T}}\exp\left\{- \frac{(x+uT)^2}{4\nu T} \right\}\\
\end{align}

平均 $-u T$ 、分散$2\nu T$ のガウス分布が得られる。
　よって、時刻 $t=0$ におけるトレーサーは、$x =-u T$ 付近に存在していたと推定できる。

最後に

　数値予報などでよく使われる、4次元変分法について紹介した。天気予報の他にも、地震波のシミュレーションとデータ同化を組み合わせることで。地下のプレートの構造を推定するなどの研究が行われている。
　
今後は、偏微分方程式を用いて、データ同化の実装を行いたい。

参考文献

記事は、以下の書籍を参考にしています。

データ同化　観測・実験とモデルを融合するイノベーション
淡路敏之・蒲地政文・池田元美・石川洋一編著