import numpy as np
from math import exp, sqrt, sin, cos, pi
"""
オイラー法（Heun法）を用いて1階常微分方程式を解くモジュール。
このモジュールは、一般的な1階常微分方程式 `dx/dt = f(t, x)` をHeun法（修正オイラー法）で数値的に解く機能を提供します。
特定の例として、`dx/dt = -x*x` という方程式の数値解を計算し、厳密解と比較します。
:doc:`diffeq_heun_usage`
"""
# dx/dt = force(x,t)
# define function to be integrated
[ドキュメント]
def force(t, x):
"""
積分対象となる関数 `dx/dt = f(t, x)` を定義する。
この関数は `f(t, x) = -x*x` を返す。これは例として使用される特定の微分方程式の右辺である。
:param t: float: 現在の時刻。
:param x: float: 現在の状態変数。
:returns: float: `dx/dt` の値。
"""
return -x*x
# solution: x = 1 / (C + t), C = 1 for x(0) = 1.0
[ドキュメント]
def fsolution(t):
"""
定義された微分方程式の厳密解を計算する。
初期条件 `x(0) = 1.0` の場合、厳密解は `x(t) = 1 / (1 + t)` で与えられる。
この関数は、数値計算の結果と比較するための基準値を提供する。
:param t: float: 時刻。
:returns: float: 時刻 `t` における厳密解 `x` の値。
"""
return 1.0 / (1.0 + t)
[ドキュメント]
def diffeq_heun(force, t0, x0, dt):
"""
Heun法（修正オイラー法）を用いて、次の時刻の状態変数を計算する。
Heun法は、オイラー法よりも精度の高い常微分方程式の数値解法である。
現在の点での勾配と、オイラー法で予測された次の点での勾配の平均を用いて、次のステップを計算する。
1. 現在の状態 `(t0, x0)` からオイラー法で次のステップを予測 (`k0`)。
2. 予測された次のステップ `(t0+dt, x0+k0)` を用いて勾配 `k1` を計算。
3. `k0` と `k1` の平均値を使用して、より正確な次の状態 `x1` を計算する。
:param force: callable: `dx/dt = force(t, x)` の形式で定義される微分方程式の右辺関数。
第一引数に時刻 `t`、第二引数に状態変数 `x` をとる。
:param t0: float: 現在の時刻。
:param x0: float: 現在の状態変数。
:param dt: float: 時間ステップ幅。
:returns: float: 次の時刻 `t0 + dt` における状態変数 `x1` の値。
"""
k0 = dt * force(t0, x0)
k1 = dt * force(t0+dt, x0+k0)
x1 = x0 + (k0 + k1) / 2.0
return x1
#===================
# parameters
#===================
x0 = 1.0
dt = 0.1
nt = 501
iprint_interval = 20
#===================
# main routine
#===================
[ドキュメント]
def main(x0, dt, nt):
"""
メインルーチンとして、Heun法による微分方程式の数値計算を実行し、結果を出力する。
指定された初期条件、時間ステップ、計算ステップ数に基づき、`diffeq_heun` 関数を繰り返し呼び出して数値解を求める。
計算された数値解と厳密解を比較し、指定された間隔 (`iprint_interval`) でコンソールに出力する。
:param x0: float: 微分方程式の初期状態変数 `x(0)` の値。
:param dt: float: 時間ステップ幅。
:param nt: int: 総計算ステップ数。
:returns: None
"""
print("Solve first order diffrential equation by Heun method")
print("{:^5} {:^12} {:^12}".format('t', 'x(cal)', 'x(exact)'))
for i in range(nt):
t = i * dt
x1 = diffeq_heun(force, t, x0, dt)
xexact = fsolution(t)
if i == 1 or i % iprint_interval == 0:
print("t={:5.2f} {:12.6f} {:12.6f}".format(t, x0, xexact))
x0 = x1
if __name__ == '__main__':
main(x0, dt, nt)