Selasa, 28 Agustus 2012

[bs] Delta Function and It's Differential

delta-function didefinisikan sebagai :


dengan phi(x) adalah test-function. delta-function dapat juga dipikirkan sebagai fungsi yang memiliki sifat : 


Kurva delta-function tampak dalam gambar di bawah ini.


Kita juga boleh membayangkan delta-function sebagai limit suatu fungsi sederhana seperti di bawah ini :


karena,


dan

untuk epsilon > 0, sehingga



lihat Gambar di atas. delta-function biasa dipakai untuk mengungkapkan distribusi potensial listrik dalam atom, kemudian untuk mengungkapkan sebuah impuls (sebuah gaya yang bekerja pada waktu yang singkat).


Untuk membuat sebuah fungsi dalam Matlab, seperti yang biasa kita lakukan, adalah dengan memakai perintah function, namun kali ini kita tidak langsung memberinya input berupa angka melainkan sebuah variabel. Oleh karena ini kita juga akan menggunakan perintah syms untuk mendeklarasikan nama variabel.

file 1 : delta-function

function d=delta(x,epsi,p)
% delta dirac
d=(1/p)*(epsi./(epsi^2+x^2));


file 2 : compiler untuk menampilkan hasil

clear all
syms x epsi p
Y = delta(x,epsi,p)


Save dengan nama, sebagai contoh, d3. Setelah di eksekusi kita akan memperoleh keluaran berupa fungsi yang masih mengandung variabel, yakni x, epsi dan p.


Karena keluaran masih berupa fungsi, maka kita perlu melakukan substiusi nilai-nilai agar dapat mendapatkan kurva fungsi delta. Yakni dengan menggunakan perintah subs :

clear all
syms x epsi p % nama-nama variabel dalam fungsi

awal  = -2;   % batas kiri nilai-x
akhir =  2;   % batas kanan nilai-x
n = 0;
X    = awal;  % Nilai yang akan di substitusikan ke variable x.
Epsi1 = 0.01; % Nilai yang akan disubstitusikan ke variable epsi
Epsi2 = 0.02;
Epsi3 = 0.03;

Pi   = 3.14;  % Nilai yang akan disubstitusikan ke variabel p.

while X<=akhir; 
   n = n+1; 
   Y1 = subs(delta(x,epsi,p),{x,epsi,p},{X,Epsi1,Pi});
   Y2 = subs(delta(x,epsi,p),{x,epsi,p},{X,Epsi2,Pi});
   Y3 = subs(delta(x,epsi,p),{x,epsi,p},{X,Epsi3,Pi});
   Xaxis(n)  = X ;
   Y1axis(n) = Y1;
   Y2axis(n) = Y2;
   Y3axis(n) = Y3;
   X=X+0.02*(akhir-awal);
   end
plot(Xaxis,Y1axis,'b.-',Xaxis,Y2axis,'r.-',Xaxis,Y3axis,'g.-')
legend('epsilon = 0.01','epsilon = 0.02','epsilon = 0.03')
xlabel('X')
ylabel('delta function')


Fungsi Y (delta-function) kini, setelah disubstitusikan nilai akan dapat menghasilkan kurva :



Kemudahan yang ditawarkan perintah syms ini kita dapat melakukan proses matematis seperti : diferensial dan integral secara langsung tanpa harus menyelesaikannya secara eksak (di atas kertas) terlebih dahulu. Pada contoh ini kita akan melakukan proses diferensial. Selanjtnya, untuk memperoleh hasil numerik dari fungsi yang mengalami diferensial ini kita gunakan perintah substitusi. Contoh algoritmanya sebagai berikut.

clear all

syms x epsi p
awal  = -2;
akhir =  2;
n = 0;
X    = awal;
Epsi1 = 0.01;
Epsi2 = 0.02;
Epsi3 = 0.03;

Pi   = 3.14;

while X<=akhir; 
   n = n+1; 
   Y1 = subs(diff(delta(x,epsi,p),x),{x,epsi,p},{X,Epsi1,Pi});
   Y2 = subs(diff(delta(x,epsi,p),x),{x,epsi,p},{X,Epsi2,Pi});
   Y3 = subs(diff(delta(x,epsi,p),x),{x,epsi,p},{X,Epsi3,Pi});
   Xaxis(n)  = X ;
   Y1axis(n) = Y1;
   Y2axis(n) = Y2;
   Y3axis(n) = Y3;
   X=X+0.02*(akhir-awal);
   end
