ODE sample by using odeint package in scipy.

odeint_scipy.py

# -*- coding: utf-8 -*-
 
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint


def main():

    x_ini = 0.0 # x(0) = 0
    v_ini = 0.0 # v(0) = d/dt * x(0) = 0
    g = 9.8

    def eom(X,t):
        '''
        X=[x,v]
        dx/dt = F(t,x,v) = v
        dv/dt = G(t,x,v) = -g
        '''
        return [X[1],g]

    h = 0.1 # ステップ幅(s)
    T = 1000 # ステップ回数

    # 理論値(ここは比較のため)
    X_ideal =[x_ini]
    for i in range(T):
        tt = h*(i+1)
        X_ideal.append(-(g*tt**2)/2)

    # ここから
    t = np.arange(0,T*h+h/10,h) # 時間幅
    c = odeint(eom,np.array([x_ini,v_ini]),t)

    # 描画
    a = c[:,0]
    b = np.array(X_ideal)
    plt.plot(t,b,color='green',label='X_correct')
    plt.plot(t,a,color='black',label='X_calculate')
    plt.legend(bbox_to_anchor=(0.01, 0.01), loc='lower left', borderaxespad=0)
    plt.xlabel("Time [s]")
    plt.show()

if __name__ == '__main__':
    main()