import csv
import numpy as np
from math import exp, sqrt, sin, cos, pi
"""
Solve simulataneous second order diffrential equations by Verlet method
概要: 2次元の2階常微分方程式をVerlet法で解くモジュール。
詳細説明: 2次元の単振動方程式 d2x/dt2 = -x, d2y/dt2 = -y をVerlet法で数値積分し、解析解と比較します。
結果はCSVファイルに出力されます。
関連リンク: :doc:`diffeq2nd_2d_verlet_usage`
"""
# d2x/dt2 = force(t, x, y), d2y/dt2 = force(t, x, y)
# define function to be integrated
[ドキュメント]
def force(t, x, y):
"""
概要: 2階常微分方程式の右辺（加速度項）を定義します。
詳細説明: 現在の例では d2x/dt2 = -x, d2y/dt2 = -y に対応します。
:param t: float: 時間
:param x: float: x座標
:param y: float: y座標
:returns: tuple[float, float]: x方向の力（加速度）、y方向の力（加速度）
"""
return -x, -y
# solution: x(0.0) = 0.0, y(0.0) = 2.0, vx(0.0) = 1.0, vy(0.0) = 0.0
[ドキュメント]
def xysolution(t):
"""
概要: 2階常微分方程式の解析解を提供します。
詳細説明: 初期条件 x(0)=0, y(0)=2, vx(0)=1, vy(0)=0 に対応する解析解 x(t)=sin(t), y(t)=2cos(t) を返します。
:param t: float: 時間
:returns: tuple[float, float]: 時刻tにおけるx座標の解析解、y座標の解析解
"""
return sin(t), 2.0 * cos(t)
#===================
# parameters
#===================
outfile = 'diffeq2nd_2D_Verlet.csv'
x0, vx0 = 0.0, 1.0
y0, vy0 = 2.0, 0.0
dt = 0.01
nt = 501
iprint_interval = 20
#===================
# main routine
#===================
[ドキュメント]
def main(x0, y0, vx0, vy0, dt, nt):
"""
概要: 2次元の2階常微分方程式をVerlet法で数値積分し、結果をファイルに保存します。
詳細説明: 指定された初期条件と時間ステップ、総ステップ数に基づき、Verlet法を用いて運動をシミュレーションします。
計算結果は標準出力とCSVファイルに出力され、解析解との比較も行われます。
初期ステップはHeun法で計算されます。
:param x0: float: 時刻 t=0 における初期x座標
:param y0: float: 時刻 t=0 における初期y座標
:param vx0: float: 時刻 t=0 における初期x方向速度
:param vy0: float: 時刻 t=0 における初期y方向速度
:param dt: float: 時間ステップ幅
:param nt: int: 総時間ステップ数
:returns: None
"""
print("Solve simulataneous second order diffrential equations 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)', 'y(cal)', 'y(exact)'
])
print("{:^5} {:^12} {:^12} {:^12} {:^12}".format('t', 'x(cal)', 'x(exact)', 'y(cal)', 'y(exact)'))
# Solve the 1st data by Heun method
kvx0, kvy0 = force(0.0, x0, y0)
kvx1, kvy1 = force(0.0+dt, x0+kvy0, y0+kvy0)
vx1 = vx0 + dt * (kvx0 + kvx1) / 2.0
vy1 = vy0 + dt * (kvy0 + kvy1) / 2.0
kx0 = dt * vx0
kx1 = dt * vx1
ky0 = dt * vy0
ky1 = dt * vy1
x1 = x0 + (kx0 + kx1) / 2.0
y1 = y0 + (ky0 + ky1) / 2.0
xexact, yexact = xysolution(0.0)
print("t={:5.2f} {:12.6f} {:12.6f} {:12.6f} {:12.6f}".format(0.0, x0, xexact, y0, yexact))
fout.writerow([0.0, x0, xexact, y0, yexact])
for i in range(1, nt):
t1 = i * dt
fx1, fy1 = force(t1, x1, y1)
x2 = 2.0 * x1 - x0 + dt * dt * fx1
y2 = 2.0 * y1 - y0 + dt * dt * fy1
vx1 = (x2 - x0) / 2.0 / dt
vy1 = (y2 - y0) / 2.0 / dt
xexact, yexact = xysolution(t1)
if i == 1 or i % iprint_interval == 0:
print("t={:5.2f} {:12.6f} {:12.6f} {:12.6f} {:12.6f}".format(t1, x1, xexact, y1, yexact))
fout.writerow([t1, x1, xexact, y1, yexact])
x0 = x1
vx0 = vx1
y0 = y1
vy0 = vy1
x1 = x2
y1 = y2
f.close()
if __name__ == '__main__':
main(x0, y0, vx0, vy0, dt, nt)