plot(Xaxis,Y1axis,'b.-',Xaxis,Y2axis,'r.-',Xaxis,Y3axis,'g.-')
legend('epsilon = 0.01','epsilon = 0.02','epsilon = 0.03')
xlabel('X')
ylabel('delta function')



Kurva differensial delta-function tampak dalam gambar di bawah ini.


Kekurangan dari penggunaan perintah syms ini adalah pemakaian waktu yang menjadi lebih panjang bila dibandingkan ketika kita langsung memberikan input berupa angka dan mengerjakan perhitunga diferensial secara eksak (di atas kertas) terlebih dahulu.

Selasa, 21 Agustus 2012

[bs] Review on A chaotic system with only one stable equilibrium

By: 
Xiong Wang and Guanrong Chen

Department of Electronic Engineering, City University of Hong Kong, Hong Kong SAR, China.

If you are given a simple three-dimensional autonomous quadratic system that has only one stable equilibrium, what would you predict its dynamics to be, stable or periodic? Will it be surprising if you are shown that such a system is actually chaotic? Although chaos theory for three-dimensional autonomous systems has been intensively and extensively studied since the time of Lorenz in the 1960s, and the theory has become quite mature today, it seems that no one would anticipate a possibility of finding a three-dimensional autonomous quadratic chaotic system with only one stable equilibrium. The discovery of the new system, to be reported in this Letter, is indeed striking because for a three-dimensional autonomous quadratic system with a single stable node-focus equilibrium, one typically would anticipate non-chaotic and even asymptotically converging behaviors. Although the new system is of non-hyperbolic type, therefore the familiar ˇ Si’lnikov homoclinic criterion is not applicable, it is demonstrated to be chaotic in the sense of having a positive largest Lyapunov exponent, a fractional dimension, a continuous broad frequency spectrum, and a period-doubling route to chaos.

Sumber artikel : download paper on arxiv

Setelah kita berlatih mengenai chaotic map dalam pembahasan sebelumnya. Ayo kita coba memahami  problem yang berserakan dalam jurnal-jurnal ilmiah. Here we go. Paper Xiong Wang dan Guanrong Chen membicarakan suatu sistem yang berbentuk :


dengan a adalah parameter yang ditentukan kemudian. Mengingat kita masih beginer (hehe), kita coba aja memproduksi hasil-hasil yang telah diperoleh tersebut. Sejauh mana kita bisa memfollow mereka ? That is our problem ;). Jangan lupa, x dot artinya dx/dt, y dot adalah dy/dt dan z dot adalah dz/dt.

Adapun susunan program yang bisa kita buat untuk melihat phase portrait dalam sistem persamaan di atas adalah sebagai berikut.

Untuk kasus a = 0:

File 1 : Persamaan sistem :

function dxdt=china(t,x);
% Xiong Wang dan Guanrong Chen
a = 0;
dxdt=[ x(2).*x(3)+a;
x(1).^2-x(2);
 1-4.*x(1)];

File 2 : compiler untuk menjalankan iterasi

clear all           % a clear start
t=0:0.01:100;       % time window
y0=[0.25;0.4;3];    % initial conditions
[t,Y] = ode45('china',t,y0);
%figure
%plot(t,Y(:,1),'LineWidth',2.5)       % t vs X-axis
%plot(Y(:,1),Y(:,3),'LineWidth',2.5)  % X-axis vs Z-axis
%plot(Y(:,1),Y(:,2),'LineWidth',2.5)  % X-axis vs Y-axis
%plot(Y(:,2),Y(:,3),'LineWidth',2.5)  % Y-axis vs Z-axis
plot3(Y(:,1),Y(:,2),Y(:,3),'LineWidth',2.5) % X-Y-Z-axis
%xlabel('time')
xlabel('X-AXIS')
ylabel('Y-AXIS')
zlabel('Z-AXIS')
view([5,21])  %diaktifkan hanya untuk plot 3 dimensi
grid on
%figure
%comet3(Y(:,1),Y(:,2),Y(:,3))



