オブジェクト指向を用いた線形常微分方程式の数値計算（Python）2-2：Euler陽解法クラス

# 1. はじめに
　常微分方程式の離散化方法の一つであるEuler陽解法を解説し、Pythonで数値計算を行うためのEuler陽解法クラスを実装します。

　同シリーズの記事はこちらです。
- [オブジェクト指向を用いた線形常微分方程式の数値計算（Python）1：線形常微分方程式クラス](https://qiita.com/itotomball/items/8c25ffa93070603af794)
- [オブジェクト指向を用いた線形常微分方程式の数値計算（Python）2-1：差分法基底クラス](https://qiita.com/itotomball/items/344789710d0ba1b9ebcd)

# 2. Euler陽解法
　$N$ 階線形常微分方程式は、式$(1)$のような一階連立常微分方程式に書き換えられました。
$$ \phi_k^\prime = \rho_k(x,\phi_0,\ldots,\phi_{N-1}),\qquad (k=0,\ldots,N-1). \tag{1}$$
差分法では、この左辺の一階微分を差分で置き換えます。
$$ \phi^\prime \quad\rightarrow\quad \frac{\phi^{n+1}-\phi^n}{\Delta x}. \tag{2} $$
右辺の関数 $\rho_k$ を評価する方法は数値計算スキームによって異なりますが、今回紹介するEuler陽解法では現在ステップ $x_n$ における値を用います。
```math
\begin{align}
\\
\rho_k(x,\phi_0,\ldots,\phi_{N-1}) \quad&\rightarrow\quad \rho_k^n :=\rho_k(x_n,\phi_0^n,\ldots,\phi_{N-1}^n). \tag{3}
\end{align}
```
式$(2)$, $(3)$の置き換えを行うことで、Euler陽解法の公式が求まります。
$$ \phi_k^{n+1} = \phi_k^n + \rho_k^n \Delta x. \tag{4}$$
$\rho_k^n$ の具体形は以下のとおりでした。
```math
\left\{
\begin{align}
&\rho_0^n = \phi_1^n, \\
&\rho_1^n = \phi_2^n, \\
&\qquad\vdots \\
&\rho_{N-1}^n = r^n - (c_0^n\phi_0^n + \ldots +c_{N-1}^n\phi_{N-1}^n).\\
\end{align}
\right.
\tag{5}
```
これを式$(4)$に代入すると、以下の公式を得ます。
```math
\left\{
\begin{align}
  \phi_k^{n+1} &= \phi_k^n + \phi_{k+1}^n \Delta x,&(&k=0,\ldots,N-2),\\
  \phi_{N-1}^{n+1} &= \phi_{N-1}^n + \left[ r^n - (c_0^n\phi_0^n + \ldots +c_{N-1}^n\phi_{N-1}^n) \right] \Delta x,&(&k=N-1). 
\end{align}
\right.
\tag{6}
```

# 3. 実装
### 3.1 クラスの継承
　前回実装した `FiniteDifferenceMethod` クラスを継承して、`ExplicitEuler` クラスを定義します。ここでは、前回と同じファイルに書いていきます。
``` finite_difference_method.py
""" Class of Finite Difference Method for ODE.
"""
class FiniteDifferenceMethod(metaclass=ABCMeta):
    # 実装部省略

class ExplicitEuler(FiniteDifferenceMethod):
    """ Class of explicit Euler method.
    """
```
ここで、`ExplicitEuler` クラスのコンストラクタは定義しません。`FiniteDifferenceMethod` クラスのコンストラクタをそのまま使用します。

### 3.2 ステップの進行
 前回、`FiniteDifferenceMethod` クラスの仮想関数として定義した `advance_step()` を、継承先のクラス `ExplicitEuler` でオーバーライドします。式$(6)$をもとに、次のように実装します。
``` finite_difference_method.py
    # 中略

class ExplicitEuler(FiniteDifferenceMethod):
    """ Class of explicit Euler method.
    """
    def _compute_explicit_rhs(self, phi:np.ndarray,
            x:float) -> np.ndarray:
        """ Compute right hand side of ODE explicitly.

        Parmaters
        ---------
        phi : array-like, shape=(order,)
            Function value at step n.
        x : float
            Current value of x.

        Return
        ------
        np.ndarray, shape=(order,), dtype=float
            Right hand side values of ODE system.
        """
        order = self.get_order()
        rhs = np.zeros(order)

        # phi_{N-1}^\prime = r^n - sum_{k=0}^{N-1} c_{k}^n phi_{k}^n
        coeff = self.get_coeff(x)
        rhs[order-1] = self.get_rhs(x) - np.dot(coeff, phi)

        # phi_{k}^\prime = phi_{k+1}^n
        rhs[0:order-1] = phi[1:order]

        return rhs


    def advance_step(self, phi:np.ndarray, x:float, x_inc:float) -> None:

        rho = self._compute_explicit_rhs(phi, x)
        phi[:] += x_inc * rho[:]
```
ここで、`_compute_explicit_rhs()` メソッドを作り、右辺の関数の評価値を求める部分を分けています。後に Runge-Kutta 法のクラスを実装する際にメソッドを再利用できるようにするため、このようにしています。

# 4. 付録
### ソースコード
``` finite_difference_method.py
""" Class of Finite Difference Method for ODE.
"""
# import libraries
from __future__ import annotations
from abc import ABCMeta, abstractmethod
from types import FunctionType
import numpy as np
from . import OrdinaryDifferentialEquation

# Base class of finite difference scheme
class FiniteDifferenceMethod(metaclass=ABCMeta):
    """ Base class of finite difference methods for ODE.
    """
        # 省略

# Explicit Euler scheme
class ExplicitEuler(FiniteDifferenceMethod):
    """ Class of explicit Euler method.
    """
    def _compute_explicit_rhs(self, phi:np.ndarray,
            x:float) -> np.ndarray:
        """ Compute right hand side of ODE explicitly.

        Parmaters
        ---------
        phi : array-like, shape=(order,)
            Function value at step n.
        x : float
            Current value of x.

        Return
        ------
        np.ndarray, shape=(order,), dtype=float
            Right hand side values of ODE system.
        """
        order = self.get_order()
        rhs = np.zeros(order)

        # phi_{N-1}^\prime = r^n - sum_{k=0}^{N-1} c_{k}^n phi_{k}^n
        coeff = self.get_coeff(x)
        rhs[order-1] = self.get_rhs(x) - np.dot(coeff, phi)

        # phi_{k}^\prime = phi_{k+1}^n
        rhs[0:order-1] = phi[1:order]

        return rhs


    def advance_step(self, phi:np.ndarray, x:float, x_inc:float) -> None:

        rho = self._compute_explicit_rhs(phi, x)
        phi[:] += x_inc * rho[:]

```

# 5. 参考文献
- [差分法 - Wikipedia](https://ja.wikipedia.org/wiki/%E5%B7%AE%E5%88%86%E6%B3%95)
- [オブジェクト指向を用いた線形常微分方程式の数値計算（Python）1：線形常微分方程式クラス - Qiita](https://qiita.com/itotomball/items/8c25ffa93070603af794)
- [オブジェクト指向を用いた線形常微分方程式の数値計算（Python）2-1：差分法基底クラス - Qiita](https://qiita.com/itotomball/items/344789710d0ba1b9ebcd)