cms.ode.diffeq2nd_euler のソースコード

import sys
import csv
import numpy as np
from math import exp, sqrt, sin, cos, pi


"""
  Solve second order diffrential equation by Euler method
"""


#===================
# parameters
#===================
outfile = 'diffeq2nd_euler.csv'
x0 = 0.0
v0 = 1.0
dt = 0.01
nt = 501
iprint_interval = 20


if __name__ == "__main__":
    if (len(sys.argv)) >= 2:
        dt = float(sys.argv[1])
    if (len(sys.argv)) >= 3:
        nt = int(sys.argv[2])
    if (len(sys.argv)) >= 4:
        iprint_interval = int(sys.argv[3])


# d2x/dt2 = force(x,t)
# define function to be integrated


def force(t, x):
    return -x


# solution: x(0.0) = 0.0 v(0.0) = 1.0


def xsolution(t):
    return sin(t)



#===================
# main routine
#===================


def main(x0, v0, dt, nt):
    print("Solve second order diffrential equation by Euler 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)'))
    for i in range(nt):
        t0 = i * dt
# dx/dt = v
# dv/dt = force(t, x)
        v1 = v0 + dt * force(t0, x0)
        x1 = x0 + dt * v0
        xexact = xsolution(t0)
        if i == 1 or i % iprint_interval == 0:
            print("t={:5.2f}  {:12.6f}  {:12.6f}  {:12.6f}".format(t0, x0, xexact, v0))
        fout.writerow([t0, x0, xexact, v0])
        x0 = x1
        v0 = v1

    f.close()


if __name__ == '__main__':
    main(x0, v0, dt, nt)