. Dengan perkembangan teknologi (khususnya komputer) batas antara fisika komputasi dan fisika teori menjadi semakin tidak nampak seiring dengan semakin tidak nampaknya perbedaan hasil yang diperoleh secara komputasi dan analitik. Ini merupakan salah satu alasan yang menyebabkan mengapa fisikawan teori dapat berperan sebagai fisikawan komputasi dan sebaliknya (contoh Feynman, Kenneth G Wilson, Steven Koonin dan sebagainya). Ramalan yang secara akurat dapat diperoleh dari hasil komputasi terhadap beberapa masalah fisika menunjukkan bahwa fisika komputasi tidak lagi hanya sekedar alat visualisasi atau simulasi proses fisis agar nampak sederhana. Penghargaan Nobel Fisika untuk K.G. Wilson dalam teori fenomena kritis (critical phenomena) dan Nobel Kimia untuk Pope dalam pengembangan teori fungsional kerapatan untuk molekul menunjukkan peranan fisika komputasi yang semakin mantap sebagai bentuk pendekatan ketiga (selain fisika teori dan fisika eksperimen) dalam mempelajari fisika. Fisikawan eksperimen perlu memiliki ketrampilan penanganan peralatan agar pengukuran besaran fisis yang dilakukannya dapat memberikan hasil. Pemahaman mengenai karakteristik alat yang dihadapi sangat membantu dalam melacak sumber-sumber kesalahan yang mungkin timbul atau bahkan dapat membantu melokalisir permasalahan. Fisikiwan teori perlu memiliki ketrampilan matematis agar dapat memperoleh penyelesaian terhadap persamaan matematis yang dihadapi. Disinipun diperlukan pemahaman yang mendalam tentang watak-watak fungsi dan operasi matematis agar mampu melalukan trik atau teknik tertentu untuk meyederhanakan persoalan matematis yang rumit. Hal yang semacam juga diperlukan dalam fisika komputasi. Fisikawan komputasi perlu memiliki ketrampilan untuk mengubah persamaan matematis atau hukum fisika ke dalam bentuk diskrit yang sesuai. Dalam hal ini diperlukan bentuk diskrit karena segala informasi fisis mengenai sistem nantinya akan tersaji dalam bentuk angka-angka atau nilai-nilai numerik. Metode pengubahan persamaan matematis ke bentuk diskrit beserta berbagai metode penyelesaiannya ini biasa disebut metode numerik. Pemahaman mengenai metode numerik itu saja belum mencukupi karena hasil nilai-nilai numerik yang diharapkan baru akan muncul setelah diproses oleh komputer. Dalam hal ini diperlukan bahasa komputer atau piranti tertentu (biasanya perangkat lunak berbentuk paket siap pakai) sehingga pemakai dapat berkomunikasi dengan komputer dan memerintahkannya untuk melakukan proses komputasi seperti yang diharapkan. Pemahaman tentang watak-watak alur langkah komputasi (algoritma) serta watak komputer itu sendiri tentunya sangat membantu dalam melakukan trik dan manipulasi untuk optimasi proses komputasi. Bidang fisika komputasi yang memerlukan tidak hanya pemahaman fisika tetapi juga metode numerik dan bahasa pemrograman menimbulkan kesan seolah-olah sebagai pengetahuan yang memerlukan multi disiplin ilmu. Untuk masa sekarang (paling tidak di Indonesia) kesan ini memang tidak salah. Seseorang yang belajar fisika di perguruan tinggi tidak otomatis mendapatkan kuliah metode atau analisis numerik dan bahasa pemograman. Keduanya biasanya diberikan di program studi lain seperti matematika atau ilmu komputer. Hal ini lain jika dibandingkan dengan kuliah fisika matematika (diperlukan dalam bidang fisika teori) dan praktikum/eksperimen fisika (diperlukan dalam bidang fisika eksperimen) yang tersedia dalam program studi fisika. Di masa depan keadaan ini mungkin akan berubah oleh beberapa sebab. Pertama karena komputer sudah memasyarakat dan user friendly sehingga interaksi dengan komputer baik melalui bahasa pemrograman atau piranti lain sudah dikuasai oleh setiap pengguna sejak dini. Kedua, metode numerik sendiri tidak terpisah lagi dengan fisika karena begitu intensif digunakan dalam banyak masalah fisika atau juga karena metode numerik itu sendiri dikembangkan melalui analogi proses fisika. 2. Beberapa aspek dasar Dalam bagian (1) telah disinggung bahwa inti dari proses komputasi adalah mengubah persamaan matematika ke dalam bentuk diskritnya sedemikian hingga perilaku sistem dan penyelesaiannya dapat diwakili oleh nilai-nilai numerik yang terlibat. *) Pemahaman persamaan dan arti fisisnya. Bagi pemula prinsip ini kadang digunakan secara tidak hati-hati : asalkan persamaan matematika yang menggambarkan sisten fisika sudah diubah ke bentuk diskrit maka nilai hasil keluaran yang diperoleh akan diterima apa adanya sebagai penyelesaian sistem tersebut. Perlu ditekankan bahwa meskipun fisika komputasi terlibat dengan angka-angka numerik, tapi bukan nilai angka-angka itu sendiri yang penting. Yang lebih penting adalah apakah nilai-nilai tersebut secara benar dapat menggambarkan sistem fisika yang ditinjau. Dengan kata lain fisika komputasi adalah bidang yang berhadapan dengan masalah fisika (seperti halnya fisika teori dan fisika eksperimen) dan bukan sekedar masalah angka/matematika. Pemahaman ini memberikan konsekuensi pada beberapa aspek dasar yang perlu diperhatikan, yang pada akhirnya akan berguna untuk mendapatkan informasi yang benar terhadap sistem yang ditinjau dan bahkan berguna dari sisi efisiensi dan ketelitian proses komputasi. Aspek yang perlu diperhatikan antara lain seperti yang diuraikan sebagai berikut. *) Satuan universal Berbagai besaran fisis di dalam fisika dapat dikaitkan dengan satuannya masing-masing. Mengetahui satuan yang digunakan untuk besaran fisis tertentu akan dapat diketahui keadaan sistem bagi besaran tersebut. Ketika besaran panjang dinyatakan dalam satuan angstrom dan besaran tenaga dalam satuan eV maka kira-kira sistem yang ditinjau adalah sistem mikroskopis non-relativistik. Sebaliknya ketika satuan panjang dinyatakan dalam meter dan tenaga dalam Joule maka yang terlintas adalah sistem makroskopis. Untuk penyelesaian analitis, pengambilan satuan sesuai dengan sistem fisisnya tersebut tidak begitu berarti pada perolehan hasil akhir karena selalu dapat dilakukan penyederhanaan pada langkah perhitungannya. Namun dalam komputasi pengambilan satuan yang sesuai tersebut sangat mempengaruhi hasil mengingat yang terlibat dalam komputasi hanyalah nilai-nilai numerik dan bukan simbol besaran fisis. Karena keterbatasan presisi komputer, suatu nilai yang sebenarnya tidak nol seperti nilai massa elektron m = 9.1 x 10-31 kg dapat diperlakukan sebagai nilai nol. Hal ini sangat berbahaya apabila muncul sebagai nilai pembagi pada pertengahan proses komputasi yang berakibat pada gagalnya proses komputasi. Cara lain untuk menjamin agar nilai yang terlibat dalam komputasi tidak terlalu besar atau terlalu kecil sesuai sistem fisisnya adalah dengan menormalisir persamaan matematik yang terlibat dengan cara mendefinisikan satuan universal bagi sistem sedemikian hingga semua besaran fisis yang terlibat menjadi tidak bersatuan. Keuntungan cara ini adalah selain persamaan yang terlibat menjadi berbentuk sederhana juga dimunkinkannya diperoleh ketelitian proses komputasi yang tinggi mengingat angka numerik yang terlibat berorde besar sesuai batas ketelitian komputer. Sumber diambil dari : Dr. Pekik Nurwantoro 21 Agustus 2000 BAB II. GALAT (KESALAHAN ) DALAM KOMPUTASI Perlu disadari bahwa komputasi numerik adalah komputasi yang mengikuti suatu algoritma pendekatan (aproksimasi) untuk menyelesaikan suatu persoalan. Dengan demikian bidang komputasi numerik memiliki kemungkinan kesalahan/galat sebagaimana terdapat di dunia eksperimen fisika. Beberapa sumber galat pada komputasi numerik adalah: 1. Round-off error : kesalahan akibat pembuatan angka 2. Truncation error : kesalahan akibat pemotongan suku pada deret aproksimasi, misalnya suatu rumus rumit diganti dengan rumus yang lebih sederhana. 3. Range error : kesalahan yang terjadi karena nilai hasil komputasi melampaui batas angka yang diperbolehkan oleh komputer, misalnya sangat kecil atau sangat besar. 1) Galat Pembulatan Pembulatan bilangan sering dilakukan dalam bidang komputasi numerik, misalnya mengurangi cacah digit pada suatu nilai hampiran dengan membuang/mengabaikan beberapa digit terakhir. Aturan pembulatan yang sering diterapkan adalah sebagai berikut: *) Bila digit yang dibulatkan kurang dari 5, digit depannya tidak berubah. *) Bila digit yang dibulatkan lebih dari atau sama dengan 5, digit depannya ditambah 1 nilainya. Contoh : dibulatkan menjadi 3.35 3.344 dibulatkan menjadi 3.34 Pengulangan pembulatan tidak disarankan dalam komputasi numerik karena akan memperbesar galat. Sebagai contoh, 3.3446 jika dibulatkan tiga angka dibelakang koma menjadi 3.345 dan jika dibulatkan lagi dua angka dibelakang koma menjadi 3.35. Galat pembulatan pertama sebesar 0.0004, galat pembulatan kedua sebesar 0.0054. Terlihat galat semakin besar sehingga pembulatan hanya boleh dilakkukan sekali. Kesalahan yang timbul akibat pembulatan pada digit ke-N di belakang koma, nilainya selalu 10-N / 2. Kesalahan pembulatan bisa terjadi karena komputasi terkait dengan angka-angka yang nilainya berbeda jauh, misalnya : hasil seharusnya 1, tetapi oleh komputer bisa jadi hasilnya karena pembilang ataupun penyebut nilainya kecil sekali sehingga dianggap nol. 2) Galat Pemotongan Galat pemotongan terjadi ketika suatu rumus komputasi disederhanakan. Biasanya suatu rumus rumit diganti dengan deret Taylor, lalu deret dipotong hanya sampai pada suku tertentu, suku selanjutnya diabaikan. Contoh : Misal x=1.5, nilai eksak cos (1.5) = 0.070737 sedangkan jika dihitung sampai suku ke-4 saja dihasilkan nilai komputasi cos(1.5) = 1 - (1.5)2/2! + (1.5)4/4! - (1.5)6/6! = 0.070187. Dengan demikian terdapat galat pemotongan sebesar 0.000550. 3) Galat batasan angka Setiap komputer memiliki keterbatasan dalam jangkauan representasi angka, misalnya angka presisi-tunggal (single precision) sekitar 10+37 dan presisi ganda (double precision) sekitar 10+308 . Melampaui batas jangkauan presisi tunggal dalam komputasi fisis dapat terjadi dengan mudah, misalnya dalam komputasi jari-jari atom Bohr, sebagai berikut: Walaupun hasil akhir tidak melampui batas single precision, tetapi pembilang dan penyebut sudah jelas melampaui sehingga komputasi bisa mengalami round-off dan hasil akhirnya tidak tepat. Masalah serupa juga dijumpai ketika dilakukan komputasi pada nilai yang sangat besar, misalnya komputasi faktorial : n! = n(n-1)(n-2)(n-3)…1 ketika n besar, misalnya n=200, aka nilai 200! Dapat melampaui besaran presisi ganda (double precision) . Jalan keluar untuk mengatasi masalah ini adalah dengan mencari rumus pendekatan. Sebagai contoh, faktorial dengan n besar dapat diganti dengan komputasi logaritmik, di mana : log(n!) = log (n) + log (n-1) + log (n-2) + log (n-3) + … Referensi : Suarga, 2007, Fisika Komputasi solusi problema fisika dengan MATLAB, Andi Yogyakarta Sahid, 2005, Pengantar Komputasi numerik, Andi Yogyakarta. BAB III. TURUNAN NUMERIK Pada bab ini akan dijelaskan metode numerik untuk menaksir nilai turunan suatu fungsi. Suatu fungsi f, baik diketahui rumusnya secara eksplisit maupun dalam bentuk data titik-titik yang dilalui kurvanya, dapat dihampiri dengan sebuah fungsi lain yang lebih sederhana. Suatu polinomial p merupakan pilihan yang paling mudah sebagai hampiran suatu fungsi f karena setiap polinomial dapat dengan mudah diturunkan. Turunan polinomial p, yakni p’(x) digunakan sebagai hampiran untuk f’(x) untuk sembarang nilai x. Secara geometris hal ini ekivalen dengan menghampiri gradien garis singgung pada kurva f di x dengan gradien garis singgung pada kurva p di x. Rumus-rumus (metode-metode) turunan numerik bermanfaat di dalam pengembanagn algoritma untuk menyelesaikan masalah nilai awal pada persamaan diferensial biasa dan parsial. 1. Turunan tingkat satu *) Rumus selisih maju dua titik (beda maju) Misalkan f(x) adalah sebuah fungsi riil satu variabel, deret Taylor untuk f(x) disekitar x = x0 adalah: (1) dengan c adalah sebuah bilangan antara x dan x0. Misalkan x = x0+h, dengan h>0, maka persamaan (1) dapat diltulis ulang sebagai: untuk suatu c [x0, x0+h] (2) Atau (3) Jika h mengecil, akan memberi taksiran untuk nilai f’(x0). Ini berarti : (4) disebut dengan rumus selisih maju dua titik (beda maju ) dengan galat sebesar O(h). Gambar 1 Perhatikan ruas kanan persamaan (4) adalah gradien tali busur yang melalui titik-titik (x0, f(x0)) dan (x0+h, f(x0+h)) sedangkan f’(x0) merupakan gradien garis singgung di titik (x0, f(x0)). Jadi gradien garis busur merupakan hampiran gradien garis singgung untuk nilai h yang cukup kecil sebagaimana ditunjukkan pada Gambar 1. Dengan simulasi dapat ditunjukkan bahwa gradien garis-garis bususr akan semakin menyamai gradien garis singgung jika h semakin kecil dan akhirnya gradien kedua garis sama apabila h 0. Agar lebih jelas lagi, perhatikan contoh berikut. Carilah turunan fungsi f(x)=ln x di x =1. Penyelesaian : Tabel berikut memberikan nilai-nilai turuan untuk beberapa nilai pilihan h. H ln(1+h) f’(1) 0.1 0.01 0.001 0.0001 0.00001 0.000001 0.0953101798043 0.00995033085317 0.000999500333083 0.0000999950000333 0.0000099999500004 0.0000009999994999 0.9531017980432 0.995033085317 0. 999500333083 0. 999950000333 0. 99999500004 0. 9999994999 Dari tabel terlihat bahwa semakin kecil nilai h yang digunakan (h mendekati nol), hampiran semakin mendekati nilai turunan yang sesungguhnya, yaitu = 1. 2. Rumus Selisih mundur dua titik (beda mundur) Rumus selisih maju menggunakan nilai fungsi di x0 dan x0 +h. Rumus serupa untuk mencari f’(x0) juga dapat diperoleh dengan menggunakan nilai-nilai fungsi di x0 dan x0-h, yakni : (5) disebut rumus selisih mundur dua titik dengan galat sebesar O(h). Visualisi selisih mundur dua titik untuk mendapatkan hampiran turunan fungsi di suatu titik diilustrasikan pada Gambar 2. Gambar2 3. Rumus Selisish pusat dua titik (beda pusat) Dari rumus selisih maju dua titik, diketahui : dan dari rumus selisih mundur dua titik, diketahui : Jika keduanya dijumlahkan diperoleh : Jadi, (6) disebut rumus selisih pusat dua titik dengan galat (O) h2 . Gambar3 Rumus selisih pusat dua titik menggunakan fakta bahwa gradien garis busur yang melalui titik-titik (x0-h, f(x0-h)) dan (x0+h, f(x0+h)) merupakan hampiran gradien garis singgung di x0 sebagaimana ditunjukkan pada Gambar 3. Dengan simulasi dapat ditunjukkan bahwa gradien garis busur semakin mendekati gradien garis singgung jika nilai h semakin kecil menuju nol. Rumus beda pusat memiliki galat yang lebih kecil dibandingkan rumus beda maju dan beda mundur. Ini dapat dibuktikan dengan contoh berikut : Contoh : Misalkan . Hitunglah hampiran f’(1) dengan menggunakan rumus beda pusat dan beda maju untuk nilai-nilai h =0.2 , 0.02 , 0.002, 0.0002. Penyelesaian : Untuk h=0.2 dengan menggunakan rumus beda maju diperoleh: dengan menggunakan rumus beda pusat diperoleh: Diketahui secara analitik , 2.71828182845905. Tabel berikut memuat hasil simulasi pencarian turunan dengan beda maju dan beda pusat untuk berbagai nilai h. h Beda maju Galat beda maju Beda pusat Galat beda pusat 0.2 0.02 0.002 0.0002 3.009175471388 2.745646775263 2.721001923382 2.718553674765 0.290893642928 0.027364946803 0.002720094922 0.000271846306 2.73643998561 2.71846305087 2.71828364064 2.71828184658 0.018158157151 0.000181222412 0.000001812188 0.000000018121 Dari tabel di atas tampak jelas bahwa galat metode beda pusat memberikan hasil yang lebih baik. BAB IV. PENCARIAN AKAR-AKAR FUNGSI NON LINIER, f(x) = 0 Banyak masalah-masalah fisika yang melibatkan fungsi/persamaan kontinu yang tidak linier. Persamaan non linier adalah persamaan di mana variabel di dalamnya pada umumnya tidak linier, misalnya dalam bentuk polinomial atau dalam bentuk perkalian atau pembagian beberapa variabel. Pada persamaan non linier dengan satu variabel, tujuan analisis pencarian akar –akar adalah mencari nilai variabel, mencari nilai x agar f(x)=0. Pada sistem persamaan non linier akan dijumpai lebih dari satu variabel yang terkait secara non linier dalam beberapa persamaan, kemudian akan dicari nilai dari masing-masing variabel yang membuat semua persamaan non linier bernilai nol. Terdapat banyak cara untuk menyelesaikan persamaan non linier seperti metode bagi dua (bisection), pisisi palsu (regula false) dan Newton Raphson. Pada bab ini akan dibahas metode Bagi dua dan metode Newton Raphson, karena dua metode ini merupakan metode fundamenal pada kasus pencarian akar-akar fungsi. non linier. *) Metode Bagi dua (bisection) Suatu fungsi kontinu pada interval tertutup [a b] sedemikian hingga f(a) dan f(b) berlawanan tanda, maka terdapat suatu akar persamaan f(x) = 0 pada interval [a b]. Akan dicari akar r (a,b) yang memenuhi f(r) = 0. Dari fakta ini muncullah metode pengapitan akar, berusaha mendapatkan interval kecil yang memuat suatu akar. Salah satu metode pengapitan akar adalah metode bagi dua. Metode bagi dua menggunakan konsep : Interval yang memuat akar dibagi menjadi dua subinterval sama panjang, kemudian dipilih subinterval mana yang memuat akar, dan selanjutnya sub interval ini dibagi dua lagi. Demikian seterusnya sampai diperoleh sebuah subinterval yang memuat akar yang dicari dan interval ini memiliki lebar tidak lebih dari nilai tertentu (sangat kecil mendekati nol). ALGORITMA BAGI DUA 1. Pilih dua titik a dan b sehingga f(a)*f(b)<0 2. Tetapkan toleransi error TOL dan maximum langkah iterasi N. 3. for i =1 sampaidengan N d = a+ (b-a)/2 cetak hasil (‘akar = “, d) if (f(p)=0) or ((b-a)/2
Sabtu, 22 Juni 2013
fisika komputasi
KEDUDUKAN FISIKA KOMPUTASI
1. Pendahuluan
Fisika Komputasi merupakan bidang yang mengkaji masalah fisika berdasarkan hasil tinjauan komputasi numerik. Sebagai perbandingan, fisika teori melandaskan bidang kajian fisika berdasarkan analisis matematis analitik sedangkan fisika eksperimen melandaskannya pada interpretasi hasil-hasil pengukuran beberapa besaran fisis yang terkait TOL
then a = d
else b=d
4. selesai
*) Metode Newton Raphson
Metode diperoleh dari ekspansi deret Taylor di sekitar x. Andaikan pada awalnya diberikan nilai x1 sehingga terjadi simpangan h = (x1 - x) dari akar yang dicari, nilai x dapat ditulis sebagai : x =x1-h.
Ekspansi deret Taylor di sekitar x memberikan :
karena nilai fungsi di titik akar x adalah f(x)=0, maka:
atau
dengan O(h2) adalah galat komputasi.
Karena h = x1-x maka
jadi .
Jika ditulis dalam bentuk iteratifnya :
Disebut dengan metode newton rapson untuk mencari akar-akar persamaan non-linier. Jika dipenuhi (telah konvergen, ) maka jelaslah bahwa akar x telah diperoleh sebesar xi+1.
ALGORITMA NEWTON RAPHSON
1. Tetapkan nilai awal x0
2. Tetapkan besar toleransi error TOL dan maximum langkah iterasi N.
3. for i =1: N
a. hitung f(x0) dan f’(x0)
b. hitung delta = f(x0) / f’(x0)
c. hitung x = x0 - delta
d. if ( delta < TOL)
tulis (‘akar yang di cari =’, x)
break
end if
else x0=x
e. if (i >= maximum langkah) tulis (‘gagal mencari akar’)
end if
Keunggulan metode Newton Raphson dibandingkan metode lain :
Hanya dibutuhkan satu nilai coba awal, sedangkan metode bagi dua membutuhkan dua nilai coba yang harus mengapit akar yang dicari.
Kelemahan : metode Newton Raphson membutuhkan informasi turunan pertama fungsi f’(x). Hal ini akan menjadi masalah jika fungsi terlalu kompleks sehingga turunan pertama fungsi tersebut sulit dicari.
Contoh kasus:
Sebuah bola pejal homogen terbuat dari suatu bahan dengan kerapatan seragam. Bola tersebut terbenam sebagian dalam air setinggi d. Berdasarkan hukum Archimedes : berat bola sebanding dengan gaya angkat Fa :
Massa bola :
Hukum Archimedes :
Volume bola yang terbenam di air setinggi d = Va, dapat dicari dengan kaedah volume benda putar (ingat kalkulus integral), sedemikian sehingga :
dimana ρ = ρb ;
Nilai d dapat dicari secara analitik dengan konsep pencarian akar-akar fungsi f(d) =0 ;
Masalah ini dapat juga diselesaikan secara komputasi . Nilai d dapat dicari dengan konsep pencarian akar-akar fungsi f(d) =0 dengan berbagai metode, diantaranya adalah metode bagi dua.
Jika diandaikan ρa = 1; ρb = ρ = 0.638; r =10 cm,
d = …?
Untuk mendapatkan nilai coba yang tepat, perlu dilihat bentuk fungsi nonlinier dari masalah ini. Hal ini dapat dibantu dengan melukis fungsi nonliniernya yakni :
dengan memasukkan ρa = 1; ρb = ρ = 0.638; r =10 cm, diketahui bentuk fungsi matematisnya seperti tampak pada gambar (tampak bahwa nilai d terletak di sekitar 10 < d < 20).
*) Berikut ini adalah metode bagi dua untuk menyelesaikan masalah ini ( mendapatkan nilai d) dengan bahasa MATLAB. Nilai coba pengapit akar adalah [a b] = [10 15].
clc;clear;
r=10; rho=0.638;
% nilai coba [a b]
a=10;
b=15;
% batas toleransi
tol=0.00001;
% banyaknya langkah
N=50;
hasil=[];
% kontrol loop
for i=1:N
d=a+(b-a)/2; % metode bagi dua
fd=pi*(d^3-3*d^2*r+4*r^3*rho)/3; % substitusi nilai d ke persamaan nonlinier
fa=pi*(a^3-3*a^2*r+4*r^3*rho)/3; % substitusi nilai a(batas kiri) ke persamaan nonlinier
hasil=[hasil;i a b d fd];
if (fd==0)||((b-a)=maxstep
fprintf('gagal mencapai akar hingga iterasi ke-%g', i);
disp('gagal dalam %g langkah', maxstep);
end
--------------------------------------------------------------------------------------------------------------------
Definisi fungsi ditulis dalam M-File terpisah bernama FNEWTRAPHSON(x)
--------------------------------------------------------------------------------------------------------------------
function [fx,f1x]=FNEWTRAPHSON(x,r,rho)
% definisi fungsi non-linier dan turunan pertamanya
r=10; rho=0.638;
fx=pi*(x^3-3*x^2*r+4*r^3*rho)/3;
f1x=pi*(3*x^2-6*x*r)/3;
return
-------------------------------------------------------------------------------------------------------------------
Hasil running program newton_raphson.m
taksiran awal x0=5
step x0 x delta
1 5.000000 13.564444 -8.564444
2 13.564444 11.761945 1.802499
3 11.761945 11.861318 -0.099373
4 11.861318 11.861502 -0.000184
5 11.861502 11.861502 -0.000000
akar x= 11.8615 pada iterasi ke-5
Berikut adalah penyelesaian masalah di atas dengan menggunakan bahasa Java.
import java.io.*;
public class akar {
public static final double PHI=3.14;
public static void main(String arg[])throws IOException{
int i,N;
double x0,x,dx,fx,ftx,tol,tol2,r,rho,delta,delta2;
// inisialisasi
x0=0.0;
// inputkan nilai coba
BufferedReader inputan = new BufferedReader(new InputStreamReader(System.in));
inputan = new BufferedReader(new InputStreamReader(System.in));
try{
System.out.print("masukkan nilai coba x0 =");
x0=Double.parseDouble(inputan.readLine());}
catch(Exception e){}
// batas iterasi dan toleransi yang diinginkan
N=100;
tol=0.00001;
tol2=tol*tol;
r=10.0;
rho=0.638;
dx=0.01;
// pencarian akar dimulai dengan loop WHILE
i=0;
while (iN){
System.out.println("gagal mencari nilai akar ");
}
}
}
public static float abs(float delta){
return delta;
}
public static double fgs_x(double x, double r,double rho) {
return (PHI*(x*x*x-3*x*x*r+4*r*r*r*rho)/3); // definisi fungsi
}
/* turunan fungsi diperoleh dari turunan numerik */
public static double fgs_dpn(double x, double r,double rho,double dx) {
return (PHI*((x+dx)*(x+dx)*(x+dx)-3*(x+dx)*(x+dx)*r+4*r*r*r*rho)/3);
}
public static double fgs_blk(double x, double r,double rho,double dx) {
return (PHI*((x-dx)*(x-dx)*(x-dx)-3*(x-dx)*(x-dx)*r+4*r*r*r*rho)/3);
}
}
===========================================================
Hasil running program newtonrapson.java
Kesimpulan :
Dari dua metode yang telah digunakan, terbukti bahwa metode Newton Raphson memberikan hasil yang lebih cepat dibanding metode bagi dua.
BAB V. INTEGRASI NUMERIK (satu dimensi)
Integrasi numerik disebut juga sebagai quadrature merupakan alat utama bagi para insinyur dan saintis untuk mendapatkan perkiraan jawaban suatu masalah fisis yang terkait dengan integral dan secara analitik sangat sulit diselesaikan. Quadrature satu dimensi pada prinsipnya memiliki konsep yang sederhana, yaitu bagaimana mengevaluasi integral suatu fungsi:
(1)
jika dipandang dari sudut persamaan diferensial aka mencari nilai integral I adalah sama dengan menyelesaikan persamaan diferensial dengan syarat batas y(a) = 0.
Metode yang umum digunakan dalam menghitung integral numerik adalah Newton-Cotes Formula, dimana batas antara a dan b dibagi ke dalam bagian yang lebih kecil (lebar langkah h) sedemikian rupa sehingga notasi integral dapat diganti dengan notasi penjumlahan(sigma), yaitu :
untuk loop tertutup
untuk loop terbuka
dimana fungsi f(x) adalah fungsi yang terintegralkan (kontinu) seperti Gambar berikut.
Ada dua cara dasar yang populer pada formula Newton-Cotes, yaitu Trapezoida-rule dan Simpson rule.
*) Aturan Trapezoid
Sesuai dengan namanya, integrasi numerik dengan aturan trapesium menggunakan cara menjumlahkan trapesium-trapesium kecil sebanyak N buah.
Hampiran penyelesaian integrasi numerik aturan trapesium dilakukan dengan langkah-langkah berikut:
1. Partisi interval [a, b] menjadi N subinterval berbentuk [xi, xi+1] sedemikian hingga .
2. Buat trapesium [xi, xi+1] dengan sisi-sisi yang sejajar f(xi), dan f(xi+1). Misalkan lebar subinterval adalah xi = xi+1 - xi. Maka luas trapesium yang terbentuk adalah (f(xi+ f(xi+1)) xi/2.
3. Jumlah luas trapesium-trapesium tersebut merupakan hampiran integral yang diinginkan, yakni:
(2)
Jika subinterval tersebut memiliki lebar sama, misalkan h, perhitungan di atas akan lebih mudah, dapt ditulis ulang sebagai:
(3)
dengan h =(b-a)/n, x0 = a, xn = b, xi+1 = xi +h , fi =f(xi)
*) Aturan Simpson
Selain metode trapesoid, metode fundamental pada integrasi numerik adalah menggunakan aturan Simpson. Aturan Simpson memiliki formula:
(4)
dengan n bilangan genap, h =(b-a)/n, x0 = a, xn = b, xi+1 = xi +h , fi =f(xi) dan
0 i n-1.
Contoh: Carilah nilai integral berikut :
Secara analitik hasilnya adalah :
Berikut ini diberikan dua implementasi atruran trapezoid dan aturan Simpson dengan bahasa MATLAB.
---------------------------------------------------------------------------------------------------------------------
% trapesium.m
% untuk mencari nilai integral
clc;clear;help trapesium;
a=1;
h=0.25;
b=9;
n=(b-a)/h;
s=0;
for i=1:n-1
x=a+i*h;
s=s+f(x);
end
I_trap=h*(f(a)+f(b))/2+h*s
---------------------------------------------------------------------------------------------------------------------
-----------------------------------------------------------------------
% simpson.m
% untuk mencari nilai integral
clc;clear;help simpson1;
a=1;
N=100;
b=9;
h=(b-a)/(2*N);
s1=0;
s2=0;
for i=1:N
x=a+h*(2*i-1);
s1=s1+f(x);
end
for i=1:N-1
x=a+h*2*i;
s2=s2+f(x);
end
I_simp=h*(f(a)+f(b)+4*s1+2*s2)/3
---------------------------------------------------------------------------------------------------------------------
Fungsi f(x) ditulis terpisah, yakni di Mfile f.m sebagai berikut:
------------------------------------------------------------------
function y = f(x)
y=1/x;
------------------------------------------------------------------
Jika file trapesium.m dikompile hasilnya :
trapesium.m
untuk mencari nilai integral
I_trap =
2.1972
Sedangkan file simpson.m dikompile hasilnya :
simpson.m
untuk mencari nilai integral
I_simp =
2.1972
Berikut ini contoh penyelesaian metode trapesium dan simpson dengan bahasa pemrograman Java.
================================================================
/*INTEGRASI NUMERIK DENGAN METODE TRAPESIUM */
public class trapesium {
public static void main(String args[]){
int i,N;
float a,h,b,s,x,fx,fa,fb,Trpsm;
N=100; // banyaknya segmen integrasi
a=1; // batas bawah integrasi
b=9; // batas atas integrasi
h=(b-a)/N; // lebar segmen
s=0; // inisiasi penjumlahan
fa=fs_x(a); fb=fs_x(b); // nilai fungsi pada batas
// hitung nilai fungsi tiap segmen
for(i=1;i<=N-1;i++){
x=a+i*h;
fx=fs_x(x); // panggil nama fungsi
s=s+fx; }
Trpsm=h*(fa+fb)/2+s*h; // RUMUS TRAPESIUM
System.out.println(" hasil integral f(x) dengan metode trapesium = "+Trpsm);
}
// definisi fungsi
public static float fs_x(float x){
return 1/x; // definisi fungsi yg hendak diintegralkan
}
}
================================================================
Hasil running program trapesium.java
run:
hasil integral f(x) dengan metode trapesium = 2.197751
BUILD SUCCESSFUL (total time: 0 seconds)
================================================================
/*INTEGRASI NUMERIK DENGAN METODE SIMPSON */
public class simpson {
public static void main(String arg[]){
int i,N;
float a,h,b,s,s1,s2,x,fx,fa,fb,spm;
N=100; // banyaknya segmen integrasi
a=1; // batas bawah integrasi
b=9; // batas atas integrasi
h=(b-a)/(2*N); // lebar segmen
s1=0; // inisiasi penjumlahan bagian ganjil
s2=0; // inisiasi penjumlahan bagian genap
fa=fs_x(a); fb=fs_x(b); // nilai fungsi batas bawah dan batas atas
/* hitung nilai-nilai bagian ganjil */
for(i=1;i<=N;i++){
x=a+h*(2*i-1);
fx=fs_x(x);
s1=s1+fx;
}
/* hitung nilai-nilai bagian genap */
for (i=1;i<=N-1;i++){
x=a+h*2*i;
fx=fs_x(x);
s2=s2+fx;
}
spm=h*(fa+fb+4*s1+2*s2)/3; // rumus simpson
System.out.println(" hasil integral 3x dengan metode simpson = "+spm);
}
/* definisi fungsi */
public static float fs_x(float x){
return 1/x;
}
================================================================
Hasil running program simpson.java
run:
hasil integral 3x dengan metode simpson = 2.1972246
BUILD SUCCESSFUL (total time: 1 second)
BAB VI. PERSAMAAN DIFERENSIAL BIASA
Sebagian besar masalah fisis (fisika) menyangkut kajian suatu sistem selama periode waktu tertentu. Hal ini berarti masalah yang ditinjau menggunakan bentuk persamaan diferensial dengan waktu sebagai variabel bebas. Bagi orang matematika, kajian persamaan diferensial merupakan salah satu bagian tercantik, karena kajian ini dapat diaplikasikan di fisika, kimia, teknik, biologi dan ekonomi. Sebagai contohnya: masalah-masalah sistem pegas, rangkaian induktansi resistor-kapasitor, pemuaian lempeng logam , reaksi kimia, ayunan pendulum, gerak putar massa mengitari benda lain, pertumbuhan bakteri, dan lain-lain dapat dimodelkan dalam bentuk persamaan diferensial.
Pada bab ini akan dibahas beberapa meto de fundamental untuk menyelesaikan persamaan diferensial biasa (Ordinary differential equations/ ODE) yaitu metode Euler dan metode Runge Kutta 4. Dua metode tersebut nantinya akan diterapkan pada contoh penyelesaikan masalah gerak jatuh bebas yang menyertakan gaya hambatan udara.
Suatu masalah yang terkait dengan persamaan diferensial orde satu dapat ditulis sebagai:
(1)
dimana f(x,y) adalah fungsi yang diberikan oleh masalah yang bersangkutan. Selanjutnya hendak dicari bentuk eksplisit dari y(x). Untuk mendapatkan bentuk ekslisit y(x) dibutuhkan informasi awal y(0) sehingga kajian ini juga disebut sebagai masalah syarat awal. Dari informasi awal y(0) dan f(x,y) selanjutnya dapat dicari nilai y(x) secara komputasi dengan metode Euler dan metode Runge-Kutta4.
*) Metode Euler
Konsep metode Euler berasal dari pemotongan deret Taylor hingga suku orde pertama :
(2)
Jika syarat awal f(x) diketahui dan turunan pertama fungsi f’(x) juga diketahui maka nilai hampiran f(x+x) berikutnya dapat dicari.
Dengan metode Euler, persamaan (2) dapat dirumuskan sebagai:
i = 0,1,2,... (3)
dengan , dan lebar langkah . Dari informasi awal y(0) dan f(x,y) selanjutnya dapat dicari nilai y(x). Ini berarti nilai y(x) sudah berhasil diperoleh untuk sebarang nilai x yang kita tentukan.
Metode Euler disebut juga metode satu langkah karena informasi satu langkah sebelumnya digunakan untuk menghitung hampiran sekarang. Dilihat dari pemotongan deret Taylor tampak bahwa galat metode Euler sebesar O(h).
*) Metode Rungge Kutta4
Oleh karena metode Euler memiliki galat sebesar O(h) maka diperlukan metode lain yang memberikan galat yang lebih baik. Pada bagian ini akan ditunjukkan metode Rungge Kutta4. Namun sebelumnya akan ditinjau sekilas metode RungeKutta2 (metode Heun).
Ide metode Runge Kutta adalah menghitung nilai f(x,y) pada beberapa titik di dekat kurva penyelesaian yang dipilih dengan metode tertentu di dalam interval (xi , xi +h) dan mengkombinasikann nilai-nilai ini sedemikian sehingga diperoleh keakuratan yang lebih baik pada hampiran berikutnya, yi+1.
Pada metode RungeKutta2 nilai yi+1 dihitung dengan langkah :
(4)
untuk i = 0, 1, 2,3,...dan syarat awal (x0, y0) diketahui.
Jika metode RungeKutta2 menggunakan dua langkah perhitungan, maka metode RungeKutta4 menggunakan 4langkah perhitungan untuk memperoleh hampiran nilai selanjutnya, yaitu :
(5)
*) Contoh kasus fisika yang dapat diselesaikan dengan metode Euler dan RunggeKutta4:
Pada mata kuliah mekanika, telah ditinjau kasus gerak jatuh bebas yang menyertakan pengaruh gaya hambatan udara. Penyelesaian analitik yang telah ada nantinya dapat digunakan sebagai pembanding penyelesaian secara komputasi. Dua metode komputasi (metode Euler dan metode Rungge Kutta4) akan dipakai untuk mencari solusi masalah tersebut.
Pada kondisi riil, suatu partikel bermassa m yang jatuh bebas di udara akan mengalami gaya tarik gravitasi sehinga bergerak dipercepat menuju tanah. Jika hambatan udara cukup besar (misalnya pada kasus penerjun payung) maka pada suatu saat kecepatan partikel tersebut akan konstan (gaya hambat udara seimbang dengan gaya berat). Pada kasus ini gaya hambat tersebut biasanya dinyatakan dengan suatu fungsi linier yang gayut kecepatan f(v) sebagai berikut :
F = mg – f(v) (6)
F = mg – bv (7)
Dimana b adalah koefisien hambatan.
Dari persamaan (7) diperoleh :
(8)
(9)
Persamaan diferensial orde dua tersebut (9) dapat dipecah menjadi dua persamaan differensial orde satu dalam bentuk :
(10)
(11)
Dengan syarat awal v = 0 saat t = 0, penyelesaian analitik untuk kecepatan v dari persamaan (10) di atas, adalah :
(12)
Penyelesaian analitik untuk kasus ini nanti akan digunakan sebagai pembanding metode komputasi (Euler dan Rungke kutta4).
Pada kasus ini hanya akan dicari penyelesaian kecepatan gerak saja, jadi hanya akan digunakan persamaan (10), sementara penyelesaian posisi gerak (persamaan (11)) tidak disinggung.
a) Metode Euler :
Persamaan (10) dapat dinyatakan sebagai :
dimana g dan c = b/m adalah suatu konstanta.
Menggunakan metode Euler,
Dengan syarat awal v(0) = 0 ; ukuran langkah h = ti+1 – ti; dimana n=0,1,2,…,N akan dicari vi = v(ti).
Misal: ingin dicari kecepatan gerak selama 5 sekon dengan jumlah langkah komputasi N=20. Jika massa benda, m =1 kg, konstanta gaya hambat, b= 1, percepatan gravitasi bumi, g = 9,8 m/s2. Dengan diketahi syarat awal v(0) = 0, implementasi kasus ini dengan MATLAB ditulis sebagai berikut :
-----------------------------------------------------------------------------------------------------------
%komp_uin.m
clc; clear; help komp_uin;
% input
m=1; b=1; g=9.8; gamma=b/m;
v0=0; N=20; t0=0; tn=5;
h=(tn-t0)/N;
for i=1:N
t(i)=t0; v(i)=v0;
t(i+1)=t(i)+h;
v(i+1)=v(i)+h*(g-gamma*v(i));
%set menjadi nilai awal
t0=t(i+1); v0=v(i+1);
%analitik
v_anltk(i+1)=(m/b)*g*(1-exp(-gamma*t(i+1)));
end
disp('[t v v_anltk]')
disp([t' v' v_anltk'])
plot(t,v,'o',t,v_anltk,'*')
xlabel('waktu t');ylabel('kecepatan');
title('kecepatan gerak jatuh bebas dengan gaya hambat')
------------------------------------------------------------------------
Hasil running program :
[t v v_anltk]
0 0 0
0.2500 2.4500 2.1678
0.5000 4.2875 3.8560
0.7500 5.6656 5.1708
1.0000 6.6992 6.1948
1.2500 7.4744 6.9923
1.5000 8.0558 7.6133
1.7500 8.4919 8.0970
2.0000 8.8189 8.4737
2.2500 9.0642 8.7671
2.5000 9.2481 8.9956
2.7500 9.3861 9.1735
3.0000 9.4896 9.3121
3.2500 9.5672 9.4200
3.5000 9.6254 9.5041
3.7500 9.6690 9.5695
4.0000 9.7018 9.6205
4.2500 9.7263 9.6602
4.5000 9.7448 9.6911
4.7500 9.7586 9.7152
5.0000 9.7689 9.7340
b) Metode Rungge Kutta4
Dari persamaan (10) diketahui turunan orde pertama :
dimana g dan c = b/m adalah suatu konstanta.
Metode Runge-Kutta4 untuk mencai nilai kecepatan digunakan langkah-langkah:
v(i)=0;
%komp_rk4.m
clc; clear; help komp_uin_rk4;
% input
m=1; b=1;
g=9.8; gamma=b/m;
v0=0;
N=20;
t0=0;
tn=5;
h=(tn-t0)/N;
for i=1:N
t(i)=t0; v(i)=v0;
% kecepatan
t(i+1)=t(i)+h;
v_1(i)=g-gamma*v(i);
v_2(i)=g-gamma*(v(i)+(h/2)*v_1(i));
v_3(i)=g-gamma*(v(i)+(h/2)*v_2(i));
v_4(i)=g-gamma*(v(i)+h*v_3(i));
v(i+1)=v(i)+(h/6)*(v_1(i)+2*v_2(i)+2*v_3(i)+v_4(i));
%set menjadi nilai awal
t0=t(i+1);
v0=v(i+1);
%analitik
v_anltk(i+1)=(m/b)*g*(1-exp(-gamma*t(i+1)));
end
disp('[t v v_anltk]')
disp([t' v' v_anltk'])
plot(t,v,'o',t,v_anltk,'*')
xlabel('waktu t');ylabel('kecepatan');
title('kecepatan gerak jatuh bebas dengan gaya hambat')
Hasil running program :
komp_rk4.m
[t v v_anltk]
0 0 0
0.2500 2.1677 2.1678
0.5000 3.8559 3.8560
0.7500 5.1707 5.1708
1.0000 6.1946 6.1948
1.2500 6.9921 6.9923
1.5000 7.6132 7.6133
1.7500 8.0969 8.0970
2.0000 8.4736 8.4737
2.2500 8.7670 8.7671
2.5000 8.9955 8.9956
2.7500 9.1734 9.1735
3.0000 9.3120 9.3121
3.2500 9.4200 9.4200
3.5000 9.5040 9.5041
3.7500 9.5695 9.5695
4.0000 9.6205 9.6205
4.2500 9.6602 9.6602
4.5000 9.6911 9.6911
4.7500 9.7152 9.7152
5.0000 9.7340 9.7340
Dari dua metode di atas, tampak bahwa metode Runge-Kutta4 hasilnya jauh lebih baik dan sesuai dengan hasil analitik (paling tidak untuk ketelitian empat angka di belakang koma).
Langganan:
Posting Komentar (Atom)
Tidak ada komentar:
Posting Komentar