BAB III METODE PENELITIAN
3.1. Perancangan Sistem dan Blok Diagram Sistem Model penelitian yang akan dilakukan adalah model penelitian pengembangan. Untuk mempermudah dalam memahami sistem yang akan dibuat dapat dijelaskan melalui blok diagram pada Gambar 3.1.
Citra Sinar X Tulang Belakang
Pre-processing
Segmentasi menggunakan CPM
Citra hasil
Gambar 3.1. Blok Diagram
Penelitian yang dilakukan memiliki 2 tahap utama, yaitu bagian preprocessing dan segmentasi menggunakan CPM. 1. Pre-processing Pre-processing yaitu proses untuk filtering citra. filtering digunakan untuk memperjelas bagian gambar tulang punggung yang biasanya kurang terlihat jelas, dikarenakan tertutup tulang dada dan organ-organ tubuh yang 22
23
berada di sekitar daerah dada. Pre-processing yang digunakan yaitu Modified Tophat filter dan Gaussian cropping. Modified Tophat filter digunakan untuk mengurangi noise pada daerah tulang belakang. Gaussian cropping digunakan untuk menghilangkan dan mengurangi noise dari bagian-bagian citra sinar x yang tidak diperlukan. 2. Segmentasi menggunakan charged particle model (CPM) Segmentasi dilakukan untuk memperoleh segmentasi kelengkungan tulang belakang pada citra tulang belakang yang sudah selesai difilter. Setelah didapatkan citra hasil filtering, maka selanjutnya inisialisasi partikel-partikel di sekitar tulang belakang. Partikel-partikel ini yang nantinya akan bergerak berdasarkan gaya total yang diperoleh dari gaya Lorentz dan Coulomb dari citra.
3.2. Perancangan Perangkat Lunak Pada perancangan CPM ini, penulis menggunakan compiler Microsoft Visual C++ 2008 dan library yang digunakan adalah OpenCV v1.1, library ini menyediakan fungsi-fungsi yang akan digunakan untuk pengolahan citra. Dalam penulisan dan pembuatan program, akan meliputi bagian-bagian penting dalam setiap langkah-langkah per bagian sesuai dengan algoritma atau logika sekuensial dari awal sampai output. Flowchart program secara global dijelaskan pada gambar 3.2.
24
Start
Input Citra Sinar X Tulang Belakang
Modified Tophat filter
Gaussian Cropping
Hitung gaya Lorentz Inisialisasi partikel Iterasi dimulai dari 0
Y Iterasi <1000
Hitung gaya Coulomb Hitung gaya total Update posisi partikel Iterasi bertambah 1
T Citra hasil segmentasi menggunakan CPM
End
Gambar 3.2. Flowchart program Penjelasan flowchart: 1. Inputkan citra sinar x tulang belakang yang akan diproses. 2. Citra yang telah diinputkan selanjutnya diproses menggunakan Modified Tophat filter 3. Citra hasil Modified Tophat filter selanjutnya diproses menggunakan Gaussian cropping.
25
4. Hitung gaya Lorentz dari citra hasil Gaussian cropping, digunakan sebagai medan gaya negatif. 5. Selanjutnya proses inisialisasi partikel-partikel yang nantinya akan bergerak menuju tepi dari tulang belakang yang berada di citra. 6. Proses iterasi dimulai dari 0 dan akan berulang sebanyak 1000 iterasi. 7. Hitung gaya Coulomb dari partikel-partikel, kemudian kalkulasikan dengan gaya Lorentz untuk mendapatkan gaya total, dari gaya total akan didapatkan posisi baru dari tiap partikel. 8. Citra hasil segmentasi akan didapatkan setelah langkah ke-7 berulang sebanyak iterasi yang ditentukan.
3.2.1. Perancangan Modified Tophat Filter Tahap awal dari proses pre-processing adalah proses filtering untuk memperjelas citra tulang belakang dan mengurangi noise di sekitar tulang belakang menggunakan Modified Tophat filter. Rumus Modified Tophat filter adalah : ∘
modTH= Dimana :
•
f
× Ɣ ...........................................................................(3.1)
= citra gambar
b = struktur elemen ∘
= operator morfologi closing
•
= operator morfologi opening
Ɣ
= attenuation factor
Alur proses Modified Tophat filter bisa dilihat pada flowchart dibawah ini :
26
Start
Input Citra Sinar X (f)
Inisialisasi struktur elemen (B) = disk(5) dan Attenuation Factor (Ɣ) = 11
TH = f ∘ B (opening) CL = f • B (closing)
TH = TH – f
ModTH =
xƔ
Citra hasil modified tophat filter
End
Gambar 3.3. Flowchart Modified Tophat filter Jika diprogram menggunakan rumus Modified Tophat filter, maka dapat ditulis code program sebagai berikut : int disk = 5, factor = 11; norMat(f, f32Image, 1); IplConvKernel *B =cvCreateStructuringElementEx(2*disk+1, 2*disk+1, disk, disk, CV_SHAPE_ELLIPSE, 0); //B = disk(5) cvMorphologyEx(f32Image,TH,tempMat,B,CV_MOP_OPEN,1);//opening cvMorphologyEx(f32Image,CL,tempMat,B,CV_MOP_CLOSE,1);//closing cvSub(TH,f32Image,TH); //TH=TH-f cvDiv(TH, CL, modTH, factor); //ModTH=(TH/CL)*atten_factor norMat(modTH,modTH); cvCopyImage(modTH,result);
Hasil yang diperoleh ditunjukan pada gambar 3.4 :
27
Gambar 3.4. Hasil Modified Tophat filter Perancangan Modified Tophat filter dimulai dengan melakukan operasi opening , yaitu proses erosi dilanjutkan dengan dilasi pada citra gambar dengan struktur elemen bertipe disk dengan ukuran 5. Hasil dari proses opening kemudian dikurangkan dengan citra awal dari citra gambar sehingga menyisakan tepi dari objek saja. Selanjutnya dilakukan proses closing menggunakan citra awal, yaitu proses dilasi yang dilanjutkan dengan erosi dengan menggunakan struktur elemen berukuran sama. Proses opening menyebabkan citra menjadi lebih mengembang, sebaliknya proses closing menyebabkan citra menyusut. Kemudian hasil dari opening dikurangi citra awal dibagi dengan hasil closing untuk menghilangkan sisa noise selanjutnya dikalikan dengan attenuation factor untuk memperjelas bentuk tepi.
28
3.2.2. Perancangan Gaussian Cropping Tahap kedua dari proses pre-processing adalah dengan menghilangkan dan mereduksi noise pada bagian-bagian dari citra sinar x yang tidak diperlukan khususnya bagian kanan dan kiri dari citra tulang belakang menggunakan Gaussian cropping. Rumus Gaussian cropping adalah : (
Gaussian =
(
) )
...........................................................................(3.2)
Keterangan : a = koordinat piksel pada citra b = koordinat center of peak dari fungsi Gaussian c = standar deviasi
Gambar 3.5. Kurva fungsi Gaussian (Sugianto, 2013) Alur dari penggunaan Gaussian cropping bisa dilihat pada gambar 3.6 :
29
Start
Input Citra hasil ModTH
Inisialisasi garis Simpan koordinat dari setiap garis ( , ) Y=0
T Y
Y T Y == Nilai piksel = 0
Y (
Gaussian =
(
) )
Nilai piksel = nilai piksel * gaussian
Citra hasil segmentasi menggunakan CPM
End
Gambar 3.6 Flowchart Gaussian cropping Jika diprogram menggunakan rumus Gaussian cropping, maka dapat ditulis kode program sebagai berikut :
30
int counter=0; int jmlh_baru=jmlh; //baris atas for (int i=1;i<=koordinat[1][2];i++) for (int j=0;j<=edge->width;j++) edge->imageData[edge->widthStep*i + j*edge->nChannels]=0; //baris tengah int index=1; for (int y=koordinat[1][2];y<=koordinat[jmlh_baru][2];y++){ if (y==koordinat[index][2]){ y=koordinat[index][2]; float nilai=0; for (int x=0;x<=edge->width;x++){ double gaussian=0; gaussian=((koordinat[index][1]-x) *(koordinat[index][1]-x)) / ((2*10)*(2*10)); gaussian=exp(-(gaussian)); nilai=(uchar)edge -> imageData[edge->widthStep*y + x*edge->nChannels]; edge->imageData[edge->widthStep * y + x* edge->nChannels]=nilai* gaussian; } index=index+1; }else{ for (int x=0;x<edge->width;x++){ edge->imageData[edge->widthStep * y + x*edge->nChannels] =(uchar)edge ->imageData [edge->widthStep*(y-1) + x*edge->nChannels]; } } } //baris bawah for (int y=koordinat[jmlh_baru][2];y<=edge->height-1;y++) for (int x=0;x<=edge->width;x++) edge->imageData[edge->widthStep*y + x*edge->nChannels]=0; cvDestroyWindow("Gaussian Filter"); cvNamedWindow("Gaussian Filter",1); cvShowImage("Gaussian Filter",edge); //kosongkan SLL while (!InitContour.empty()){ InitContour.pop_back(); } std::cout<<"kosong"<<std::endl;
Kode program untuk proses inisialisasi diletakkan pada function terpisah, ditulis sebagai berikut: void gausian(IplImage * img) { cvCopyImage(img,temp); int k=0; int l=0; for (int x=0;x
31
cvShowImage( "srcImage", temp ); k=l; l=l+1; } } void filter( int event, int x, int y, int flags, void* ptr) { if(event == CV_EVENT_LBUTTONDOWN) { gausian((IplImage *)ptr); jmlh=jmlh+1; InitContour.push_back(cvPoint(x,y)); koordinat[jmlh][1]=x; koordinat[jmlh][2]=y; q=1; } if (q==1) { if (event==CV_EVENT_MOUSEMOVE) { gausian((IplImage *)ptr); jmlh=jmlh+1; InitContour.push_back(cvPoint(x,y)); koordinat[jmlh][1]=x; koordinat[jmlh][2]=y; } } }
Diawali dengan membuat inisialisasi pada citra tulang belakang menggunakan gerakan mouse dan menyimpan setiap titik yang dilewati oleh mouse. Untuk setiap titik tersebut akan menjadi titik puncak dari fungsi Gaussian sehingga semakin ke kanan dan ke kiri dari koordinat x dari titik puncak akan menghasilkan nilai Gaussian yang mendekati nilai 0, nilai Gaussian tersebut kemudian dikalikan dengan nilai piksel awal untuk menghasilkan nilai piksel yang baru. Penentuan nilai variabel c disesuaikan dengan lebar dari tulang belakang. Hasil yang diperoleh ditunjukan oleh Gambar 3.7.
32
Gambar 3.7 Hasil Gaussian cropping dengan nilai c=10
3.2.3. Perancangan Charged Particle Model (CPM) Perancangan CPM meliputi perancangan negative field yang diaplikasikan sebagai gaya Lorentz citra sinar x dan partikel positif yang diaplikasikan sebagai gaya Coulomb. Setelah gaya Lorentz dan Coulomb didapat maka diteruskan dengan perhitungan pergerakan partikel. Partikel bergerak sebanyak 1000 iterasi, dimana tiap iterasinya dihitung nilai Coulomb, gaya total, kecepatan partikel, percepatan partikel, dan posisi partikel yang baru.
A. Perancangan Negatif Field Gaya Lorentz pada citra berfungsi sebagai medan gaya negatif. Hasil perhitungan dari medan gaya negatif dan partikel positif ini yang menunjukkan arah partikel positif untuk bergerak. Perhitungan gaya Lorentz menggunakan
33
rumus yang ada pada persamaan (2.7). Alur dari perancangan gaya Lorentz sesuai dengan flowchart pada gambar 3.8.
Gambar 3.8. Flowchart menghitung gaya Lorentz Penjelasan flowchart: 1. Citra input didefinisikan sebagai variabel f. 2. Hitung Gradient Map dari f dan simpan dengan variabel (gx,gy). 3. Setelah mendapat (gx,gy), langkah selanjutnya adalah mencari nilai magnitude didefinisikan sebagai variabel e. 4. Dari variabel e langkah selanjutnya adalah menghitung energi potensial. 5. Hitung Gradient Map dari energi potensial hasil dari langkah sebelumnya, simpan ke variabel (Ex,Ey).
34
6. Langkah terakhir adalah menghitung gaya Lorentz (FLx,FLy). Yang jika diimplementasikan dalam bentuk program ditulis seperti di bawah ini : IplImage *src; CvMat *CH, *fx, *fy, *X, *Y, *x, *y; char gam[100], lok[100]; double ip,is,in,itot; //load image grayscale --> src strcpy(lok, "C:/2/fr10"); strcpy(gam, lok); strcat(gam, ".jpg"); src = cvLoadImage(gam, CV_LOAD_IMAGE_GRAYSCALE); CH = cvCreateMat(m, n, CV_32FC1); fx = cvCreateMat(m, n, CV_32FC1); fy = cvCreateMat(m, n, CV_32FC1); X = cvCreateMat(m, n, CV_32FC1); Y = cvCreateMat(m, n, CV_32FC1); x = cvCreateMat(1, n, CV_32FC1); y = cvCreateMat(1, m, CV_32FC1); //menghitung gradient map terhadap sumbu x citra input --> for(int j=0; j<m; j++){ for(int i=0; i
imageData [src->widthStep*j + (i+1)*src->nChannels]); is = ((double)src -> imageData [src->widthStep*j + i*src->nChannels]); in = ((double)src -> imageData [src->widthStep*j + (i-1)*src->nChannels]); if(i==0){ itot = ip-is; }else if(i==n-1){ itot = is-in; }else{ itot = ((ip - is) + (is - in)) / 2; } cvmSet(fx, j, i, itot); } } //menghitung gradient map terhadap sumbu y citra input --> for(int j=0; j<m; j++){ for(int i=0; i imageData [src->widthStep*(j+1) i*src->nChannels]); is = ((double)src -> imageData [src->widthStep*j + i*src->nChannels]); in = ((double)src -> imageData [src->widthStep*(j-1) i*src->nChannels]); if(j==0){ itot = ip-is; }else if(j==m-1){ itot = is-in; }else{ itot = ((ip - is) + (is - in)) / 2; } cvmSet(fy, j, i, itot);
fx
fy
+
+
35
} } //menghitung magnitude --> CH for (int j=0; j<m; j++){ for (int i=0; i
Program diatas diletakkan pada void main(). Pada penjelasan flowchart nomor 4 dan 5 diatas dijadikan satu menjadi persamaan (2.6). Kemudian persamaan (2.6) disubstitusikan dengan persamaan (2.5) didapatkan persamaan (2.7) yang sudah mencakup langkah 4 sampai 6 pada penjelasan flowchart diatas, kode program tersebut diletakkan pada dua function berikut: double LorentzForceX (CvMat* X, CvMat* Y, CvMat* CH, double x, double y, int I, int J) double LorentzForceY (CvMat* X, CvMat* Y, CvMat* CH, double x, double y, int I, int J)
Kode program dari function diatas adalah sebagai berikut: for (int s=0; s<m; s++){ for (int t=0; t
36
a = cvmGet(X, b = cvmGet(Y, cvmSet(Rx, s, cvmSet(Ry, s,
s, s, t, t,
t); t); a-x); b-y);
} }//ri-Rk......(1) for (int s=0; s<m; s++){ for (int t=0; t
37
double h = cvmGet(R, s, I); double i = f/h; fh=fh+i; }//jumlahkan semua nilai tiap piksel pada kolom ke=i //baris ke-(J+1) sampai m Ex=Ex+ac+fh; //jumlahkan kedua nilai di atas ditambah dengan nilai //Ex sebelumnya return(Ex); //nilai Ex sebagai nilai akhir dari gaya LorentzX //Ey untuk LorentzY
Perhitungan nilai gaya Lorentz diatas diujicobakan pada citra berukuran 10x20 piksel seperti yang ditunjukkan pada Gambar 3.9.
Gambar 3.9. Citra berukuran 10x20 piksel Nilai tiap-tiap piksel yang didapat dari hasil perhitungan gaya Lorentz dari gambar 3.9 terhadap sumbu x adalah sebagai berikut : 88.0136 130.631 172.102 175.5 171.356 167.42 167.521 165.129 164.387 163.753 163.246 162.942 162.849 163.046 163.744 165.171 164.966 146.107 92.8302 58.3159
132.995 173.818 146.647 104.518 88.0211 82.3215 80.2915 79.4959 77.8843 77.1663 76.7536 76.5766 76.718 77.4008 79.3162 84.7014 101.394 138.453 123.263 68.691
143.097 99.2419 79.4903 54.2632 41.8488 35.8796 32.6223 31.0238 29.1481 28.1702 27.7391 27.6421 27.9146 28.8055 30.9802 36.1198 47.5173 63.6729 62.201 54.4202
68.1244 -15.3685 3.68569 -3.33751 -9.56371 -14.0223 -16.3481 -18.3915 -19.0991 -19.8747 -20.2475 -20.2936 -20.0239 -19.2402 -17.3622 -12.6919 -9.91632 -42.5529 -26.7152 14.6163
0.614066 -31.8323 -55.6242 -81.0894 -91.4246 -93.4487 -95.0389 -95.754 -97.118 -97.6683 -97.8299 -97.8277 -97.6736 -97.2874 -96.3358 -92.9475 -84.8399 -81.7438 -32.6177 -12.2017
-20.3432 -0.30634 -9.41755 -0.77267 2.34226 4.26558 5.21113 5.40245 5.17375 5.5132 5.65636 5.70693 5.69847 5.52523 4.7558 1.47897 -4.52043 4.36892 -13.2994 -17.6183
-102.779 -31.5865 -55.1135 -42.7589 -34.1241 -28.6939 -26.5801 -26.406 -24.7204 -24.3757 -24.1892 -24.0953 -24.2137 -24.7741 -26.4325 -31.3752 -37.3974 -12.5231 -22.4757 -43.2968
-161.597 -129.011 -103.836 -72.8272 -58.9244 -52.4371 -49.958 -47.6149 -46.008 -45.1108 -44.6319 -44.4342 -44.5416 -45.1849 -47.0403 -52.0606 -63.6232 -77.3288 -78.8592 -67.7901
-118.691 -154.236 -128.764 -96.6623 -81.4581 -76.9728 -75.5929 -71.8976 -70.2287 -68.8029 -68.071 -67.7495 -67.7276 -68.149 -69.5292 -73.5178 -85.8848 -111.616 -98.6632 -57.9619
-75.7633 -103.792 -127.789 -127.571 -121.924 -120.273 -118.83 -117.938 -116.133 -115.023 -114.374 -113.974 -113.761 -113.765 -114.091 -114.751 -114.369 -104.245 -69.6802 -47.2275
38
Nilai tiap-tiap piksel yang didapat dari hasil perhitungan gaya Lorentz dari gambar 3.9 terhadap sumbu y adalah sebagai berikut : 69.3786 68.9034 48.236 19.4998 6.97394 5.12057 1.89259 -0.66823 -2.36739 -3.5954 -4.66227 -6.27698 -8.20501 -10.5414 -13.5661 -17.4475 -31.4903 -59.2013 -70.3938 -63.9405
118.156 114.166 55.4428 5.24815 1.29011 3.23435 1.54925 -2.40577 -2.41663 -3.91917 -5.00715 -6.39431 -8.00505 -9.76799 -11.5279 -13.2601 -26.8202 -92.4452 -128.038 -98.5591
169.682 78.5412 -13.3538 -31.7258 -8.84695 -1.32576 -0.44343 -3.40954 -4.78575 -4.62996 -5.29038 -6.37328 -7.5529 -8.4788 -7.98079 -2.67965 18.1646 -43.2519 -146.115 -145.558
142.362 42.3371 21.2765 -36.5684 -13.6314 -4.90363 -2.40664 -5.54383 -5.15023 -5.13204 -5.51207 -6.28524 -7.04176 -7.17894 -5.24634 2.58831 -0.04895 -80.5076 -135.668 -131.769
109.183 -0.29054 1.12856 -59.8018 -18.9553 -7.72763 -4.35352 -6.1188 -6.42459 -5.78748 -5.65172 -6.19681 -6.66924 -6.25712 -2.62297 13.7872 -5.10149 -55.4201 -80.9469 -116.266
105 -10.8891 -10.0006 -64.9066 -20.8897 -7.26096 -4.45807 -7.4167 -6.86254 -5.609 -5.64231 -6.14605 -6.57642 -6.14903 -2.76488 13.5791 -3.19249 -45.2283 -70.8441 -109.296
130.439 7.87279 -14.5801 -48.2962 -16.2796 -4.7884 -3.76507 -6.34561 -6.32679 -5.38198 -5.52327 -6.13213 -6.75987 -6.82892 -5.00239 3.20758 4.42238 -50.0695 -97.0844 -111.365
135.911 47.5936 -28.2751 -31.915 -9.53139 -0.12884 -2.84381 -5.55092 -4.77875 -5.19421 -5.28627 -6.10292 -7.07312 -7.87144 -7.76452 -4.47956 12.033 -30.5813 -107.008 -106.031
86.336 73.2699 29.0034 -2.06007 -0.00075 2.57686 0.146992 -2.28078 -3.88331 -4.09997 -4.86929 -5.98612 -7.32102 -8.82903 -10.5384 -13.0181 -19.0984 -61.8331 -88.9458 -68.9028
49.7954 44.832 27.9903 11.5047 6.02884 3.91392 1.70904 -0.33896 -2.05617 -3.16162 -4.36989 -5.74519 -7.34494 -9.29164 -11.9138 -16.1817 -24.6875 -41.6405 -50.2074 -44.7975
B. Perancangan Partikel Positif Partikel positif yang digunakan merupakan partikel yang diberi gaya Coulomb sehingga partikel bermuatan positif. Partikel positif akan saling tolak menolak dengan partikel positif lainnya, dan tarik menarik dengan medan negatif citra. Alur dari perancangan gaya Coulomb sesuai dengan flowchart pada gambar 3.10.
39
Start
Baca citra input f
[sx , sy] = inisialisasi partikel (f)
= =
− − − −
CoulombForce = (Fcx , Fcy)
End
Gambar 3.10. Flowchart Menghitung Gaya Coulomb Flowchart pada Gambar 3.10 jika diimplementasikan dalam bentuk program ditulis seperti berikut : int s=0, sx[100], sy[100]; double CoulombX[100], CoulombY[100]; for (int i=0;i<s;i++){ CoulombX[i]=0; CoulombY[i]=0; for (int j=0;j<s;j++){ if (j != i){ if ((sx[i]-sx[j] != 0) && (sy[i]-sy[j] != 0)){ CoulombX[i] = CoulombX[i] + ((sx[i]-sx[j]) / pow(Absolut(sx[i] - sx[j]), 3)); CoulombY[i] = CoulombY[i] + ((sy[i]-sy[j]) / pow(Absolut(party[i] - party[j]), 3)); }else if (partx[i] - partx[j] == 0){ CoulombY[i] = CoulombY[i] + ((sy[i]-sy[j]) / pow(Absolut(party[i] - party[j]), 3)); }else if (party[i] - party[j] == 0){ CoulombX[i] = CoulombX[i] + ((sx[i]-sx[j]) / pow(Absolut(sx[i] - sx[j]), 3)); } } } }
40
Kode program untuk proses inisialisasi diletakkan pada function terpisah, ditulis sebagai berikut: void showContent(IplImage * img){ cvCopyImage(img,temp); for (int i=0;i<=part;i++){ cvCircle(temp, cvPoint(sx[i],sy[i]), 2, CV_RGB(255,0,0)); } cvShowImage("CPM", temp); } void on_mouse( int event, int x, int y, int flags, void* ptr){ if (event == CV_EVENT_LBUTTONDOWN){ sx[part] = x; sy[part] = y; printf("%i %i %i\n",part,sx[part],sy[part]); part++; showContent((IplImage *)ptr); } }
Untuk mendapatkan nilai gaya Coulomb dimulai dengan input partikel yang dilakukan dengan cara klik pada citra lokasi piksel yang diinginkan. Koordinat piksel partikel disimpan pada variabel sx[s] untuk koordinat sumbu x dan sy[s] untuk koordinat sumbu y, dengan s adalah indeks partikel. Setelah proses inisialisasi, perhitungan dilakukan berdasarkan sumbu x dan y. Hasil dari perhitungan sumbu x disimpan dalam variabel CoulombX, dan CoulombY untuk menyimpan hasil dari perhitungan sumbu y.
C. Perancangan Perhitungan Pergerakan Partikel Gaya total didapatkan dengan menjumlahkan gaya Lorentz dan gaya Coulomb. Sesuai dengan rumus yang ada pada persamaan (2.10). Berdasarkan penelitian yang dilakukan Jalba dkk, partikel yang bergerak memiliki energi potensial sehingga disebut partikel dinamis, maka rumus untuk menghitung gaya total menjadi Persamaan (3.3). ( )=
( )+
( )−
............................................................. (3.3)
41
( )
=
( )
= −∇∅( )|
= ( ( ))
.................................................................. (3.4)
Dimana ∅ adalah energi potensial partikel, kecepatan, dan
adalah posisi,
adalah
adalah percepatan. Bila dibandingkan dengan Persamaan (2.3), ( )=−
Persamaan (3.3) memiliki tambahan,
, yaitu pengurangan
energi yang diperlukan partikel agar mencapai keadaan setimbang, digunakan untuk mengurangi energi potensial (Jalba dkk, 2004). Persamaan (3.3) dan (3.4) di atas belum bisa untuk disubstitusi dikarenakan
( ( )) belum diketahui,
namun masih bisa untuk mencari turunan pertamanya, yang artinya bisa menggunakan deret Taylor agar persamaan tersebut bisa digunakan. Untuk memudahkan, tiap partikel menggunakan interval waktu yang sama yaitu ∆ . Oleh karena itu, pada setiap interval waktu yang bekerja pada tiap partikel mendapatkan percepatan partikel
=
+ ∆ gaya
dihitung sesuai dengan persamaan (3.3) untuk . Kemudian untuk mendapatkan kecepatan,
percepatan, dan posisi dalam tiap interval waktu bisa menggunakan metode Euler. Metode ini dapat mengakibatkan error yang terakumulasi tiap interval waktu. Oleh karena itu, dikembangkan metode alternatif yang efisien yang didasarkan pada deret Taylor yang memanfaatkan nilai-nilai yang dihitung sebelumnya untuk mendapatkan percepatan dan kecepatan (yaitu metode bertingkat). Untuk mencari kecepatan, percepatan, dan posisi dalam satu interval waktu menggunakan persamaan seperti berikut: ( +∆ ) = ( )+∆
( )+ ∆
( +∆ )= ( )+∆
( )+ ∆
( )+ ∆ ( )
( )
................................... (3.5)
...................................................... (3.6)
42
Karena implementasi persamaan
( )/
pada komputasi susah
dilakukan, dapat digunakan (kuadrat) Lagrange polynomial ( ), dalam satuan waktu { ,
−∆ ,
− 2∆ } (Jalba dkk, 2004). Mengambil turunan dari
polynomial, dan mensubstitusi = , maka didapat persamaan sebagai berikut. ( )
=
( )
=
∆
(3 ( ) − 4 ( − ∆ ) + ( − 2∆ )) ................ (3.7)
Berdasarkan Persamaan (3.3) sampai (3.7) diatas maka bila diasumsikan massa dari partikel positif adalah 1 maka implementasi dalam program adalah sebagai berikut : percx[i][ite] = (w1 * CoulombX[i]) + (w2 * lorentx) - (beta * kecx[i][ite-1]); percy[i][ite] = (w1 * CoulombY[i]) + (w2 * lorenty) - (beta * kecy[i][ite-1]); kecx[i][ite] = kecx[i][ite-1] + (deltat * percx[i][ite-1]) + (0.5 * deltat * deltat * (((3 * percx[i][ite-1]) - (4 * percx[i][ite-2]) + percx[i][ite-3]) / (2 * deltat))); kecy[i][ite] = kecy[i][ite-1] + (deltat * percy[i][ite-1]) + (0.5 * deltat * deltat * (((3 * percy[i][ite-1]) - (4 * percy[i][ite-2]) + percy[i][ite-3]) / (2 * deltat))); posx[i][ite] = posx[i][ite-1] + (deltat * kecx[i][ite-1]) + (0.5 deltat * deltat * percx[i][ite-1]) + (0.167 * deltat * deltat deltat * (((3 * percx[i][ite-1]) - (4 * percx[i][ite-2]) + percx[i][ite-3]))); posy[i][ite] = posy[i][ite-1] + (deltat * kecy[i][ite-1]) + (0.5 deltat * deltat * percy[i][ite-1]) + (0.167 * deltat * deltat deltat * (((3 * percy[i][ite-1]) - (4 * percy[i][ite-2]) + percy[i][ite-3])));
* *
* *
