ジャパニーズ・アトラクタを描きたくなったので,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.
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)])
![](https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgQid_S204whJ8yT23G_jd67mDjpTMXy3dciScre-BY4SfNDrUNlXdvhGwHLAwHVkJv0536pjtUc4eohe5eMZ_pEEPxajNDBh-Uh2NjxKpk_Ob2VSTtaFtUHqOeUKjhu7ao3iLfLbkQM-A/s640/ja.png)
なんかいい感じ.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
後でもう少し早い方法でやる.
0 件のコメント:
コメントを投稿
コメントの追加にはサードパーティーCookieの許可が必要です