1. 概要
この記事では、常微分方程式の差分法について説明し、対応する基底クラスの実装を行います。前回の記事はこちらです。
2. 差分法
2.1 線形常微分方程式
前回の内容をおさらいします。以下の線形常微分方程式を対象としています。
$$ f^{(N)} + c_{N-1} f^{(N-1)} + \ldots + c_1 f^{(1)} + c_0 f = r. \tag{1}$$
$\phi_k=f^{(k)}$ $(k=0,\ldots,N-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}
次節以降で表記を簡単にするため、式$(2)$の右辺を $\rho_k(x,\phi_0,\ldots,\phi_{N-1})$ と置いて、次のように表すことにします。
$$ \phi_k^\prime = \rho_k(x,\phi_0,\ldots,\phi_{N-1})\quad (k=0,\ldots,N-1) \tag{3}$$
2.2 差分近似
独立変数 $x$ の領域の離散化を行い、 $n$ ステップ目の値を $x_n$ と表します。多くの場合ではステップ幅は一定にとられるため、$x_n=x_0+n\Delta x$ と表せますが、単調に増加する数列になっていれば一定幅でなくとも構いません。各点に対応する関数の値を $\phi^n_k=\phi_k(x_n)$ とします。次ステップの関数値を求めるため、差分法では一階微分を次のように差分で置き換えます。
$$ \phi^\prime \rightarrow \frac{\phi^{n+1}-\phi^n}{\Delta x} \tag{4}$$
この式を用いると、式$(3)$から差分方程式を得ます。
$$ \frac{\phi_k^{n+1}-\phi_k^n}{\Delta x} = \rho_k \tag{5}$$
この方程式を解くと次ステップの関数値 $\phi_k^{n+1}$ が求まります。ですが、ここで右辺の関数をどのステップの値で評価するかが問題であり、Euler法やCrank-Nicolson法など様々な数値計算スキームが考案されています。各手法の具体的な内容については個別の記事で解説します。この記事では、差分法の数値解法全般に共通する部分を抽出して基底クラスを実装します。
3. 実装
3.1 ライブラリ・モジュールのインポート
このクラスでは抽象関数を作成するので、抽象基底クラスを実装するための標準ライブラリである abc から ABCMeta と abstractmethod をimportします（abc は標準ライブラリなのでインストールする必要はありません）。また、前回実装した ordinary_differential_equation.py から OrdinaryDifferentialEquation クラスをimportしておきます。 annotations, FunctionType は関数のアノテーションに使うためです。
from __future__ import annotations
from abc import ABCMeta, abstractmethod
from types import FunctionType
import numpy as np
from ordinary_differential_equation import OrdinaryDifferentialEquation
3.2 コンストラクタ
いずれの差分スキームにおいても、解きたい問題の微分方程式の形を定義する OrdinaryDifferentialEquation オブジェクトが必要です。よって、基底クラス FiniteDifferenceMethod のコンストラクタでは OrdinaryDifferentialEquation オブジェクトを引数に取り、privateなインスタンス変数 self.__ode に代入します。
class FiniteDifferenceMethod:
""" Base class of finite difference methods for ODE.
"""
def __init__(self, ode:OrdinaryDifferentialEquation):
""" Constructor.
Parameters
----------
ode : OrdinaryDifferentialEquation object
ODE to solve using FDM.
"""
self.__ode = ode
3.3 ステップの進行（抽象関数）
差分法では差分方程式を解き、次ステップの関数の値を求めるメソッドが共通して必要になります。しかし、現段階ではどの差分スキームを使用するか決まっていません。そこで、抽象メソッドととして定義だけをしておき、派生クラスで各スキームの実装を行うようにします。
抽象メソッドにするには、@abstractmethod デコレータをつけます。@abstractmethod をつけることで、派生先のクラスで advance_step メソッドが未実装でないかチェックしてくれます。具体的な実装は一切書きませんが、メソッド内に文がないとエラーを出すので、docstringを書いておくと良いです（passでも大丈夫です）。引数はそれぞれ、phi:現ステップにおける未知関数の値、x:現ステップの独立変数の値、x_inc:独立変数の増分値です。
class FiniteDifferenceMethod:
# 中略
@abstractmethod
def advance_step(self, phi:np.ndarray, x:float, x_inc:float) -> None:
""" Advance a step to update state of function.
Parameters
----------
phi : array-like
Function values at current step.
x : float
Current value of variable.
x_inc : float
Increment of variable.
"""
3.4 getter関数
派生クラスの差分スキームの計算で使用するために、インスタンス変数 __ode のgetter関数を呼ぶメソッドを定義しておきます（self.__ode.get_order() などと何回も書くのは面倒なので）。
class FiniteDifferenceMethod:
# 中略
def get_order(self,) -> int:
""" Get the order of ODE.
"""
return self.__ode.get_order()
def get_coeff(self, x:float=None)\
-> np.ndarray | list[FunctionType]:
""" Get coefficient functions in ODE.
"""
return self.__ode.get_coeff(x)
def get_rhs(self, x:float=None) -> float | FunctionType:
""" Get right hand side function in ODE.
"""
return self.__ode.get_rhs(x)
4. 付録
ソースコード全文
""" Classes of Finite Difference Method for ODE.
"""
from __future__ import annotations
from abc import ABCMeta, abstractmethod
from types import FunctionType
import numpy as np
from ordinary_differential_equation import OrdinaryDifferentialEquation
class FiniteDifferenceMethod(metaclass=ABCMeta):
""" Base class of finite difference methods for ODE.
"""
def __init__(self, ode:OrdinaryDifferentialEquation):
""" Constructor.
Parameters
----------
ode : OrdinaryDifferentialEquation object
ODE to solve using FDM.
"""
self.__ode = ode
@abstractmethod
def advance_step(self, phi:np.ndarray, x:float, x_inc:float) -> None:
""" Advance a step to update state of function.
Parameters
----------
phi : array-like
Function values at current step.
x : float
Current value of variable.
x_inc : float
Increment of variable.
"""
def get_order(self,) -> int:
""" Get the order of ODE.
"""
return self.__ode.get_order()
def get_coeff(self, x:float=None)\
-> np.ndarray | list[FunctionType]:
""" Get coefficient functions in ODE.
"""
return self.__ode.get_coeff(x)
def get_rhs(self, x:float=None) -> float | FunctionType:
""" Get right hand side function in ODE.
"""
return self.__ode.get_rhs(x)