import numpy as np
from math import exp, sqrt, sin, cos, pi
"""
概要: 1階常微分方程式を4段4次のルンゲ・クッタ法で解くモジュール。
詳細説明:
時間発展する系のシミュレーションに用いられるルンゲ・クッタ法の実装例を提供する。
具体的には、初期値問題 dx/dt = f(t, x) に対して数値解を求め、厳密解と比較する。
このモジュールは、オイラー法、ホイン法、そして4段4次のルンゲ・クッタ法の実装を含んでいる。
関連リンク:
:doc:`diffeq_rungekutta4_usage`
"""
#===================
# parameters
#===================
x0 = 1.0
dt = 0.1
nt = 10000001
iprint_interval = 1000000
# dx/dt = dxdt(x,t)
# define function to be integrated
[ドキュメント]
def dxdt(t, x):
"""
概要: 1階常微分方程式の右辺 `dx/dt = f(t, x)` を定義する。
詳細説明:
この関数は、積分されるべき微分方程式の右辺 `f(t, x)` を表す。
この例では `dx/dt = -x*x` を実装している。
:param t: float: 時刻。
: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):
"""
概要: 定義された微分方程式の厳密解を返す。
詳細説明:
`dx/dt = -x*x` の厳密解は `x = 1 / (C + t)` で与えられる。
初期条件 `x(0)=1.0` を満たす場合、定数 `C` は `1.0` となるため、
この関数は `x = 1 / (1.0 + t)` を返す。
:param t: float: 時刻。
:returns: float: 時刻 `t` における厳密解の値。
"""
return 1.0 / (1.0 + t)
[ドキュメント]
def diffeq_euler(diff1func, t0, x0, dt):
"""
概要: オイラー法を用いて1ステップの数値積分を実行する。
詳細説明:
微分方程式 `dx/dt = f(t, x)` に対し、次の時刻 `t0 + dt` における `x` の値を、
`x(t0+dt) = x0 + dt * f(t0, x0)` の式を用いて計算する。
これは最も基本的な数値積分法であり、精度は低い。
:param diff1func: function: `dx/dt = f(t, x)` を定義する関数。
引数として `(t, x)` を取り、`float` を返す。
:param t0: float: 現在の時刻 `t`。
:param x0: float: 現在の従属変数 `x` の値。
:param dt: float: 時間ステップ幅。
:returns: float: `t0 + dt` における `x` の予測値 `x1`。
"""
k1 = dt * diff1func(t0, x0)
x1 = x0 + k1
return x1
[ドキュメント]
def diffeq_heun(diff1func, t0, x0, dt):
"""
概要: ホイン法（改良オイラー法）を用いて1ステップの数値積分を実行する。
詳細説明:
微分方程式 `dx/dt = f(t, x)` に対し、予測子・修正子のアプローチを用いて
次の時刻 `t0 + dt` における `x` の値を計算する。
具体的には、まずオイラー法で予測値 `x0+k0` を計算し、それを用いて勾配 `k1` を修正。
最終的な `x(t0+dt)` は `x0 + (k0 + k1) / 2.0` で求められる。
オイラー法よりも精度が高い。
:param diff1func: function: `dx/dt = f(t, x)` を定義する関数。
引数として `(t, x)` を取り、`float` を返す。
:param t0: float: 現在の時刻 `t`。
:param x0: float: 現在の従属変数 `x` の値。
:param dt: float: 時間ステップ幅。
:returns: float: `t0 + dt` における `x` の予測値 `x1`。
"""
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):
"""
概要: 1階常微分方程式を4段4次のルンゲ・クッタ法で数値積分し、結果を出力するメインルーチン。
詳細説明:
初期値問題 dx/dt = f(t, x) に対して、指定された時間ステップ数 `nt` だけ
4段4次のルンゲ・クッタ法を適用し、数値解と厳密解を比較しながらコンソールに出力する。
ルンゲ・クッタ法は高い精度を持つ数値積分法である。
計算の最初のステップ (`x1` の計算) にはホイン法を使用している。
指定された `iprint_interval` ごとに結果を出力する。
:param x0: float: 初期時刻 `t=0` における従属変数 `x` の初期値。
:param dt: float: 時間ステップ幅。
:param nt: int: シミュレーションの総ステップ数。
:returns: None: 結果を標準出力に表示するため、明示的な戻り値はない。
"""
print("Solve first order diffrential equation by four-stage fourth-order Runge-kutta method")
print("{:^10} {:^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={:10.2f} {:16.10e} {:16.10e}".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+dt, x0+k1)
k3 = dt * dxdt(t0+2.0*dt, x0+2.0*k2)
x2 = x0 + 1.0 / 3.0 * (k0 + 2.0 * k1 + 2.0 * k2 + k3)
xexact = fsolution(t1)
if i % iprint_interval == 0:
print("t={:10.2f} {:16.10e} {:16.10e}".format(t1, x1, xexact))
t0 = t1
x0 = x1
x1 = x2
if __name__ == '__main__':
main(x0, dt, nt)