ジャパニーズ・アトラクタを描きたくなったので，SageMathでやってみた．SageMathには一階の常微分方程式系を解くメソッドが用意されている．つまり，次の未知関数$x_1,\ldots,x_n$に関する微分方程式系 \[ \begin{aligned} \frac{dx_1}{dt} & = f_1(t,x_1,\ldots,x_n),\\ \frac{dx_2}{dt} & = f_2(t,x_1,\ldots,x_n),\\ \cdots\\ \frac{dx_n}{dt} & = f_n(t,x_1,\ldots,x_n) \end{aligned} \] を解くことができる．まずは四次のRunge–Kutta法で解いてくれるdesolve_system_rk4を使ってみた．（desolveはDifferential Equation SOLVEだろうか．）使い方は

t = var('t')
x1 = var('x1')
...
xn = var('xn')

と変数を宣言した後，

P = desolve_system_rk4([f_1,...,f_n],[x1,...,xn],ivar=t,ics=[t_0,x_0(t_0),...,x_n(t_0)],end_points=t'_0,step=s)

とする．各引数は

• 第一引数は各微分方程式の右辺．関数を具体的に書く．
• 第二引数は未知関数の変数名．
• ivarは各関数$x_i$の変数．
• icsは初期値からなるリスト．ics[0]は$t$の初期値で，それを$t_0$とするとics[i]は$x_i(t_0)$．
• end_pointsは$t$の最後の値．
• stepは微分を差分に直す時の隣同士の$t$の差．省略すると0.1．

である．戻り値はリストで，$t$の初期値を$t_0$，stepを$s$とすると，$t_n = t_0 + ns$として，P[n][0]には$t_n$が，$P[n][i]$には$x_i(t_n)$が入っている．戻り値から$x_1$のグラフを描くには，例えば$n = 3$の場合には

list_plot([(t1,y1) for t1,y1,y2,y3 in P])

とする．（未知関数として使っているx1,x2,x3とわけるためにy1,y2,y3とした．）list_plotはリストをとるので，Pから適当にリストを生成して渡すことになる，上は各P[n]から(P[n][0],P[n][1])をとって渡している．点同士を線で結ぶならば，

list_plot([(t1,y1) for t1,y1,y2,y3 in P], plotjoined=True) 

とする．

例えば \[ \begin{aligned} \frac{dx}{dt} & = y,\\ \frac{dy}{dt} & = -x \end{aligned} \] を，初期値$t = 0$，$x(0) = 0$，$y(0) = 1$で解いてみる．（答えは言うまでも無く$x = \sin(t)$，$y = \cos(t)$である．

t = var('t')
x = var('x')
y = var('y')
P = desolve_system_rk4([y,-x],[x,y],ivar=t,ics=[0,0,1],end_points=10)
list_plot([(t1,x1) for t1,x1,y1 in P],plotjoined=True)

で$\sin$の関数が出てくる．

list_plot([(x1,y1) for t1,x1,y1 in P],plotjoined=True)

だと円．

問題のジャパニーズ・アトラクタだけど，Wikipediaにあるパラメータだとうまくでてくれないので（何が悪いのだろう？）こちらを参考にして$k = 0.05$，$B = 7$で試してみた．

t = var('t')
x = var('x')
y = var('y')
k = 0.05
b = 7
P = desolve_system_rk4([y,-k*y - x*x*x + b*cos(t)],[x,y],ics=[0,0,0],ivar=t,end_points=30000,step=0.1047196)
list_plot([(P[60*i][1],P[60*i][2]) for i in range(len(P)/60)])

なんかいい感じ．60個に一つ，$(x(t),y(t))$を出力している．stepは$\pi / 30$なので，$2\pi$に一つ．ちなみにかなり時間がかかる．

CPU times: user 27min 9s, sys: 24.8 s, total: 27min 34s
Wall time: 29min 4s

後でもう少し早い方法でやる．