雨粒の落下: 数値計算・機械学習・量子計算の３手法で運動方程式を解く

## 1. はじめに

機械学習や量子計算の理解を深めるために、それぞれの手法で同じ微分方程式を解いてみました。今回対象にした微分方程式は "雨粒の落下" を記述する運動方程式です。物理現象としての"雨粒の落下"は、その背後にある ”空気抵抗を含む運動方程式" により支配されています。この運動方程式に対して、3つの解法

- 数値計算（ルンゲクッタ法）
- 機械学習
- 量子計算

でアプローチしてみました。

今回実装したソースコードは [こちらのGitリポジトリ](https://github.com/p-circle-kuh/rain-drops-falling) にあります。


## 2. 雨粒の落下の運動方程式

まず雨粒の落下について整理します。地球上で雨粒が落下する際、主に二つの力が働きます。一つは重力、もう一つは空気抵抗です。これらの力のバランスによって、雨粒の落下速度や落下時間が決まることになります。

運動方程式はこの力のバランスを数式化したもので、微分方程式となります。細かいことを考慮すると、雨粒の形状や空気の密度・粘性などさまざまな要素が影響しますが、ここでは単純なモデル化として雨粒は質点・空気抵抗の大きさは係数kで表現されるとして扱うことにします。そうすると、運動方程式は、具体的には以下のように表されます。

```math
m \frac{dv}{dt} + k v - m g = 0
```
ここで t は時間、v は鉛直下向きを正とした雨粒の落下速度、m は雨粒の質量、g は重力加速度です。この方程式を解くことで、雨粒がどのように落下するのかを予測することが可能となります。

以降では、この運動方程式を3つの異なる手法で解いていきます。



## 3. 解法1: 数値計算（ルンゲクッタ法）

数値計算の解法として、比較的一般的な手法であるルンゲクッタ法を用いて、"雨粒の落下"の運動方程式の解を求めてみます。

ルンゲクッタ法は、微分方程式の解を逐次的に計算するための手法で、各ステップでの微小変化を用いて次の状態を計算します。その特徴は、単純なオイラー法に比べて高精度な近似解を得られる点にあります。Wikipedia等にも詳しい説明が載っているので、ここでは手法の詳細は割愛します。

実際に計算してみて得られた解は下図のようになりました。

![rain-drops-falling_Runge-Kutta.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1219280/75e496f6-b46b-0f7b-748c-9ca0a24f940a.png)

## 4. 解法2: 機械学習

機械学習の解法として、Hamiltonian Neural Networks (HNN) を用いました。

HNNは、Neural NetworksでHamiltonianをモデル化し、Hamiltonの運動方程式を満たすように学習を行います。詳細は著者らの[論文](https://proceedings.neurips.cc/paper_files/paper/2019/file/26cd8ecadce0d4efd6cc8a8725cbd1f8-Paper.pdf)をご覧ください。

ルンゲクッタ法の結果に乱数でノイズを加えたものを真値として、学習データ数25、検証データ数25を使用しました。実際に計算してみて得られた解は下図のようになりました。

![rain-drops-falling_HNN.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1219280/9ba45d4f-e63b-57f7-60df-6e11c69fce0c.png)

## 5. 解法3: 量子計算
量子計算の解法として、Xanadu社が開発している量子プログラミング開発環境[Pennylane](https://www.xanadu.ai/products/pennylane)と[Strawberry fields](https://strawberryfields.ai/)を使用した量子機械学習により、問題を解くことを考えます。Strawberry fieldsは光量子コンピュータのフレームワークです。光量子コンピュータによる計算については、[資料１](https://qiita.com/ryuNagai/items/b106f81c17224d512189)[資料2](https://arxiv.org/abs/1804.03159)などをご覧ください。

まず、計算の見通しを良くするための変数変換を行い、$\ell$の微分方程式に書き換えます。今回作成する量子回路は、この$\ell$を推定するための量子回路になります。

```math
\ell = g - kv
```
```math
\frac{d\ell}{dt} =  - ku
```

量子回路は[文献](https://browse.arxiv.org/pdf/2012.12220.pdf)、[コード1](https://github.com/XanaduAI/quantum-neural-networks)[コード2](https://pennylane.ai/qml/demos/quantum_neural_net)、[フォーラムでのディスカッション](https://discuss.pennylane.ai/t/variational-circuits-for-solving-differential-equations/1778)の情報から下記のような構成にしました。シングルqumodeに各種ゲートを作用させていき、最後に測定を行います。
その結果得られるx(position)の期待値が量子回路が出力する推定値になります。この推定値から算出される損失関数をもとに、推定値と正解値の差異が出来るだけ小さくなるよう、古典コンピュータで（量子回路の）最適なパラメータを探索していきます。
古典コンピュータによる計算を必要としているため、この様な手法は量子・古典ハイブリッド機械学習などと呼ばれます。

![量子回路図2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1219280/646bf34e-0c22-f79b-ee4d-a9f526332aab.png)



解析解に微小なノイズを加えた100個のデータを学習データとした場合、$\ell$の出力は下記のようになりました。

![ell_result.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1219280/4390c3f4-1aa8-6e9a-b8ef-26c9849227c6.png)

この結果をもとに、台形公式で速度を計算すると下記のようになり、
ルンゲクッタ法による数値計算の結果より、やや高めの値を出力しています。

![v_result.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1219280/94d0a3a4-8b37-25d4-a008-f5a3823d3117.png)

## 6. まとめと今後の挑戦

- 目標とした3つの方法で、雨粒の自由落下の微分方程式を解くことが出来ました。
- 機械学習・量子計算に関して新たな手法やライブラリを見つけて使うことができ、理解を深めることができました。
- 機械学習と量子計算では真値からのズレがあるので、今後は誤差の定量化と要因分析をして、精度改善につなげていきたいです。