常微分方程式の数値解法(1) - 流体解析のあれこれ

SciPyのodeintを用いて以下のOrdinary Differential Equations (ODE)を解き，文献[1]の数値解を再現する．

$\displaystyle y_1^{\prime} = y_2$

$\displaystyle y_2^{\prime} = \eta \left( 1 - y_1^2 \right) y_2 - y_1$

ここで， $\eta = 100.0$ ，初期値 $y_1\left( 0 \right) = 2.0$ ， $y_2 \left( 0 \right) = 0.0$ である．

[1] Petzold, L., Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations, SIAM, J. Sci. Stat. Comput., 5(1), 136-148 (1983)

#
# Van der Pol's equation in
# Petzold, L.,
# Automatic selection of methods for solving stiff and nonstiff systems
# of ordinary differential equations,
# SIAM J. Sci. Stat. Comput., 4(1), 136-148 (1983)
#
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# constant
eta = 100.0

# right-hand side
def func(y, x, eta):
    return [y[1],
            eta*(1 - y[0]**2)*y[1] - y[0]]

# initial condition
y0 = [2.0, 0.0]

# output interval
x = np.arange(0.0, 1000.0, 1.0)

# solve ode
y = odeint(func, y0, x, args=(eta,))

# output
plt.xlabel('x')
plt.ylabel('y')
plt.plot(x, y[:,0])
#plt.plot(x, y[:,1])
#plt.show()
plt.savefig('Petzold1983-Fig1.png')