Rossler Attractor

Tulisan ini menunjukkan bagaimana cara menyelesaikan atau mencari solusi persamaan diferensial menggunakan ode45 di Octave, Python dan Julia.

Octave

close all
clear
clc

a = 0.1;
b = 0.1;
c = 14;

Rossler = @(t,x) [-x(2)-x(3); x(1)+a*x(2); b+(x(1)-c)*x(3)];

xo = 1.1930;
yo = 14.768;
zo = 0.0686;

N = 50000;
step_size = 0.01;
tspan = (0:N-1)*step_size;

tic
[t,xyz] =ode45(Rossler,tspan,[xo,yo,zo]);
toc

x = xyz(:,1);
y = xyz(:,2);
z = xyz(:,3);

save('Rossler_Octave.mat','x','y','z')

Julia

Untuk Julia, jangan lupa menginstall "ODE" dengan mengetik Pkg.add("ODE") di Julia Command Line anda.

using ODE

a = 0.1;
b = 0.1;
c = 14;

Rossler(t,x) = [-x[2]-x[3]; x[1]+a*x[2]; b+(x[1]-c)*x[3]];

xo = 1.1930;
yo = 14.768;
zo = 0.0686;

N = 50000;
step_size = 0.01;
tspan = (0:N-1)*step_size;

tic()
(t,pos) = ode45(Rossler,[xo,yo,zo],tspan,points=:specified)
toc()

x = map(v->v[1],pos)
y = map(v->v[2],pos)
z = map(v->v[3],pos)

using MAT

file = matopen("Rossler_Julia.mat", "w")
write(file,"x",x)
write(file,"y",y)
write(file,"z",z)
close(file)

Python

from scipy.integrate import odeint
import numpy as np
import time as tm

a = 0.1
b = 0.1
c = 14.0

def Rossler(x,t): 
  return [-x[1]-x[2], x[0]+a*x[1], b+(x[0]-c)*x[2]];

xo = 1.1930
yo = 14.768
zo = 0.0686

N = 50000
step_size = 0.01
t = np.arange(0.0, N*step_size, step_size)

tic = tm.time()
xyz = odeint(Rossler,[xo,yo,zo],t)
toc = tm.time()

print('Elapsed time: %.4f seconds' % (toc-tic))

import scipy.io
save_to_mat = {'x':xyz[:,0],'y':xyz[:,1],'z':xyz[:,2]}
scipy.io.savemat('Rossler_Python.mat',save_to_mat)

Plotting

Hasil dapat dilihat dengan menjalan kode octave berikut

load Rossler_Octave.mat

figure
plot3(x,y,z,'r','LineWidth',1.05)
title('Lorenz attractor calculated by Octave')
xlim([-20 25])
ylim([-25 20])
zlim([-20 40])
view(28,28)
xlabel('x')
ylabel('y')
zlabel('z')
grid minor
box off

load Rossler_Julia.mat

figure
plot3(x,y,z,'r','LineWidth',1.05)
title('Lorenz attractor calculated by Julia')
xlim([-20 25])
ylim([-25 20])
zlim([-20 40])
view(28,28)
xlabel('x')
ylabel('y')
zlabel('z')
grid minor
box off

load Rossler_Python.mat

figure
plot3(x,y,z,'r','LineWidth',1.05)
title('Lorenz attractor calculated by Python')
xlim([-20 25])
ylim([-25 20])
zlim([-20 40])
view(28,28)
xlabel('x')
ylabel('y')
zlabel('z')
grid minor
box off