この記事に関して
ここでは一階線形常微分方程式の数値解を求める方法について、量子コンピュータ理論を用いて説明します。
常微分方程式
一般的に常微分方程式は以下の形で表されます。
F(t, x(t), x'(t), \cdots, x^{(n)}(t)) = 0 \\
\left(x^{(i)}(t):=\frac{\mathrm{d}^i}{\mathrm{d}t^i}x(t)\right)
例えば
x' = x + t\ , \ \left(x^{(3)}\right)^{2} - x' + t = 0
などが挙げられます。
この内
qt + \sum_{i=0}^{n}p_ix^{(i)}(t) = r
の形で表されるものを線形、それ以外を非線形と呼びます。
数値解法
今回は上で説明したものの中で一階の常微分方程式を考えます。
一階常微分方程式
\frac{\mathrm{d}x}{\mathrm{d}t}(t) = f(x)\ , \ x(t_0) = x_0\\
f: \mathbb{R} \rightarrow \mathbb{R}
今回は $t = t_n$ のときの $x(t_n) = x_n$ を求めていきます。
また $f$ は次の条件を考えます。
ヘルダー条件
ヘルダー空間 $F^{r, \rho}$ を以下で定義します。
F^{r, \rho}
:= \left\{f: \mathbb{R} \rightarrow \mathbb{R}\ |\ f \in C^r(\mathbb{R}) \ ,\
\left|f^{(i)}(x)\right|\leq D_i, i \leq r \ ,\\
\left|f^{(r)}(x)- f^{(r)}(y)\right|\leq H \left|x-y\right|^{\rho}, x, y \in \mathbb{R}, 0 < \rho \leq 1
\right\}
ここで各 $D_i, H$ は正の定数です。
$f$ は $r$ 階微分可能で各導関数は有界になります。
微分方程式の解の存在と一意性はリプシッツ条件というものを満たさないと成立しません。
ただし上の空間を考えると $r+\rho\geq 1$ のときリプシッツ条件を満たします。
ですので今回 $f$ は $F^{r, \rho}$ 上の関数とします。
続き
微分方程式を求めやすく変形します。
積分を用いて以下のように $x(t_n)$ を表すことができます。
x(t_n) = x(t_0) + \int_{t_0}^{t_n}f(x(t))dt\ , \
(x(t_0) = x_0)
後半の積分を離散和で近似するためにまずこの積分区間を細かくします。
$t \in [t_0,t_n]$ として区間を $n$ 等分します。
[t_0,t_n] = \sum_{i=0}^{n-1}[t_i, t_{i+1}]\ ,\ \left(t_i = t_0 + \frac{t_n-t_0}{n}i\right)
これをさらに $m$ 等分します。
[t_i, t_{i+1}] = \sum_{j=0}^{m-1}\left[t_j^i, t_{j+1}^i\right]\ ,\ \left(t_j^i = t_i + \frac{t_{i+1}-t_i}{m}j\right)
よって以下のように書き直せます。
x(t_n) = x(t_0) + \sum_{i=0}^{n-1}\sum_{j=0}^{m-1}\int_{t_j^i}^{t_{j+1}^i}f(x(t))dt\ , \
(x(t_0) = x_0)
特に $x(t_i), x(t_{i+1})$ に着目すると
x(t_{i+1}) = x(t_i) + \sum_{j=0}^{m-1}\int_{t_j^i}^{t_{j+1}^i}f(x(t))dt\ , \
(x(t_0) = x_0)
となります。以下はこの式を元に考えていきます。
次に上の積分を求めるために $f(x), x(t)$ をそれぞれ近似することを考えます。
$x(t)$ の $t = t_j^i$ 周りのテイラー展開を $r$ 次まで考えると、
x(t) \approx \sum_{p=0}^{r}\frac{1}{p!}x^{(p)}(t_j^i)(t-t_j^i)^{p} =: l_j^i(t)
また $f(x)$ についても $l_j^i(t_j^i)$ 周りでテイラー展開を $r$ 次まで考えると
f(x) \approx \sum_{q=0}^{r}\frac{1}{q!}f^{(q)}(l_j^i(t_j^i))(x-l_j^i(t_j^i))^{q} =: w_j^i(x)
よって元の式は以下のように近似できます。
x(t_{i+1}) = x(t_i) + \sum_{j=0}^{m-1}\int_{t_j^i}^{t_{j+1}^i}w_j^i(l_j^i(t))dt
+ \sum_{j=0}^{m-1}\int_{t_j^i}^{t_{j+1}^i}\left\{f(x(t)) - w_j^i(l_j^i(t))\right\}dt \\
\approx x(t_i) + \sum_{j=0}^{m-1}\int_{t_j^i}^{t_{j+1}^i}w_j^i(l_j^i(t))dt
$t = t_j^i + u\overline{h}$ とおくと
x(t_{i+1}) \approx x(t_i) + \sum_{j=0}^{m-1}\overline{h}\int_{0}^{1}w_j^i(l_j^i(t_j^i + u\overline{h}))du\ , \
\left(\overline{h} = \frac{t_{i+1}-t_i}{m}\right)
区間を $[0,1]$ にすることができました。
この積分をさらに $v$ 分割します。
x(t_{i+1}) \approx x(t_i) + m\frac{\overline{h}}{mv}\sum_{j=0}^{m-1}\sum_{k=0}^{v-1}w\left(l\left(t_j^i + \frac{k}{v}\overline{h}\right)\right)
これで離散和で書き直すことができました。
ここで $\mathcal{J}_i$ を以下で定義しておきます。
\mathcal{J}_i := \frac{\overline{h}}{mv}\sum_{j=0}^{m-1}\sum_{k=0}^{v-1}w\left(l\left(t_j^i + \frac{k}{v}\overline{h}\right)\right)
スケーリング
次の量子数値積分を用いるために $f$ が $[0,1]$ に値をとるように式変形します。
g\left(l\left(t_j^i + \frac{k}{v}\overline{h}\right)\right) = \frac{f\left(l\left(t_j^i + \frac{k}{v}\overline{h}\right)\right) - f_{\rm{min}}}{\Delta f}\\
f_{\rm{min}} = \min_{0 \leq j \leq m , 0 \leq k \leq v} f\left(l\left(t_j^i + \frac{k}{v}\overline{h}\right)\right)\ ,\
\Delta f = f_{\rm{max}} - f_{\rm{min}}
$f_{\rm{min}} \leq f \leq f_{\rm{max}}$ より、これで $g$ を $[0,1]$ 区間に制限することができました。
よって $\mathcal{J}_i$ をスケーリングした $\mathcal{G}_i$ は以下のようになります。
\mathcal{J}_i \approx \overline{h}f_{\rm min} + \overline{h}\mathcal{G}_i \\
\mathcal{G}_i := \frac{1}{mv}\sum_{j=0}^{m-1}\sum_{k=0}^{v-1} g\left(l\left(t_j^i + \frac{k}{v}\overline{h}\right)\right)\ , \
(0 \leq \mathcal{G}_i \leq 1)
結局 $x(t_{i+1})$ は以下のように近似できます。
x(t_{i+1}) \approx x(t_i) + m\overline{h}f_{\rm min} + m\overline{h}\mathcal{G}_i
最後に $\mathcal{G}_i$ を量子数値積分を使って求めます。
量子数値積分
詳しい内容はここで述べています。
以下の $S(f)$ を量子振幅推定を用いて求めるアルゴリズムです。
S(f) := \sum_{x\in J^d}^{}p(x)f(x)\ , \ J = \{0,1, \cdots, N-1\}
今回は $\mathcal{G}_i$ を量子数値積分で求めます。
$m \rightarrow 2^m, v \rightarrow 2^v$ と置き直します。
\mathcal{G}_i = \sum_{x=0}^{2^m-1}\sum_{y=0}^{2^v-1}\frac{1}{2^{m+v}}g\left(l\left(t_i + x\overline{h}+y\frac{\overline{h}}{2^v}\right)\right) \xrightarrow{\rm Rewrite} \sum_{x=0}^{2^m-1}\sum_{y=0}^{2^v-1}\frac{1}{2^{m+v}}g(x,y)
よって以下のように作用素 $\mathcal{A}_i$ を定義できます。
\mathcal{A}_i\left|0\right>_{m+v+1}
= \sum_{x=0}^{2^m-1}\sum_{y=0}^{2^v-1}\sqrt{\frac{1}{2^{m+v}}}\left|x\right>_m\left|y\right>_v
\left(\sqrt{g(x,y)}\left|0\right>+\sqrt{1-g(x,y)}\left|1\right>\right)
この右辺を書き直すと
\mathcal{A}_i \left|0\right>_{m+v+1}
= \sqrt{\mathcal{G}_i}\left|\Psi_0\right>_{m+v+1} + \sqrt{1-\mathcal{G}_i}\left|\Psi_1\right>_{m+v+1}
これを $\sqrt{\mathcal{G}_i} = \sin\theta$ と考えると振幅推定から $\sin^2\theta = \mathcal{G}_i$ が求められます。
各 $t_i$ で同様に上記の手順を行うことで $x(t_n)=x_n$ を求められます。
以上で微分方程式の数値解を求めることができました。
まとめ
今回は一階常微分方程式の数値解を求める方法について説明しました。
参考文献
Frank Gaitan
Finding flows of a Navier–Stokes fluid through quantum computing
https://www.nature.com/articles/s41534-020-00291-0#MOESM1
量子数値積分
https://qiita.com/KeiichiroHiga/items/dae9b5395366f0f4fc2c
Gilles Brassard, Peter Høyer, Michele Mosca, Alain Tapp
Quantum Amplitude Amplification and Estimation
https://arxiv.org/pdf/quant-ph/0005055.pdf
Bolesław Kacewicz
Almost optimal solution of initial-value problems by randomized and quantum algorithms
https://www.sciencedirect.com/science/article/pii/S0885064X06000148