import numpy as np
from math import exp, sqrt, sin, cos, pi
"""
シンプソン法を用いて1階常微分方程式を解くモジュール。
このモジュールは、指定された1階常微分方程式をオイラー法、ホイン法、
シンプソン法（Runge-Kutta法にシンプソン則の重みを適用したもの）を用いて
数値的に解くための関数群を提供します。また、特定の微分方程式に対する
解析解も含まれており、数値解との比較が可能です。
:doc:`diffeq_simpson_usage`
"""
# dx/dt = dxdt(x,t)
# define function to be integrated
[ドキュメント]
def dxdt(t, x):
"""評価対象となる微分方程式 dx/dt = -x^2 の右辺を計算します。
詳細説明:
時刻 `t` と変数 `x` の値を受け取り、微分方程式の右辺である `-x*x` を計算します。
:param t: 時刻 `t`。
:type t: float
:param x: 変数 `x` の値。
:type x: float
:returns: 微分方程式の右辺 `dx/dt` の値。
:rtype: float
"""
return -x*x
# solution: x = 1 / (C + t), C = 1 for x(0) = 1.0
[ドキュメント]
def fsolution(t):
"""微分方程式 dx/dt = -x^2 の解析解を計算します。
詳細説明:
初期条件 x(0)=1.0 の場合の解析解 `x = 1 / (1 + t)` を返します。
:param t: 時刻 `t`。
:type t: float
:returns: 時刻 `t` における解析解 `x` の値。
:rtype: float
"""
return 1.0 / (1.0 + t)
[ドキュメント]
def diffeq_euler(diff1func, t0, x0, dt):
"""オイラー法を用いて1階常微分方程式を1ステップ進めます。
詳細説明:
`x(t+dt) = x(t) + dt * f(t, x(t))` の式に基づき、
次の時刻における `x` の値を計算します。
:param diff1func: 評価対象となる微分方程式の右辺を計算する関数 (`f(t, x)` の形式)。
:type diff1func: callable
:param t0: 現在の時刻 `t`。
:type t0: float
:param x0: 現在の変数 `x` の値。
:type x0: float
:param dt: 時間刻み幅。
:type dt: float
:returns: 次の時刻 `t + dt` における `x` の値。
:rtype: float
"""
k1 = dt * diff1func(t0, x0)
x1 = x0 + k1
return x1
[ドキュメント]
def diffeq_heun(diff1func, t0, x0, dt):
"""ホイン法（修正オイラー法）を用いて1階常微分方程式を1ステップ進めます。
詳細説明:
オイラー法で求めた予測値と、その予測値を用いて再計算した傾きの平均を取ることで、
より精度の高い次の時刻の `x` の値を計算します。
:param diff1func: 評価対象となる微分方程式の右辺を計算する関数 (`f(t, x)` の形式)。
:type diff1func: callable
:param t0: 現在の時刻 `t`。
:type t0: float
:param x0: 現在の変数 `x` の値。
:type x0: float
:param dt: 時間刻み幅。
:type dt: float
:returns: 次の時刻 `t + dt` における `x` の値。
:rtype: float
"""
k0 = dt * diff1func(t0, x0)
k1 = dt * diff1func(t0+dt, x0+k0)
x1 = x0 + (k0 + k1) / 2.0
return x1
[ドキュメント]
def diffeq_simpson(diff1func, t0, x0, dt):
"""シンプソン法（に似た重みを持つRunge-Kutta法）を用いて1階常微分方程式を1ステップ進めます。
詳細説明:
3つの傾き `k0, k1, k2` を計算し、シンプソン則の重み `(1/3, 4/3, 1/3)` に似た
加重平均 `(k0 + 4*k1 + k2) / 3` を用いて、次の時刻における `x` の値を計算します。
この実装はRunge-Kutta法のバリアントとして、シンプソン則の重み付けを適用しています。
:param diff1func: 評価対象となる微分方程式の右辺を計算する関数 (`f(t, x)` の形式)。
:type diff1func: callable
:param t0: 現在の時刻 `t`。
:type t0: float
:param x0: 現在の変数 `x` の値。
:type x0: float
:param dt: 時間刻み幅。
:type dt: float
:returns: 次の時刻 `t + dt` における `x` の値。
:rtype: float
"""
k0 = dt * dxdt(t0, x0 )
k1 = dt * dxdt(t0+dt, x0+k0)
k2 = dt * dxdt(t0+2.0*dt, x0+k0+k1)
x1 = x0 + 1.0 / 3.0 * (k0 + 4.0 * k1 + k2)
return x1
#===================
# parameters
#===================
x0 = 1.0
dt = 0.1
nt = 501
iprint_interval = 20
#===================
# main routine
#===================
[ドキュメント]
def main(x0, dt, nt):
"""1階常微分方程式をシンプソン法で数値的に解き、解析解と比較して結果を表示します。
詳細説明:
シミュレーションパラメータに基づいて、時刻 `t`、数値計算結果 `x(cal)`、
および解析解 `x(exact)` をコンソールに出力します。
初期値と時間刻み幅を用いて、ループ内でシンプソン法により `x` を更新していきます。
`iprint_interval` で指定された間隔で結果を出力します。
:param x0: 初期時刻 `t=0` における `x` の初期値。
:type x0: float
:param dt: 時間刻み幅。
:type dt: float
:param nt: 総ステップ数。
:type nt: int
:returns: None
:rtype: None
"""
print("Solve first order diffrential equation by Simpson method")
print("{:^5} {:^12} {:^12}".format('t', 'x(cal)', 'x(exact)'))
t0 = 0.0
f0 = dxdt(t0, x0)
x1 = diffeq_euler(dxdt, t0, x0, dt)
# x1 = diffeq_heun(dxdt, t0, x0, dt)
xexact = fsolution(t0)
print("t={:5.2f} {:12.6f} {:12.6f}".format(t0, x0, xexact))
for i in range(1, nt):
t1 = i * dt
xexact = fsolution(t1)
if i % iprint_interval == 0:
print("t={:5.2f} {:12.6f} {:12.6f}".format(t1, x1, xexact))
x2 = diffeq_simpson(dxdt, t0, x0, dt)
t0 = t1
x0 = x1
x1 = x2
if __name__ == '__main__':
main(x0, dt, nt)