scikit-FiniteDiff を使ってみる

# Scikit-FiniteDiff を試用する

差分法の偏微分方程式ソルバー [Scikit-FiniteDiff](https://scikit-fdiff.readthedocs.io/en/latest/introduction.html) をインストールして使い始めています。[^1]
[^1]: https://joss.theoj.org/papers/10.21105/joss.01356

リンク先のイントロダクションに概要が書かれています。
`sympy` を使って離散化した方程式を作り、時間に対して Method of line の手法で解くとのこと。まずは、道具として使う観点で触っています。

## 環境構築

仮想環境を作って pip  インストールして使い始めることができる。(Windows11, Python 3.10.11)

```pwsh
> cd %Userprofile% 
> .\AppData\Local\Programs\Python\Python310\python.exe -m venv --upgrade-deps skfdiff
> .\skfdiff\Scripts\Activate.ps1
> pip install scikit-fdiff[interactive,numba]
> ipython kernel install --user --name=py310_skfdiff
```

## チュートリアル例題 ['One dimension Korteweg–de Vries'](https://scikit-fdiff.readthedocs.io/en/latest/auto_examples/plot_1D_kdv.html) を試す

```python
%matplotlib inline
```


```python
import numpy as np
# import pylab as pl
import matplotlib.pyplot as plt
from skfdiff import Model, Simulation
```

`Zarr module not found, ZarrContainer not available.` という Warning が出るけれど、先にすすむことはできる。

### モデルを組み立てて解く

例題は、水深の浅い水面を進む波（浅水波）を記述する方程式として知られる Korteweg–de Vries (KdV) 方程式ということだが、
検索して上位で見つかる例では $a=0$ に相当する式が使われていることが多い。
ソリトン解で有名。

$
\frac{\partial U}{\partial t} + U\frac{\partial U}{\partial x} = a\frac{\partial^2 U}{\partial x^2} + b\\frac{\partial^3 U}{\partial x^3}
$

これを変形して、
$
\frac{\partial U}{\partial t} = -U\frac{\partial U}{\partial x} + a\frac{\partial^2 U}{\partial x^2} + b\frac{\partial^3 U}{\partial x^3}
$

右辺を `Model` にする



```python
model = Model("-U * dxU + a * dxxU + b * dxxxU", 
              "U(x)", 
              ["a", "b"]
             )
```

最初の引数が微分方程式で、第 2 の引数が未知数と独立変数を示す。第 3 の引数はパラメータである。関数の引数名としては、それぞれ 'evolution_equations', 'unknowns', 'parameters' となっている。


次に変数の範囲を定め、初期値を作る


```python
x = np.linspace(-2, 6, 1000)

n = 20
U = np.log(1 + np.cosh(n) ** 2 / np.cosh(n * x) ** 2) / (2 * n)

initial_fields = model.fields_template(x=x, U=U, a=2e-4, b=1e-4)
```

初期値の形を見てみよう


```python
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.plot(x, U)
ax.set_xlabel("x")
ax.set_ylabel("U")
ax.set_title("Initial value")
```


![output_7_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/01726c7b-f606-3b9b-bdbb-3a9653cf04c1.png)

これを使って、「シミュレーション」`Simulation` を組み立て、数値解析を実効することができる。
このとき経過を蓄積する `container` を指定できる。


```python
simulation = Simulation(model, initial_fields, dt=0.05, tmax=10)
container = simulation.attach_container()
```


```python
simulation.run()
```

    (10.0,
     <xarray.Dataset>
     Dimensions:  (x: 1000)
     Coordinates:
       * x        (x) float64 -2.0 -1.992 -1.984 -1.976 ... 5.976 5.984 5.992 6.0
     Data variables:
         U        (x) float64 -5.501e-08 -5.347e-08 -4.966e-08 ... -0.04111 -0.04193
         a        float64 0.0002
         b        float64 0.0001)

### 結果の確認

`simulation.run()` は、最後のステップのデータだけを返すが、`container` は指定に応じて履歴データを保持できる。 
`container` オブジェクトの使い方を見ていこう。
ちなみに `simulation.attach_container()` の返り値と `simulation.container` は同じものである。container.data に、`x`, `t` と解 U が保存されている。


```python
container == simulation.container  ## True
```

```python
container.data
```

![output_10_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/9f457fa7-7584-f237-2b71-40898108beae.png)


結果のプロットを見る


```python
container.data.U[:: container.data.t.size // 8].plot(col="t", col_wrap=3, color="black")
```

![output_16_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/f7734aab-54af-d3c7-baa8-3d4d81acf905.png)


右側の短い波長で振動している部分がソリトン解ということのようだ。
ちなみに .plot の機能は、skfdiff 固有の機能ではなくて、 [xarray](https://docs.xarray.dev/en/latest/user-guide/plotting.html) の機能である。


2D データを全部プロットさせることもできる


```python
container.data.U.plot()
```


![output_25_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/105d9ba5-b65c-4739-b7ff-5d5d84e0c8f7.png)


波形の重ね書きもできる


```python
container.data.U.isel(t=slice(0, None, 50)).plot.line(hue="t")
```

![output_37_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/cc4957d0-118f-2020-c378-c981f22e1e5d.png)




```python
container.data.U.coords
```


xarray の container.data.U は、以下のような構造

![output_24_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/38979/8d6e8d7a-394d-d018-221e-d05f9c0e03d6.png)


座標に直接アクセスすることができる。

```python
container.data.U.coords
```


    Coordinates:
      * x        (x) float64 -2.0 -1.992 -1.984 -1.976 ... 5.976 5.984 5.992 6.0
      * t        (t) float64 0.0 0.05 0.1 0.15 0.2 0.25 ... 9.8 9.85 9.9 9.95 10.0


### まとめと今後

* 解のハンドリングを自分で考えるより、フレームワークを活用すると低コストでソルバーを使えそう。
* パラメータを変えて、ソリトン解の文献との比較。
* 別の微分方程式を入れて、境界条件の指定や、`hook`, 風上差分指定などの機能を見ていきたい。
* 覚書: 差分法の他のツールとして、見つけたもの
  * [py-pde](https://py-pde.readthedocs.io/en/latest/)
  * [Devito](https://www.devitoproject.org/devito_book/)
  * 一般的に言えば、任意の領域で解ける FEM の方が汎用性は高い。skfdiff の[イントロダクションページ](https://scikit-fdiff.readthedocs.io/en/latest/introduction.html#rational-and-software-history)で論じられている

## 付録
全部一式の [Jupyter notebook](https://gist.github.com/natsuwater/d90e25115b7abc21cbcb7817888f98e7) も作りました