Hasil komputasinya tampak dalam gambar di bawah ini :



dan


Untuk a = 0.006, file 1 kita modifikasi menjadi :

function dxdt=china(t,x);
% Xiong Wang dan Guanrong Chen
a = 0.006;

dxdt=[ x(2).*x(3)+a;
x(1).^2-x(2);
 1-4.*x(1)];

File 2 tidak mengalami modifikasi. Hasil komputasi memakai a = 0.006 tampak dalam gambar di bawah ini.

dan


Untuk figure-figure yang lain silakan dicoba sendiri ya. Aye juga belum bisa soalnya :D.








Senin, 20 Agustus 2012

[bs] Review on Chaotic Map

Dalam matematika, peta-chaos adalah peta (fungsi evolusi) yang menunjukkan beberapa jenis perilaku chaos. Peta ini dapat berparameterkan waktu secara diskrit ataupun kontinu. Peta diskrit biasanya mengambil bentuk fungsi iterasi. Peta chaos sering muncul dalam studi sistem dinamis.

Peta chaotic sering menghasilkan fraktal. Meskipun fraktal dapat dibangun oleh prosedur iterasi, beberapa fraktal yang dipelajari adalah berupa suatu himpunan dan bukan fungsi pemetaan yang menghasilkan mereka. Hal ini sering terjadi karena ada beberapa prosedur iteratif yang berbeda untuk menghasilkan fraktal yang sama. Daftar chaotic map :

MapTime domainSpace domainNumber of space dimensionsAlso known as
Arnold's cat mapdiscretereal2
Baker's mapdiscretereal2
Bogdanov map
Chen-Lee systemcontinuousreal3
Chossat-Golubitsky symmetry map
Circle mapdiscretereal1
Complex quadratic mapdiscretecomplex1
Complex squaring mapdiscretecomplex1
Complex Cubic map
Degenerate Double Rotor map
Double Rotor map
Duffing mapdiscretereal2
Duffing equationcontinuousreal1
Dyadic transformationdiscretereal12x mod 1 map, Bernoulli map, doubling map, sawtooth map
Exponential mapdiscretecomplex2
Gauss mapdiscretereal1mouse map, Gaussian map
Generalized Baker map
Gingerbreadman mapdiscretereal2
Gumowski/Mira map
Hénon mapdiscretereal2
Hénon with 5th order polynomial
Hitzl-Zele map
Horseshoe mapdiscretereal2
Ikeda mapdiscretereal2
Interval exchange mapdiscretereal1
Kaplan-Yorke mapdiscretereal2
Linear map on unit square
Logistic mapdiscretereal1
Lorenz attractorcontinuousreal3
Lorenz system's Poincare Return map
Lozi mapdiscretereal2
Nordmark truncated map
Pomeau-Manneville maps for intermittent chaosdiscretereal1 and 2Normal-form maps for intermittency (Types I, II and III)
Pulsed rotor
Quasiperiodicity map
Rabinovich-Fabrikant equationscontinuousreal3
Random Rotate map
Rössler mapcontinuousreal3
Shobu-Ose-Mori piecewise-linear mapdiscretereal1piecewise-linear approximation for Pomeau-Manneville Type I map
Sinai map 
Symplectic map
Standard mapKicked rotordiscretereal2Chirikov standard map, Chirikov-Taylor map
Tangent map
Tent mapdiscretereal1
Tinkerbell mapdiscretereal2
Triangle map
Van der Pol oscillatorcontinuousreal1
Zaslavskii mapdiscretereal2
Zaslavskii rotation map

Kita akan me-review beberapa chaotic map, antara lain : Lorentz systemHenon mapRossler attractor, dan Van der Pol oscillator. Mengingat Lorentz system dan Henon map sudah pernah kita bahas maka yang akan kita bahas selanjutnya adalah Rossler attractor dan Van der Pol oscillator.

Persamaan Rossler attractor :


dengan parameter: a = 0.2, b = 0.2, dan c = 5.7. Namun parameter yang sering dipakai sekarang ini adalah a = 0.1, b = 0.1 dan c = 14. Silakan dicoba sebagai latihan :D.

