import csv
import numpy as np
from math import exp, sqrt, sin, cos, pi
"""
Verlet法による2階微分方程式の数値解法モジュール。
概要:
このモジュールは、Verlet法（修正Euler法を含む初期ステップ）を使用して
与えられた2階常微分方程式を数値的に解きます。
結果はコンソールとCSVファイルに出力されます。
詳細説明:
d^2x/dt^2 = force(x, t) の形式の2階微分方程式を扱います。
初期条件 x(0)=x0, v(0)=v0 から時間発展を計算します。
出力ファイルは 'diffeq2nd_verlet.csv' です。
関連リンク:
:doc:`diffeq2nd_verlet_usage`
"""
# d2x/dt2 = force(x,t)
# define function to be integrated
[ドキュメント]
def force(t, x):
"""
2階微分方程式の右辺である力を計算します。
概要:
与えられた時間 `t` と位置 `x` における力を返します。
この例では、単振動に対応する復元力 -x を計算します。
:param t: float -- 現在の時間。
:param x: float -- 現在の位置。
:returns: float -- 計算された力 (d^2x/dt^2)。
"""
return -x
# solution: x(0.0) = 0.0 v(0.0) = 1.0
[ドキュメント]
def xsolution(t):
"""
2階微分方程式の厳密解を計算します。
概要:
初期条件 x(0)=0.0, v(0)=1.0 に対応する厳密解を返します。
この例では、単振動の厳密解である sin(t) を返します。
:param t: float -- 現在の時間。
:returns: float -- 厳密解における位置 x。
"""
return sin(t)
#===================
# parameters
#===================
outfile = 'diffeq2nd_verlet.csv'
x0 = 0.0
v0 = 1.0
dt = 0.01
nt = 501
iprint_interval = 20
#===================
# main routine
#===================
[ドキュメント]
def main(x0, v0, dt, nt):
"""
Verlet法を用いて2階微分方程式を数値的に解き、結果をファイルに出力します。
概要:
初期条件と時間ステップ、総ステップ数に基づいて、
2階微分方程式の数値解をVerlet法で計算し、結果をCSVファイルとコンソールに出力します。
詳細説明:
最初の1ステップは修正Euler法を用いて計算し、
その後はVerlet法（中心差分）を用いて位置と速度を更新します。
計算された時刻、数値解の位置、厳密解の位置、数値解の速度が出力されます。
:param x0: float -- 初期位置。
:param v0: float -- 初期速度。
:param dt: float -- 時間刻み幅。
:param nt: int -- 計算の総ステップ数。
:returns: None
"""
print("Solve second order diffrential equation by Verlet method")
print("Write to [{}]".format(outfile))
# open outfile to write a csv file
f = open(outfile, 'w')
fout = csv.writer(f, lineterminator='\n')
fout.writerow([
't', 'x(cal)', 'x(exact)', 'v(cal)'
])
print("{:^5} {:^12} {:^12} {:^12}".format('t', 'x(cal)', 'x(exact)', 'v(cal)'))
kv0 = dt * force(0.0, x0)
kv1 = dt * force(0.0+dt, x0+kv0)
v1 = v0 + (kv0 + kv1) / 2.0
kx0 = dt * v0
kx1 = dt * v1
x1 = x0 + (kx0 + kx1) / 2.0
xexact = xsolution(0.0)
print("t={:5.2f} {:12.6f} {:12.6f} {:12.6f}".format(0.0, x0, xexact, v0))
fout.writerow([0.0, x0, xexact, v0])
for i in range(1, nt):
t1 = i * dt
f1 = force(t1, x1)
x2 = 2.0 * x1 - x0 + dt * dt * f1
v1 = (x2 - x0) / 2.0 / dt
xexact = xsolution(t1)
if i == 1 or i % iprint_interval == 0:
print("t={:5.2f} {:12.6f} {:12.6f} {:12.6f}".format(t1, x1, xexact, v1))
fout.writerow([t1, x1, xexact, v1])
x0 = x1
v0 = v1
x1 = x2
f.close()
if __name__ == '__main__':
main(x0, v0, dt, nt)