[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).

download artikel

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.

[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.

[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 :

Map	Time domain	Space domain	Number of space dimensions	Also known as
Arnold's cat map	discrete	real	2
Baker's map	discrete	real	2
Bogdanov map
Chen-Lee system	continuous	real	3
Chossat-Golubitsky symmetry map
Circle map	discrete	real	1
Complex quadratic map	discrete	complex	1
Complex squaring map	discrete	complex	1
Complex Cubic map
Degenerate Double Rotor map
Double Rotor map
Duffing map	discrete	real	2
Duffing equation	continuous	real	1
Dyadic transformation	discrete	real	1	2x mod 1 map, Bernoulli map, doubling map, sawtooth map
Exponential map	discrete	complex	2
Gauss map	discrete	real	1	mouse map, Gaussian map
Generalized Baker map
Gingerbreadman map	discrete	real	2
Gumowski/Mira map
Hénon map	discrete	real	2
Hénon with 5th order polynomial
Hitzl-Zele map
Horseshoe map	discrete	real	2
Ikeda map	discrete	real	2
Interval exchange map	discrete	real	1
Kaplan-Yorke map	discrete	real	2
Linear map on unit square
Logistic map	discrete	real	1
Lorenz attractor	continuous	real	3
Lorenz system's Poincare Return map
Lozi map	discrete	real	2
Nordmark truncated map
Pomeau-Manneville maps for intermittent chaos	discrete	real	1 and 2	Normal-form maps for intermittency (Types I, II and III)
Pulsed rotor
Quasiperiodicity map
Rabinovich-Fabrikant equations	continuous	real	3
Random Rotate map
Rössler map	continuous	real	3
Shobu-Ose-Mori piecewise-linear map	discrete	real	1	piecewise-linear approximation for Pomeau-Manneville Type I map
Sinai map
Symplectic map
Standard map, Kicked rotor	discrete	real	2	Chirikov standard map, Chirikov-Taylor map
Tangent map
Tent map	discrete	real	1
Tinkerbell map	discrete	real	2
Triangle map
Van der Pol oscillator	continuous	real	1
Zaslavskii map	discrete	real	2
Zaslavskii rotation map

Kita akan me-review beberapa chaotic map, antara lain : Lorentz system, Henon map, Rossler 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.

[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 :

http://www.4shared.com/office/IB4Mhjx4/1D-FDTD_using_Matlab.html 

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.