Untuk melihat phase portrait dari Rossler attractor kita akan menggunakan metode numerik yang dinamakan Runge-Kutta. Bisa juga memakai ode45 seperti biasa, silakan dicoba sebagai latihan. Berikut ini list programnya :

clear all;
n = 0;
x = 0;     %initial value sumbu-x
y = 0;     %initial value sumbu-y
z = 0;     %initial value sumbu-z
m = -2;    %parameter
a = 0.2;   %parameter
b = 0.2;   %parameter
c = 5.7;   %parameter

while n <= 50000;
    n = n+1;    
    dt =10^m;
    x1 = x;
    y1 = y;
    z1 = z;
    
    % fungsi untuk x
    kx1 = (-y1-z1).*dt;
    kx2 = (-y1-z1+0.5*kx1).*dt;
    kx3 = (-y1-z1+0.5*kx2).*dt;
    kx4 = (-y1-z1+kx3).*dt;
    x = x1+(kx1+2*kx2+2*kx3+kx4)./6;  

    % fungsi untuk y
    ky1 = (x1+a*y1).*dt;
    ky2 = (x1+a*(y1+0.5*ky1)).*dt;
    ky3 = (x1+a*(y1+0.5*ky2)).*dt;
    ky4 = (x1+a*(y1+ky3)).*dt;
    y = y1+(ky1+2*ky2+2*ky3+ky4)./6; 
    
    % fungsi untuk z
    kz1 = (b+z1*(x1-c)).*dt;
    kz2 = (b+(z1+0.5*kz1)*(x1-c)).*dt;
    kz3 = (b+(z1+0.5*kz2)*(x1-c)).*dt;
    kz4 = (b+(z1+kz3)*(x1-c)).*dt;
    z = z1+(kz1+2*kz2+2*kz3+kz4)./6; 
    
    xaxis(n) = x;
    yaxis(n) = y;
    zaxis(n) = z;
end
    
plot3(xaxis,yaxis,zaxis,'b-')
legend('Rossler system using Runge-Kutta')
xlabel('x-axis')
ylabel('y-axis')
zlabel('z-axis')

Hasil komputasi tampak dalam gambar di bawah ini.


Persamaan Van der Pol oscillator :


dengan nilai parameter miu bervariasi dari 0.01 sampai 4.0. Untuk melihat phase portrait Van der Pol oscillator, kita akan memakai ode45 sebagai berikut.

File 1:  untuk menuliskan fungsi Van der Pol oscillator

function dxdt=vanderpol(t,x);
% Van der Pol oscillator
miu = 1.5;
dxdt=[x(2);
    miu.*(1-x(1).^2).*x(2)-x(1)]; 

File 2: compiler untuk menjalankan iterasi

clear all    % a clear start
t=[0 100];   % time window
y0=[0;0.2];  % initial conditions x0,v0
[t,Y] = ode45('vanderpol',t,y0);
figure
%plot3(t,Y(:,1),Y(:,2))
plot(Y(:,1),Y(:,2))
%plot(t,Y(:,1))
%view([3,16])
grid on
figure
%comet3(Y(:,1),Y(:,2),Y(:,3))
comet(Y(:,1),Y(:,2))
%comet(t,Y(:,1))
%xlabel('time')
%ylabel('y-axis')
%zlabel('z-axis')

Hasil komputasi tampak dalam gambar di bawah ini.


Kamis, 16 Agustus 2012

[bs] Review on 1 D FDTD using MATLAB

Kita lagi mupeng banget sama metode ini, Finite Difference Time Domain (FDTD) method . Sebagai pengantar, kita akan mencoba yang simple, level 1 dimensi. Kelak kita pengen juga mencoba yang level 3 Dimensi. Bagai rekan-rekan yang sudah bisa, share dong gimana caranya ? Atau bagi rekan-rekan yang berminat belajar bareng, yuk, kita bahas sama-sama. 


1D-FDTD using MATLAB

By Hui Loui, Student Member, IEEE

