import numpy as np
from math import exp, sqrt, sin, cos, pi
"""
三段階三次のルンゲクッタ法を用いた一次常微分方程式の数値解法モジュール。
概要:
本モジュールは、初期値問題を持つ一次常微分方程式 dx/dt = f(t, x) を
三段階三次のルンゲクッタ法を用いて数値的に解く機能を提供します。
また、厳密解との比較表示、およびオイラー法とホイン法の実装も含まれています。
詳細説明:
与えられた初期条件 x(t0) = x0 と時間ステップ dt を用いて、指定された総ステップ数 nt まで
常微分方程式を数値的に積分します。
計算された数値解は、既知の厳密解と比較され、定期的な間隔 (iprint_interval) で
時間 t、数値解 x(cal)、厳密解 x(exact) の形式で標準出力されます。
ルンゲクッタ法の初期ステップの予測には、ホイン法またはオイラー法が利用されます。
関連リンク:
:doc:`diffeq_rungekutta_usage`
"""
#===================
# parameters
#===================
x0 = 1.0
dt = 0.1
nt = 501
iprint_interval = 20
# dx/dt = dxdt(x,t)
# define function to be integrated
[ドキュメント]
def dxdt(t, x):
"""
積分対象となる常微分方程式 `dx/dt = f(t, x)` の右辺を定義します。
概要:
一次常微分方程式 `dx/dt = f(t, x)` の右辺 `f(t, x)` を計算します。
詳細説明:
この例では、`f(t, x) = -x*x` という具体的な関数を実装しています。
:param t: float: 時間 `t` の値。
:param x: float: 変数 `x` の値。
:returns: float: `t`と`x`における `dx/dt` の値。
"""
return -x*x
# solution: x = 1 / (C + t), C = 1 for x(0) = 1.0
[ドキュメント]
def fsolution(t):
"""
与えられた時間 `t` における常微分方程式の厳密解を計算します。
概要:
定義された常微分方程式の厳密解 `x` の値を時間 `t` に対して計算します。
詳細説明:
本例で対象とする `dx/dt = -x*x` (初期条件 `x(0)=1.0`) の厳密解は
`x = 1.0 / (1.0 + t)` です。この関数はこの厳密解を返します。
:param t: float: 時間 `t` の値。
:returns: float: 時間 `t` における厳密解 `x` の値。
"""
return 1.0 / (1.0 + t)
[ドキュメント]
def diffeq_euler(diff1func, t0, x0, dt):
"""
オイラー法を用いて次の時間ステップの `x` の値を予測します。
概要:
オイラー法 (Euler method) を適用して、与えられた初期値から次の時間ステップ `t0 + dt`
における `x` の値を計算します。
詳細説明:
`x_new = x_old + dt * f(t_old, x_old)` の単純な一次近似によって、
次のステップの `x` の値を予測します。精度は低いですが、最も基本的な数値積分法です。
:param diff1func: function: `dx/dt = f(t, x)` の右辺を計算する関数。
この関数は `(t, x)` を引数にとり、`f(t, x)` を返します。
:param t0: float: 現在の時間 `t` の値。
:param x0: float: 現在の変数 `x` の値。
:param dt: float: 時間ステップ幅。
:returns: float: 時間 `t0 + dt` における `x` の予測値。
"""
k1 = dt * diff1func(t0, x0)
x1 = x0 + k1
return x1
[ドキュメント]
def diffeq_heun(diff1func, t0, x0, dt):
"""
ホイン法（修正オイラー法）を用いて次の時間ステップの `x` の値を予測します。
概要:
ホイン法 (Heun's method) または修正オイラー法を適用して、
与えられた初期値から次の時間ステップ `t0 + dt` における `x` の値を計算します。
詳細説明:
オイラー法による予測値 (`k0`) を計算した後、その予測値を用いてより正確な
傾き (`k1`) を計算し、それらの平均を用いて `x` を更新します。
`x_new = x_old + (k0 + k1) / 2.0`。これにより、オイラー法よりも高い精度が得られます。
:param diff1func: function: `dx/dt = f(t, x)` の右辺を計算する関数。
この関数は `(t, x)` を引数にとり、`f(t, x)` を返します。
:param t0: float: 現在の時間 `t` の値。
:param x0: float: 現在の変数 `x` の値。
:param dt: float: 時間ステップ幅。
:returns: float: 時間 `t0 + dt` における `x` の予測値。
"""
k0 = dt * diff1func(t0, x0)
k1 = dt * diff1func(t0+dt, x0+k0)
x1 = x0 + (k0 + k1) / 2.0
return x1
#===================
# main routine
#===================
[ドキュメント]
def main(x0, dt, nt):
"""
三段階三次のルンゲクッタ法で常微分方程式を解き、結果を出力するメインルーチン。
概要:
与えられた初期値、時間ステップ幅、総ステップ数に基づき、
三段階三次のルンゲクッタ法を用いて常微分方程式の数値解を計算し、結果を表示します。
詳細説明:
この関数は、まず初期値を設定し、最初のステップの `x` 値をホイン法（またはオイラー法）で予測します。
その後、指定されたステップ数 `nt` にわたって三段階三次のルンゲクッタ法を繰り返し適用し、
`x` の値を更新していきます。
`iprint_interval` で指定された間隔で、現在の時間 `t`、計算された `x` の値、
および厳密解 `x_exact` をコンソールに出力し、数値解の精度を評価できるようにします。
:param x0: float: 初期時刻 `t=0` における `x` の初期値。
:param dt: float: 時間ステップ幅。各積分ステップで時間が `dt` だけ進みます。
:param nt: int: シミュレーションの総ステップ数。`nt` 回の積分が実行されます。
:returns: None: 計算結果を標準出力に表示します。
"""
print("Solve first order diffrential equation by three-stage third-order Runge-kutta method")
print("{:^5} {:^16} {:^16}".format('t', 'x(cal)', 'x(exact)'))
t0 = 0.0
f0 = dxdt(t0, x0)
# The next x (x1) must be predicted by Euler or Heum method
# x1 = diffeq_euler(dxdt, t0, x0, dt)
x1 = diffeq_heun(dxdt, t0, x0, dt)
xexact = fsolution(t0)
print("t={:5.2f} {:16.10f} {:16.10f}".format(t0, x0, xexact))
for i in range(1, nt):
t1 = i * dt
k0 = dt * dxdt(t0, x0)
k1 = dt * dxdt(t0+dt, x0+k0)
k2 = dt * dxdt(t0+2.0*dt, x0+4.0*k1-2.0*k0)
x2 = x0 + 1.0 / 3.0 * (k0 + 4.0 * k1 + k2)
xexact = fsolution(t1)
if i % iprint_interval == 0:
print("t={:5.2f} {:16.10f} {:16.10f}".format(t1, x1, xexact))
t0 = t1
x0 = x1
x1 = x2
if __name__ == '__main__':
main(x0, dt, nt)