import sys
import csv
import numpy as np
from math import exp, sqrt, sin, cos, pi
import matplotlib.pyplot as plt
"""オイラー法による1階常微分方程式の数値解法モジュール。
概要:
本モジュールは、1階常微分方程式 dx/dt = force(t, x) をオイラー法で数値的に解き、
その結果を厳密解と比較してCSVファイルに出力し、リアルタイムでグラフ表示します。
初期条件、時間刻み幅、総ステップ数などはコマンドライン引数で指定可能です。
詳細説明:
メイン関数 `main` は、指定された初期値 `x0`、時間刻み幅 `dt`、総ステップ数 `nt` を用いて、
`diffeq_euler` 関数を繰り返し呼び出し、各時刻での `x` の値を計算します。
同時に、`fsolution` 関数で厳密解を計算し、両者を比較します。
結果はコンソールに定期的に出力され、`diffeq_euler.csv` ファイルに保存されます。
また、matplotlib を使用して計算過程がリアルタイムで可視化されます。
関連リンク:
:doc:`diffeq_euler_usage`
"""
#===================
# parameters
#===================
outfile = 'diffeq_euler.csv'
x0 = 1.0
dt = 0.5
nt = 501
iprint_interval = 20
if __name__ == "__main__":
argv = sys.argv
n = len(argv)
if n >= 2:
x0 = float(argv[1])
if n >= 3:
dt = float(argv[2])
if n >= 4:
nt = int(argv[3])
if n >= 5:
iprint_interval = int(argv[4])
# dx/dt = force(x,t)
# define function to be integrated
[ドキュメント]
def force(t, x):
"""積分対象となる微分方程式の右辺 dx/dt を定義する。
概要:
この関数は、積分対象となる微分方程式 dx/dt の右辺の値を計算して返します。
詳細説明:
現在の実装では dx/dt = -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):
"""定義された微分方程式の厳密解を返す。
概要:
本関数は、force関数で定義された微分方程式の厳密解を計算して返します。
詳細説明:
force関数で `dx/dt = -x*x` と定義され、初期条件 `x(0)=1.0` が与えられた場合、
この微分方程式の厳密解は `x(t) = 1 / (1 + t)` となります。
この関数はその厳密解の値を計算します。
:param t: float: 時間 t。
:returns: float: 時間 t における厳密解 x(t) の値。
"""
return 1.0 / (1.0 + t)
[ドキュメント]
def diffeq_euler(force, t0, x0, dt):
"""1ステップのオイラー法を適用して次の時刻のxの値を計算する。
概要:
オイラー法を用いて、現在の状態から微小時間 dt 後の変数の値を近似的に計算します。
詳細説明:
現在の時刻 `t0` と変数 `x0` の値、そして時間刻み幅 `dt` を使用して、
`x1 = x0 + dt * force(t0, x0)` の式に基づき、次の時刻 `t0 + dt` における
変数 `x` の推定値 `x1` を計算します。
:param force: callable: 微分方程式の右辺 `dx/dt` を計算する関数。
この関数は `force(t, x)` の形式である必要があります。
:param t0: float: 現在の時刻。
:param x0: float: 現在の変数 x の値。
:param dt: float: 時間刻み幅。
:returns: float: 次の時刻 (`t0 + dt`) における x の推定値 `x1`。
"""
k1 = dt * force(t0, x0)
x1 = x0 + k1
return x1
#===================
# main routine
#===================
[ドキュメント]
def main(x0, dt, nt):
"""メインルーチンとして、オイラー法による微分方程式の数値計算を実行し、結果を出力・表示する。
概要:
指定された初期条件、時間刻み幅、総ステップ数に基づき、オイラー法で微分方程式を解き、
結果をCSVファイルに出力し、リアルタイムでグラフ表示します。
詳細説明:
本関数は、グローバル変数またはコマンドライン引数で設定された初期値 `x0`、
時間刻み幅 `dt`、総ステップ数 `nt` を用いて、微分方程式の数値積分を実行します。
計算結果は、`outfile` で指定されたCSVファイルに時間、計算値、厳密解の形式で書き込まれます。
また、matplotlib を用いて、計算値と厳密解の比較グラフが動的に更新・表示されます。
`iprint_interval` に基づき、計算途中経過がコンソールにも出力されます。
最後にユーザーがEnterキーを押すまでグラフを表示し続けます。
:param x0: float: 微分方程式の初期条件 x(0) の値。
:param dt: float: 時間刻み幅。
:param nt: int: 計算する総ステップ数。
:returns: None: 明示的な戻り値はなく、プログラムは実行後に終了します。
"""
print("Solve first order diffrential equation by Euler method")
print("Write to [{}]".format(outfile))
# prepare for graph
xt = [0.0]
yxex = [x0]
yxsim = [x0]
fig = plt.figure()
ax1 = fig.add_subplot(1, 1, 1)
# 凡例を表示させるため、とりあえずplot()を呼び出す
# 後でプロット毎にデータリストを再設定するので、lineオブジェクトを受け取っておく
line1, = ax1.plot(xt, yxex, label = 'exact')
line2, = ax1.plot(xt, yxsim, label = 'euler')
# ax1.set_xscale('log')
# ax1.set_yscale('log')
ax1.set_xlabel("t")
ax1.set_ylabel("x(t)")
ax1.legend()
# open outfile to write a csv file
f = open(outfile, 'w')
fout = csv.writer(f, lineterminator='\n')
fout.writerow([
't', 'x(cal)', 'x(exact)'
])
print("{:^5} {:^12} {:^12}".format('t', 'x(cal)', 'x(exact)'))
for i in range(1, nt):
t0 = i * dt
x0 = diffeq_euler(force, t0, x0, dt)
xexact = fsolution(t0)
xt.append(t0)
yxex.append(xexact)
yxsim.append(x0)
# graphをupdateするには、プロットデータ line1/line2 に set_data() でデータリストを設定し、plt.pause()を呼び出す
# set_data() ではグラフの表示範囲は更新されないので、データの最小・最大値で設定する
line1.set_data(xt, yxex)
line2.set_data(xt, yxsim)
ax1.set_xlim((min(xt), max(xt)))
ax1.set_ylim((min(yxsim + yxex), max(yxsim + yxex)))
plt.tight_layout()
plt.pause(0.00001)
if i == 1 or i % iprint_interval == 0:
print("t={:5.2f} {:12.6f} {:12.6f}".format(t0, x0, xexact))
fout.writerow([t0, x0, xexact])
f.close()
print("Press ENTER to exit>>", end = '')
input()
exit()
if __name__ == '__main__':
main(x0, dt, nt)