Abstract—This report presents a simple 1D implementation of the Yee FDTD algorithm using the MATLAB programming language. The fields Ex and Hy are simulated along the line X = Y = 0, i.e. propagation along the ˆz axis. Source implementation and the effects of various boundaries such as PEC, PMC, Mur on the incident/scattered/total fields are subsequently investigated. The goal of this project is to exercise the basic parts of a FDTD code in the simplest system.

Index Terms—1D FDTD, Gaussian Pulse.

Sumber artikel :


File  1 : Fungsi medan listrik sumber

function hitung = source(t);
lambda = 0.5*10^(-6);
epsl0 = 8.854*10^(-12);
mu0 = 4*pi*10^(-7);
c = 1/sqrt(epsl0*mu0);
f0 = c/lambda;
tau0 = 1/f0;
tc = 5*tau0/2;
dt = tau0/20;
dz = c*dt;
sigma = tc/(2*sqrt(2*log(2)));
m = sin(2*pi*f0*t);
g = exp(-((t-tc)^2/(2*sigma^2)));
hitung = m*g;


File 2 : fungsi gaussian envelope

function hitung = source3(t);
lambda = 0.5*10^(-6);
epsl0 = 8.854*10^(-12);
mu0   = 4*pi*10^(-7);
c     = 1/sqrt(epsl0*mu0);
f0    = c/lambda;
tau0  = 1/f0;
tc    = 5*tau0/2;
dt    = tau0/20;
dz    = c*dt;
sigma = tc/(2*sqrt(2*log(2)));
m     = sin(2*pi*f0*t);
hitung = exp(-((t-tc)^2/(2*sigma^2)));


File 3 : compiler

clear all
lambda = 0.5*10^(-6);
epsl0 = 8.854*10^(-12);
mu0   = 4*pi*10^(-7);
c     = 1/sqrt(epsl0*mu0);
f0    = c/lambda;
tau0  = 1/f0;
awal  = 0;
batas = 20*tau0;
n     = 0;
t     = awal;
dt    = tau0/20;

while t<=batas; 
   n = n+1; 
   
   X  = t; 
   Y1  = source(t);
   Y3  = source3(t); 
   Xaxis(n)  = t ;
   Y1axis(n) = Y1;
   Y3axis(n) = Y3;
   t=t+dt;
   
end

plot(Xaxis,Y1axis,'b-',Xaxis,Y3axis,'r-')
xlabel('t(s)')
ylabel('Ex(t)')
legend('Modulated Pulse','Gaussian Envelope')

Matlab untuk 1D FDTD terdiri dari dua file, antara lain :

File 1: Fungsi Pulsa Listrik


function hitung = source4(n);
lambda = 0.5*10^(-6);
epsl0 = 8.854*10^(-12);
mu0   = 4*pi*10^(-7);
c     = 1/sqrt(epsl0*mu0);
f0    = c/lambda;
tau0  = 1/f0;
tc    = 5*tau0/2;
dt    = tau0/20;
sigma = tc/(2*sqrt(2*log(2)));
m     = sin(2*pi*f0*n*dt);
g     = exp(-((n*dt-tc)^2/(2*sigma^2)));
hitung = m*g;

File 2 : FDTD


clear all
lambda = 0.5*10^(-6);
epsl0 = 8.854*10^(-12);
mu0   = 4*pi*10^(-7);
c     = 1/sqrt(epsl0*mu0);
f0    = c/lambda;
tau0  = 1/f0;
dt    = tau0/20;
dz    = lambda/20;
M     = 500;    %space
N     = 300;    %time
z     =0:dz:M*dz; 
zz    = 0:1:M;
yy    = 0:1:M-1;
Hy(1:M)  =0;
Ex(1:M+1)=0;
for n=1:N,
Ex(1)=source4(n);
Hy = Hy-(dt./mu0).*(diff(Ex)./dz);
Ex(2:M)=Ex(2:M)-(dt./epsl0).*(diff(Hy)./dz);
end
%plot(250,0,'.b',zz,Ex,'r-')
%plot(yy,Hy,'b-')
plot(250,0,'b.',zz,Ex,'r-','MarkerSize',15)
%plot(250,0,'b.',yy,Hy,'g-','MarkerSize',15)
xlabel('k (node step)')
ylabel('Ex')
%legend('t = 0 fs')

Hasil komputasi tampak dalam gambar di bawah ini.