Tulisan ini menunjukkan bagaimana cara menyelesaikan atau mencari solusi persamaan diferensial menggunakan ode45 di Octave, Python dan Julia.
Octave
close all
clear
clc
alpha = 9;
beta = 14;
m0 = -1/7;
m1 = 2/7;
F = @(x) m1*x + (1/2)*(m0-m1)*(abs(x+1)-abs(x-1));
Chua = @(t,x) [alpha*(x(2)-F(x(1)));x(1)-x(2)+x(3);-beta*x(2)];
xo = -1.1833;
yo = 0.079516;
zo = 1.0099;
N = 80000;
step_size = 0.01;
tspan = (0:N-1)*step_size;
tic
[t,xyz] =ode45(Chua,tspan,[xo,yo,zo]);
toc
x = xyz(:,1);
y = xyz(:,2);
z = xyz(:,3);
save('Chua_Octave.mat','x','y','z')
clear
clc
alpha = 9;
beta = 14;
m0 = -1/7;
m1 = 2/7;
F = @(x) m1*x + (1/2)*(m0-m1)*(abs(x+1)-abs(x-1));
Chua = @(t,x) [alpha*(x(2)-F(x(1)));x(1)-x(2)+x(3);-beta*x(2)];
xo = -1.1833;
yo = 0.079516;
zo = 1.0099;
N = 80000;
step_size = 0.01;
tspan = (0:N-1)*step_size;
tic
[t,xyz] =ode45(Chua,tspan,[xo,yo,zo]);
toc
x = xyz(:,1);
y = xyz(:,2);
z = xyz(:,3);
save('Chua_Octave.mat','x','y','z')
Julia
Untuk Julia, jangan lupa menginstall "ODE" dengan mengetik Pkg.add("ODE") di Julia Command Line anda.
using ODE
alpha = 9.0;
beta = 14.0;
m0 = -1.0/7.0;
m1 = 2.0/7.0;
F(x) = m1*x + (1/2)*(m0-m1)*(abs(x+1)-abs(x-1));
Chua(t,x) = [alpha*(x[2]-F(x[1]));x[1]-x[2]+x[3];-beta*x[2]];
xo = -1.1833;
yo = 0.079516;
zo = 1.0099;
N = 80000;
step_size = 0.01;
tspan = (0:N-1)*step_size;
tic()
(t,pos) = ode45(Chua,[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("./matfile/Chua_Julia.mat", "w")
write(file,"x",x)
write(file,"y",y)
write(file,"z",z)
close(file)
alpha = 9.0;
beta = 14.0;
m0 = -1.0/7.0;
m1 = 2.0/7.0;
F(x) = m1*x + (1/2)*(m0-m1)*(abs(x+1)-abs(x-1));
Chua(t,x) = [alpha*(x[2]-F(x[1]));x[1]-x[2]+x[3];-beta*x[2]];
xo = -1.1833;
yo = 0.079516;
zo = 1.0099;
N = 80000;
step_size = 0.01;
tspan = (0:N-1)*step_size;
tic()
(t,pos) = ode45(Chua,[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("./matfile/Chua_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
alpha = 9.0
beta = 14.0
m0 = -1.0/7.0
m1 = 2.0/7.0
def F(x):
return m1*x + (1/2)*(m0-m1)*(np.abs(x+1)-np.abs(x-1));
def Chua(x,t):
return [alpha*(x[1]-F(x[0])),x[0]-x[1]+x[2],-beta*x[1]]
xo = -1.1833
yo = 0.079516
zo = 1.0099
N = 80000
step_size = 0.01
t = np.arange(0.0, N*step_size, step_size)
tic = tm.time()
xyz = odeint(Chua,[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('./matfile/Chua_Python.mat',save_to_mat)
import numpy as np
import time as tm
alpha = 9.0
beta = 14.0
m0 = -1.0/7.0
m1 = 2.0/7.0
def F(x):
return m1*x + (1/2)*(m0-m1)*(np.abs(x+1)-np.abs(x-1));
def Chua(x,t):
return [alpha*(x[1]-F(x[0])),x[0]-x[1]+x[2],-beta*x[1]]
xo = -1.1833
yo = 0.079516
zo = 1.0099
N = 80000
step_size = 0.01
t = np.arange(0.0, N*step_size, step_size)
tic = tm.time()
xyz = odeint(Chua,[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('./matfile/Chua_Python.mat',save_to_mat)
Plotting
Hasil dapat dilihat dengan menjalan kode octave berikut
load Chua_Octave.mat
figure
plot3(x,y,z,'r','LineWidth',1.05)
title('Lorenz attractor calculated by Octave')
xlim([-2.5 2.5])
ylim([-0.5 0.5])
zlim([-8 4])
view(32,32)
xlabel('x')
ylabel('y')
zlabel('z')
grid minor
box off
load Chua_Julia.mat
figure
plot3(x,y,z,'r','LineWidth',1.05)
title('Lorenz attractor calculated by Julia')
xlim([-2.5 2.5])
ylim([-0.5 0.5])
zlim([-8 4])
view(32,32)
xlabel('x')
ylabel('y')
zlabel('z')
grid minor
box off
load Chua_Python.mat
figure
plot3(x,y,z,'r','LineWidth',1.05)
title('Lorenz attractor calculated by Python')
xlim([-2.5 2.5])
ylim([-0.5 0.5])
zlim([-8 4])
view(32,32)
xlabel('x')
ylabel('y')
zlabel('z')
grid minor
box off
figure
plot3(x,y,z,'r','LineWidth',1.05)
title('Lorenz attractor calculated by Octave')
xlim([-2.5 2.5])
ylim([-0.5 0.5])
zlim([-8 4])
view(32,32)
xlabel('x')
ylabel('y')
zlabel('z')
grid minor
box off
load Chua_Julia.mat
figure
plot3(x,y,z,'r','LineWidth',1.05)
title('Lorenz attractor calculated by Julia')
xlim([-2.5 2.5])
ylim([-0.5 0.5])
zlim([-8 4])
view(32,32)
xlabel('x')
ylabel('y')
zlabel('z')
grid minor
box off
load Chua_Python.mat
figure
plot3(x,y,z,'r','LineWidth',1.05)
title('Lorenz attractor calculated by Python')
xlim([-2.5 2.5])
ylim([-0.5 0.5])
zlim([-8 4])
view(32,32)
xlabel('x')
ylabel('y')
zlabel('z')
grid minor
box off