4.1. 2階の常微分方程式(scipy.integrate.odeint)

物理の問題として、二階の常微分方程式

\[\frac{d^2 x}{dt^2} + \gamma \frac{dx}{dt} + \omega^2 x = \epsilon f(t)\]

は頻出の問題である.odeintを用いて解く際には一階の微分方程式

\[\begin{split}\dot{x} &= p\\ \dot{p} &= -\gamma \dot{x} - \omega^2 x + \epsilon f(t)\end{split}\]

と変形して取り扱う方がよいであろう.

本例題では外力として

\[f(t) = \cos(\omega_2 t)\]

の強制振動を考える.

まず,微分方程式

# funcs.py

import numpy as np
twopi=2.0*np.pi
def NonAuto(x,t=0,gamma=1.0,omega=1.0,eps=1.0,omega2=np.sqrt(2.0)):
    """
    d2x0/dt2 + gamma*dx0/dt + (omega**2)*x0 = f(t)
    ->
    dx0/dt=x1
    dx1/dt=-gamma*x1 - (omega**2)*x0
    """
    def ft(t):
        return eps*np.cos(omega2*t)

    x0 = x[1]
    x1 = -gamma*x[1] -(omega**2)*x[0] - ft(t)
    return np.array([x0,x1])

を定義(ここでは funcs.py として保存)し, その後共通のルーチン(order2diff.py)として,

# order2diff.py

import numpy as np
import pylab
import scipy.integrate
import funcs

name = 'non-auto'

f=funcs.NonAuto 
""" parameter setting """
gamma, omega = 0.1, np.sqrt(3.0)
eps, omega2= 1.0, np.sqrt(3.0)
args = (gamma, omega, eps, omega2)

""" initial condition """
sample=5
x0 = np.linspace(0.0,1.0,sample)
x1 = np.array([0.0 for i in range(sample)])
initsets = [np.array(x) for x in zip(x0,x1)]

""" integration """
dt=0.01
t = np.arange(0, 30, dt)
trajectories = [scipy.integrate.odeint(f,x0, t, args, Dfun=None) for x0 in initsets ]

vcolors = pylab.cm.autumn_r(np.linspace(0.3, 1.0, sample)) # color parameters
for i,x in enumerate(trajectories):
    x = x.T
    pylab.plot(t,x[0],'-',lw=2) #,label='(%.f, %.f)' % (x0[i],x1[i]))
pylab.xlabel("time",fontsize=20)
pylab.ylabel("x",fontsize=20,rotation='horizontal')
#pylab.savefig(name+'.eps')
pylab.show()

を実行すると次のような結果が得られる.

fig1 fig2 fig3
左から順に \(\epsilon=0\) 減衰振動, \(\epsilon=1\) 強制振動, \(\epsilon=1, \omega=\omega_2\) 共鳴振動