ジャパニーズ・アトラクタを描きたくなったので,SageMathでやってみた.SageMathには一階の常微分方程式系を解くメソッドが用意されている.つまり,次の未知関数に関する微分方程式系
を解くことができる.まずは四次の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は各関数の変数.
- icsは初期値からなるリスト.ics[0]はの初期値で,それをとするとics[i]は.
- end_pointsはの最後の値.
- stepは微分を差分に直す時の隣同士のの差.省略すると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)
とする.
例えば を,初期値,,で解いてみる.(答えは言うまでも無く,である.
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)
での関数が出てくる.
list_plot([(x1,y1) for t1,x1,y1 in P],plotjoined=True)
だと円.
問題のジャパニーズ・アトラクタだけど,Wikipediaにあるパラメータだとうまくでてくれないので(何が悪いのだろう?)こちらを参考にして,で試してみた.
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個に一つ,を出力している.stepはなので,に一つ.ちなみにかなり時間がかかる.
CPU times: user 27min 9s, sys: 24.8 s, total: 27min 34s Wall time: 29min 4s
後でもう少し早い方法でやる.
0 件のコメント:
コメントを投稿
コメントの追加にはサードパーティーCookieの許可が必要です