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

1. 概要

　この記事では、線形常微分方程式の概要を説明し、オブジェクト指向に基づいて対応するクラスの実装を行います。シリーズ全体で、オブジェクト指向プログラミングを用いた線形常微分方程式を数値計算するためのパッケージを提供する予定です。

2. 線形常微分方程式

　変数 $x$ の未知関数 $f(x)$ の $k$ 階導関数を $f^{(k)}(x)\quad (k\in \mathbb{N})$ と表します。$N$ を自然数として、 $r(x)$ および $c_k(x)\quad (k=0,\ldots,N-1)$ を $x$ についての既知関数とすると、$N$ 階線形常微分方程式の一般形は次式で与えられます。
$$ f^{(N)}(x) + c_{N-1}(x) f^{(N-1)}(x) + \ldots + c_1(x) f^{(1)}(x) + c_0(x) f(x) = r(x). \tag{1}$$
これ以降では、関数の引数を省略します。
　$N$ 階線形常微分方程式は、$N$ 個の未知関数の一階微分のみを含む連立線形常微分方程式に書き直すことができます。$\phi_k=f^{(k)}$ $(k=0,\ldots,N-1)$ と置くと、式$(1)$は次の方程式系で表されます。

\left\{
\begin{align}
&\phi_0^\prime = \phi_1, \\
&\phi_1^\prime = \phi_2, \\
&\qquad\vdots \\
&\phi_{N-1}^\prime = r - (c_0\phi_0 + \ldots +c_{N-1}\phi_{N-1}).\\
\end{align}
\right.
\tag{2}

ここで、上付きのプライムは $x$ による微分を表しています。解析的・数値的を問わず、解を求める際は一階連立微分方程式の形に変形しておくと取り扱いやすい場合が多くあります。このシリーズの次回以降で数値解法を実装する際も、連立微分方程式の形で毎ステップの計算を行います。

3. クラスの実装

　以上を踏まえて、$N$ 階線形常微分方程式のクラスを作成します。式$(1)$ または式$(2)$より、斉次項 $r(x)$ と未知関数の係数項 $c_k(x)\quad (k=0,\ldots,N-1)$ の情報を持たせると、常微分方程式の形が決まります。よって、クラスのコンストラクタを次のように定義します。

ordinary_differential_equation.py

class OrdinaryDifferentialEquation:
    """ Class of linear ordinary differential equation.
    """
    def __init__(self, coeff_funcs:list[FunctionType],
            rhs_func:FunctionType):
        """ Constructor of OrdinaryDifferentialEquation.

        Parameters
        ----------
        order : int
            Order of ODE.
        coeff_funcs : list of function
            Coefficient functions of ODE.
            Index k -> coeffcient of k-th order derivative.
        rhs_func : function
            Right hand side function of ODE.
        """
        self.__coeff_funcs = coeff_funcs
        self.__rhs_func = rhs_func
        self.__order = len(coeff_funcs)

係数項 $c_k(x)\quad (k=0,\ldots,N-1)$ は self.__coeff_funcs に、斉次項 $r(x)$ は self.__rhs_func に格納します。両者ともアンダースコア二つ始まりの名前にすることで、private変数にしています。また、数値計算において常微分方程式の次数を取得することが多いため、あらかじめ self.__order に代入しておきます。
　常微分方程式の形を途中で変更する必要はないと思いますので、既知関数はコンストラクタでのみ設定できるようにします（セッター関数は作成しない）。一方、既知関数（の値）と微分方程式の次数は数値計算の過程で必要ですので、取得用のメソッドを作成します。

ordinary_differential_equation.py

class OrdinaryDifferentialEquation:
    # 中略

    def get_coeff(self, x:float=None)\
             -> np.ndarray | list[FunctionType]:
        """ Get coefficient functions in ODE.
        """
        if x is None:
            return self.__coeff_funcs
        else:
            return np.array(list(
                map(lambda f: f(x), self.__coeff_funcs)
            ))

    def get_rhs(self, x:float=None) -> float | FunctionType:
        """ Get right hand side function in ODE.
        """
        if x is None:
            return self.__rhs_func
        else:
            return self.__rhs_func(x)

    def get_order(self,) -> int:
        """ Get the order of ODE.
        """
        return self.__order

メソッドに $x$ の値を渡した場合は対応する関数の値を返し、引数なしでは関数自体を返すようにしています。

4. 使い方と例

　最初に、インスタンスの生成方法を説明します。例として微分方程式 $f^{\prime\prime} - (1/x) f^\prime + (1/x^2)f = x^2$ を取り上げます。

example_ode.py

c0 = lambda x: 1.0/x**2
c1 = lambda x: -1.0/x
r = lambda x: x**2
ode = OrdinaryDifferentialEquation([c0,c1], r)

第一引数には $c_k(x)$ のリストを、第二引数には $r(x)$ を渡します。$c_k(x)$ はリストとインデックスが同一になる順番で与えます。
　get_order() メソッドで微分方程式の次数を取得します。

example_ode.py

print(ode.get_order())
    # >>> 2

　get_coeff() と get_rhs() である $x$ に対する $c_k(x)$, $r(x)$ の値を取得します。

example_ode.py

print(ode.get_coeff(2.0), ode.get_rhs(2.0))
    # >>> [ 0.25 -0.5 ] 4.0

引数を渡さないときは関数オブジェクトが返ってきます。

example_ode.py

coeff = ode.get_coeff()
rhs = ode.get_rhs()
print(list(map(lambda f: f(2.0), coeff)), rhs(2.0))
    # >>> [0.25, -0.5] 4.0

5. 付録

ソースコード全文

ordinary_differential_equation.py

""" Class of linear ordinary differential equation. """

# import libraries
from __future__ import annotations
from types import FunctionType
import numpy as np

# class definition
class OrdinaryDifferentialEquation:
    """ Class of linear ordinary differential equation.
    """
    def __init__(self, coeff_funcs:list[FunctionType],
            rhs_func:FunctionType):
        """ Constructor of OrdinaryDifferentialEquation.

        Parameters
        ----------
        order : int
            Order of ODE.
        coeff_funcs : list of function
            Coefficient functions of ODE.
            Index k -> coeffcient of k-th order derivative.
        rhs_func : function
            Right hand side of ODE.
        """
        self.__coeff_funcs = coeff_funcs
        self.__rhs_func = rhs_func
        self.__order = len(coeff_funcs)

    def get_coeff(self, x:float=None)\
             -> np.ndarray | list[FunctionType]:
        """ Get coefficients in ODE.
        """
        if x is None:
            return self.__coeff_funcs
        else:
            return np.array(list(
                map(lambda f: f(x), self.__coeff_funcs)
            ))

    def get_rhs(self, x:float=None) -> float | FunctionType:
        """ Get right hand side in ODE.
        """
        if x is None:
            return self.__rhs_func
        else:
            return self.__rhs_func(x)

    def get_order(self,) -> int:
        """ Get the order of ODE.
        """
        return self.__order
 

example_ode.py

""" Example of using OrdinaryDifferentialEquation. """
from ordinary_differential_equation import OrdinaryDifferentialEquation

if __name__ == "__main__":
    a0 = lambda x: 1.0/x**2
    a1 = lambda x: -1.0/x
    r = lambda x: x**2
    ode = OrdinaryDifferentialEquation([a0,a1], r)

    print(ode.get_order())
        # >>> 2

    print(ode.get_coeff(2.0), ode.get_rhs(2.0))
        # >>> [ 0.25 -0.5 ] 4.0

    coeff = ode.get_coeff()
    rhs = ode.get_rhs()
    print(list(map(lambda f: f(2.0), coeff)), rhs(2.0))
        # >>> [0.25, -0.5] 4.0

参考文献

• 常微分方程式 - Wikipedia

後続記事はこちら

• オブジェクト指向を用いた線形常微分方程式の数値計算（Python）2-1：差分法基底クラス - Qiita