1. はじめに
常微分方程式の離散化方法の一つであるEuler陽解法を解説し、Pythonで数値計算を行うためのEuler陽解法クラスを実装します。
同シリーズの記事はこちらです。
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$ における値を用います。
\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$ の具体形は以下のとおりでした。
\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)$に代入すると、以下の公式を得ます。
\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 クラスを定義します。ここでは、前回と同じファイルに書いていきます。
""" 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)$をもとに、次のように実装します。
# 中略
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. 付録
ソースコード
""" 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[:]