Fisika Komputasi
I DERET TAYLOR 1.1 Deret Taylor Deret Taylor memegang peranan yang sangat penting dalam analisis numerik. Dengan deret Taylor kita dapat menentukan nilai suatu fungsi di titik x jika nilai fungsi di titik x0 yang berdekatan dengan titik x diketahui. Uraian deret Taylor disekitar xo dinyatakan dengan f(x) yaitu :
f ( x)
f ( n ) ( x0 ) ( x x0 ) n …………………………………………..(1.1) n!
f ( n ) ( x0 )
dn f dx n
dimana, x x0
…………………………………………………...(1.2)
f (0) ( x0 ) f ( x0 ) ………………………………………………………(1.3)
1.2 Deret Mac Laurin Deret Mac Laurin adalah bentuk khusus dari deret Taylor. Dimana dianggap bahwa titik x0 = 0 sehingga persamaan 1.1 akan berubah menjadi :
f ( x)
f ( n ) (0) n x …………………………………………………..(1.4) n!
dimana,
f ( n ) (0)
dn f dx n
x 0
……………………………………………………...(1.5)
f ( 0) (0) f (0) ………………………………………………………….(1.6) Persamaan (1.4) adalah persamaan dari deret Mac Laurin. Persamaan (1.1) biasa dituliskan dengan mensubstitusikan x dengan x-x0, sehingga :
f ( x x0 )
f ( n ) ( x0 ) n ( x) …………………………………………..(1.7) n!
Sehingga Persamaan deret Mac Laurin akan menjadi :
f ( x)
f ( n ) (0) n ( x) …………………………………………...........(1.8) n!
Contoh : Hitunglah deret Mac Laurin untuk Sin(x) © Avid-06
1
Fisika Komputasi
Jawab : f(x)
= sin(x)
;
f(0)
= 0
f1(x) = cos(x)
;
f1(0)
= 1
11
11
f (x) = -sin(x)
;
f (0) = 0
f111(x) = -cos(x)
;
f111(0) = -1
f1V(x) = sin(x)
;
f1V(0) = 0
fV(x) = cos(x)
;
fV(0)
= 1
Dari deret Mac Laurin :
f ( x)
f ( n ) (0) n ( x) n!
f ( 0) (0) x 0 f (1) (0) x1 f (11) (0) x 2 f (111) (0) x 3 f (1V ) (0) x 4 f (V ) (0) x 5 ... 0! 1! 2! 3! 4! 5!
Sehingga Sin(x) menjadi :
0.1 1.x 0.x 2 1.x 3 0.x 4 1.x 5 ... 0! 1! 2! 3! 4! 5!
x
x3 x5 x7 ... 3! 5! 7!
Latihan : 1. Buat deret Mac Laurin untuk fungsi cos(x) dan tan(x) 2. Buat deret Taylor untuk ex di sekitar x=0
1.3 Program komputer Program komputer dari fungsi sin(x) untuk deret Mac Laurin tersebut jika dinyatakan dalam bahasa Visual C++ adalah : //********************************************************* // Menghitung fungsi sinus dengan deret MacLaurin // Compiler : Visual C++ //********************************************************* #include
#include <math.h> void main() { double eps=1e-5; double x,sinx,pem,pen,s,del; int i,j,m,k,tanda;
© Avid-06
2
Fisika Komputasi
cout << "Sudut x (derajat) = "; cin >> x; x=x/57.3; cout << "sin(x) = " << sin ( x); //nilai sinus dari fungsi pustaka Visual C++ del=10; cout <<endl<<endl; cout << "Menurut Uraian deret Taylor : "<<endl; m=0; while (del > eps) { m=m+1; s=0; j=0; for (i=1;i<=2*m-1;i++) { if (i%2 ==1 ) // atau if ( fmod(i,2) == 1) { j=j+1; if (j % 2 == 0) tanda = -1; else tanda=1; pem=1; pen=1; for (k=1;k<=i;k++) { pem=pem*x; pen=pen*k; } s=s+tanda*pem/pen; } }
if (m>1 ) del= fabs(sinx-s); sinx=s; } cout << "Jumlah suku = "<<m <<endl; cout << "sin(x)
= " << s <<endl;
}
© Avid-06
3
Fisika Komputasi
Latihan : 1. Buatlah program komputer dari deret MacLaurin untuk fungsi exp(x) :
exp( x) 1 x
x 2 x3 ... 2! 3!
2. Buatlah program komputer dari deret MacLaurin untuk fungsi cos(x) :
cos( x) 1
© Avid-06
x2 x4 x6 ... 2! 4! 6!
4
Fisika Komputasi
II AKAR-AKAR PERSAMAAN Akar suatu persamaan merupakan titik potong fungsi f(x) dengan sumbu x :
f(x)
x0
x akar f(x)
Gambar 2.1. Akar dari persamaan f(x) Akar dari suatu persamaan sering dipakai untuk menyelesaikan masalah fisika seperti : -Mengetahui volume gas pada suatu tekanan tertentu P
saat P=P1 berapa nilai V1
P
1
V 1
V
Gambar 2.2. Diagram P-V gas ideal -Mengetahui jarak jatuh peluru yang ditembakkan dengan arah dan kecepatan tertentu f(x) f(x) saat f(x)=0 berapa nilai x
V(t)
x=?
x
Gambar 2.3. Gerak Peluru 2.1 Memperkirakan Letak Akar Jika diberikan fungsi f(x) dengan letak akar di x0 seperti di gambar 2.4. Untuk menentukan letak dari akar fungsi f(x) dapat ditempuh langkah - langkah sbb : © Avid-06
5
Fisika Komputasi
f(x)
x
x0
Gambar 2.4. Akar dari persamaan f(x) -tentukan sebuah titik pada sumbu x dalam daerah berlakunya fungsi f(x) misalnya x1 -tentukan titik lain misal x2, dimana x2 = x1- h -bila f(x1) x f(x2) > 0 , berarti x1 dan x2 berada pada pihak yang sama terhadap x0. -tentukan titik lain misal : x3 = x2-h -bila f(x2) x f(x3) > 0 , berarti x2 dan x3 berada pada pihak yang sama terhadap x0. f(x) x1 dan x2 sama-sama berada disebelah kanan x0
f(x2)
f(x1)
h
x0
x x2
x1
Gambar 2.5. Perkiraan akar dari persamaan f(x) -penentuan titik-titik diteruskan sampai diperoleh : f(xa) x f(xb) < 0 (negatif) f(x)
x0 xa
xa dan xb terletak berseberangan terhadap titik potong x0 xa berada di sebelah kiri x0 xb berada di sebelah kanan x0
x xb
Gambar 2.6. Perkiraan akar dari persamaan f(x) -sehingga daerah pencarian akar-akar akan dipersempit diantara xa dan xb -langkah-langkah diatas diulang pada daerah : xa x xb © Avid-06
6
Fisika Komputasi
catatan : Bila kita menentukan suatu arah dan sampai batas daerah tertentu tidak ditemukan titik potong (x0), maka kita coba sekali lagi untuk arah yang berlawanan Contoh : Tentukan titik potong f(x) terhadap sumbu x. Dimana bentuk persamaan f(x) : f(x) = x3-11x2 + 39x – 45 untuk 0 x 6,
h 0.75
n
xn
f(xn)
xn-h
1
6
9
5,25
2
5,25
1,265
4,5
3
4,5
-1,25
Dari tabel f(5,25) x f(4,5) < 0 , sehingga titik potong terletak antara x = 4,5 dan x = 5,25 2.2 Metode Newton -Rapshon Penentuan akar-akar suatu persamaan dapat dicari dengan memakai rumus iterasi Newton Rapshon yang diturunkan dari uraian deret Taylor. Jika hanya dilibatkan sampai suku dengan turunan pertama maka akan didapatkan iterasi Newton-Rapshon orde-1,dan bila yang dilibatkan sampai suku dengan turunan kedua maka disebut iterasi NewtonRapshon orde-2, dan seterusnya. Uraian deret Taylor di sekitar x0 dinyatakan :
f ( x)
f ( n ) ( x0 ) ( x x0 ) n n!
f ( 0) ( x0 ) x 0 f (1) ( x0 )( x x0 )1 f (11) ( x0 )( x x0 ) 2 ... 0! 1! 2! orde -1 orde -2
2.2.1 Metode Newton -Rapshon orde-1 Uraian deret Taylor dari fungsi f(x) di sekitar x0 sampai suku dengan turunan pertama dinyatakan : f(x) = f(x0) + f1(x0)(x-x0) Jika x disubstitusikan dengan xn+1 dan x0 disubstitusikan dengan xn maka : f(xn+1) = f(xn) + f1(xn)(xn+1-xn)………………………...……………………….(2.1) Bila xn+1 adalah akar dari f(x) maka f(xn+1) = 0. Sehingga persamaan 2.1 akan menjadi :
© Avid-06
7
Fisika Komputasi
x n 1 x n
f ( xn ) ……………………………………………………………..(2.2) f ' ( xn )
Persamaan 2.2 disebut dengan persamaan iterasi Newton-Rapshon orde-1. Nilai xn+1 dapat dicari dari xn (suatu titik dapat dicari dari titik sebelumnya). Jadi bila suatu titik awal x0 diketahui maka dengan persamaan 2.2 dapat dicari x1, x2, dan seterusnya sampai ditemukan akar dari persamaan atau fungsi yang bersangkutan (nilai x yang menyebabkan f(x) = 0). Penentuan x0 memegang peranan yang amat penting karena cepat lambatnya penentuan akar ditentukan oleh pengambilan titik awal x0. Penentuan titik awal x0 harus ditentukan sedekat mungkin dengan akar sehingga mengurangi banyaknya iterasi. Hal ini dapat dilakukan dengan metode memperkirakan letak akar (sub bab 2.1). Contoh : Tentukanlah akar f(x) = x3 – 11x2 + 39x – 45 dengan metode Newton-Rapshon orde-1. Jawab : f1(x) = 3x2 –22x + 33 misal kita ambil titik awal x0 = 6 n
xn
f(xn)
f1(xn)
xn+1
0
6,000
9,000
15,000
5,400
1
5,400
2,304
7,680
5,100
2
5,100
0,441
4,830
5,009
3
5,009
0,036
4,070
5,000
4
5,000
0,000
4,000
5,000
Dari tabel diatas dapat dilihat untuk ketelitian 3 angka dibelakang koma x 4 = x5 = 5,000. Dan nilai f(x) = 0,000. Jadi akar dari f(x) adalah x = 5,000 2.2.2 Metode Newton -Rapshon orde-2 Jika kita melibatkan sampai suku dengan turunan kedua pada uraian deret Taylor maka akan menghasilkan iterasi Newton Rapshon orde-2. Penurunan persamaannya dapat ditentukan dengan cara berikut ini. Dari persamaan umum deret Taylor (pers. 1.1)
f ( x)
f ( n ) ( x0 ) ( x x0 ) n n!
f ( 0) ( x0 ) x 0 f (1) ( x0 )( x x0 )1 f (11) ( x0 )( x x0 ) 2 0! 1! 2! dengan mensubstitusi x dengan xn+1, dan x0 dengan xn didapatkan © Avid-06
8
Fisika Komputasi
f 11 ( xn )( xn 1 xn ) 2 f ( xn1 ) f ( xn ) f ( xn )( xn1 xn ) 2! 1
Jika xn+1 akar dari f(x) maka f(xn+1) = 0
x n 1 xn
f ( xn ) f 1 ( xn )
Sehingga
f 11 ( xn )( xn 1 xn ) 2 0 f ( xn ) f ( xn )( xn 1 xn ) 2 1
f ( xn ) 1 f ( xn )
f 11 ( xn ) f ( x n ) 0 f ( xn ) f ( xn )( xn 1 xn ) ( xn1 xn ) 2 f 1 ( xn ) 1
1 f ( xn ) f 11 ( xn ) 0 f ( xn ) ( xn 1 xn ) f ( xn ) 2 f 1 ( xn ) x n 1 x n
f ( xn ) f ( x n ) f '' ( x n ) ' f ( xn ) 2 f ' ( xn )
x n 1 x n
f ( xn ) …………………………………………….(2.3) f ( x n ) f '' ( x n ) ' f ( xn ) 2 f ' ( xn )
Persamaan 2.3 merupakan persamaan iterasi Newton-Rapshon orde-2. Contoh : 1. Tentukanlah akar f(x) = x3 – 11x2 + 39x – 45 dengan metode Newton-Rapshon orde-2. Jawab : f1(x) = 3x2 –22x + 33 ; f11(x) = 6x - 22 misal kita ambil titik awal x0 = 6 N
xn
0
6,000
1
5,167
F(xn) 9,000
f1(xn) 15,000
f11(xn) 14,000
xn+1 5,167
2 3 4
© Avid-06
9
Fisika Komputasi
2. Satu mole gas ideal suhunya 3000 K, volumenya diperbesar pada suhu tetap (isotermis). Berapa liter volume gas ideal tersebut pada saat tekanannya 15 atmosfer. (tetapan gas universal 0,0823 l atm/mole 0K ) jawab: Persamaan gas ideal : PV = n R T Diket : N = 1 mole R = 0,0823 l atm / mole 0 K T = 300 0 K P = n R T / V = 1. 0,0823 . 300/ V = 24,69 / V Jika dimisalkan P = f(x) dan V = x maka f(x) = 24,69/x
P
f(x)
P1=15
f(x)1=15
V
x 1
x =?
V1=?
Karena yang dicari adalah nilai V saat P = 15 maka
P
24,69 24,69 15 atau f ( x) 15 V x
sehingga f 1 ( x)
24,69 x2
Dengan iterasi Newton-Rapshon orde-1
x n 1 x n
f ( xn ) f ' ( xn )
misal x0 = 1,5 n 0
xn
f(xn)
f1(xn)
xn+1
1,500
1 2
© Avid-06
10
Fisika Komputasi
3. Suatu celah lebarnya 0,4 mm disinari dengan cahaya dengan panjang gelombang 5900 Å. Pola difraksi yang terjadi ditangkap layar yang jaraknya 70 cm dari celah. Berapa jarak dari tengah-tengah terang pusat sampai intensitas cahaya tinggal ½ nya. Jawab : Pola difraksi pada celah ditentukan dengan :
I
I 0 Sin 2
2
dimana, I0 = intensitas cahaya yang datang β = ½ k b sin θ k = bilangan gelombang = 2π/λ ; sin θ = δ / D b = lebar celah θ = sudut antara normal pada celah dengan arah sinar yang menjadi pusat perhatian Jika I/I0 disubstitusi menjadi y dan β disubstitusi menjadi x maka :
y
sin 2 x x2
karena yang dicari adalah nilai x saat y = ½ (I / I0 = ½ atau I = ½ I0 ) Jadi persamaan yang dicari solusinya adalah :
y
sin 2 x 0,5 x2 f(x)
x1
y1=1/2
0 x=?
x
Untuk mendapatkan titik potong f(x) dengan sumbu x maka dipakai iterasi NewtonRapshon orde-1.
© Avid-06
11
Fisika Komputasi
f ( x)
sin 2 x 0,5 x2
2 sin 2 x sin 2 x f ( x) x3 x2 1
misal x0 = 1,25 n 0
xn
f1(xn)
f(xn)
xn+1
1,250
1 2
Dari perhitungan didapatkan titik potong x = 1,39155738 sehingga
1,392
1 kb sin b sin 2
5,9.10 4 sin 1,392. .0,4 2.3 Akar – Akar Dari Polinom Bentuk Umum Polinom adalah : f ( x) an x n an1 x n1 an2 x n2 ... a1 a0 ………………………..(2.4)
contoh :
f ( x) x 3 11x 2 33x 45
polinom orde 3
f ( x) x 2 2 x 8
polinom orde - 2
2.3.1 Pembagian Sintetik adalah metode untuk mencari nilai fungsi f(x) dan nilai turunan fungsi f1(x). Adapun langkah-langkahnya adalah : 1.Tuliskan koefisien-koefisien polinom dengan urutan :
an an1 an2 ... a0 bila ada suku yang tidak hadir tulis 0 untuk koefisien dari suku tsb. 2.Nilai x dimana nilai fungsi akan dicari diletakkan paling kiri dalam deretan koefisienkoefisien.
x an an1 an2 ... a0
© Avid-06
12
Fisika Komputasi
3.Selang satu baris dari deretan koefisien tersebut buat garis lurus sebagai garis penjumlahan.
x an an1 an2 ... a0 + 4.Turunkan an hingga ke bawah garis penjumlahan. Sebut an dengan bn.
x an an1 an2 ... a0 +
bn
5.Kalikan x dengan bn dan letakkan hasilnya diatas garis penjumlahan tepat di bawah an-1. Jumlahkan an-1 dengan hasil kali tadi dan letakkan hasilnya di bawah garis penjumlahan, satu kolom dengan an-1. Beri nama hasil penjumlahan tadi dengan bn-1
x
an
bn
an1 an2 ... a0 x.bn a n 1 ( x.bn )
+
bn 1 6.Kalikan x dengan bn-1 dan letakkan hasilnya diatas garis penjumlahan di bawah an-2. Jumlahkan an-2 dengan hasil kali tadi dan letakkan hasilnya di bawah garis penjumlahan, satu kolom dengan an-2. Beri nama hasil penjumlahan tadi dengan bn-2
x
an
bn
an1
an 2 ... x.bn 1
a0
bn1
a n 2 ( x.bn 1 )
+
bn 2 7.Lakukan operasi yang sama untuk koefisien-koefisien yang lain sehingga di bawah garis penjumlahan akan didapatkan deretan :
bn bn1 bn2 ... b0 8.Nilai b0 merupakan nilai fungsi f(x) yang dicari. Jika pembagian sintetik dilanjutkan untuk deretan
bn bn1 bn2 ... b0 maka akan didapatkan nilai dari turunan f1(x). Jika langkah-langkah tersebut diatas dirangkum maka akan tampak seperti berikut ini :
x
an
an1
an2
a n 3
(bn x) (bn1 x) (bn2 x) bn © Avid-06
bn1
bn2
bn3
+ 13
Fisika Komputasi
dimana,
bn
an
bn1 an1 bn x bn2 an2 bn1 x bn3 an3 bn2 x Nilai bn 3 merupakan nilai dari fungsi f(x) pada x yang bersangkutan. Sedangkan untuk menentukan nilai dari turunan dari fungsi f(x) tersebut adalah dengan meneruskan pembagian sintetik dimana nilai an , an1 , an2 , an3 diganti dengan nilai bn, bn1 , bn2 .
x
bn
bn1
bn2
(cn x) (cn1 x) cn
cn1
+
cn 2
dimana,
cn b n cn1 bn1 cn x cn2 bn2 cn1 x Nilai dari c n 2 merupakan nilai turunan dari persamaan f(x) pada x yang bersangkutan. Contoh : Misalkan kita memiliki persamaan polinom orde 3 yaitu :
f ( x) x 3 11x 2 39 x 45 Tentukan nilai fungsi dan turunannya untuk x=2 Jawab : -Secara Analitik : Untuk x = 2 maka f(2) = -3 dan f1(2) = 7 -Dengan Pembagian Sintetik maka hasilnya adalah : 2
2
1
-11
39
-45
2
-18
42
1
-9
21
-3
1
-9
21
2
–14
-7
7
1
© Avid-06
14
Fisika Komputasi
Jadi didapatkan nilai f(x) pada x=2 adalah –3, sedangkan nilai turunannya adalah 7. Hasil ini sama dengan hasil yang didapatkan secara analitik. Soal : 1. Tentukan akar dari persamaan : f(x) = x3 – 6x2 + 11x – 6 2. Tentukan akar dari persamaan : f(x) = x2 - 2x – 8 catatan : gunakan metode Iterasi Newton – Rapshon orde 1 dan orde-2
2.3.2 Akar Ganda Bila suatu fungsi mempunyai beberapa akar maka untuk mencari akar – akar yang lain dipakai iterasi Newton-Rapshon dengan titik awal yang lain juga. Jika suatu fungsi mempunyai akar ganda (lebih dari 1 akar pada titik yang sama) maka kalau hanya memakai metode Newton-Rapshon tidak akan dapat menentukan akar ganda tersebut. Contoh :
f ( x) x 3 11x 2 39 x 45 ( x 3)( x 3)( x 5)
Persamaan tersebut mempunyai akar ganda di x=3 dan akar tunggal di x=5. Untuk menyelesaikan masalah akar ganda maka kita memakai langkah-langkah berikut. Untuk suatu polinom orde n, yaitu fn(x) maka langkah-langkah mencari akar adalah : 1.Tentukan akar fn(x) paling kanan dengan metode Newton-Rapshon. Misalkan akar tersebut adalah x = x1 2.Bagi fn(x) dengan (x-x1), hasilnya polinom yang ordenya lebih rendah 1 tingkat yaitu : fn-1(x) 3.Tentukan akar paling kanan fn-1(x) dengan Newton-Rapshon. Misalkan akar tersebut adalah x=x2 4.Bagi fn-1(x) dengan (x-x2). Hasilnya adalah polinom yang ordenya lebih rendah 2 tingkat dari semula yaitu fn-2(x) 5.Langkah diteruskan sampai didapatkan polinom orde-1 contoh :
f ( x) x 3 11x 2 39 x 45 ( x 3)( x 3)( x 5) polinom orde 3
-bagilah f(x) dengan (x-5) hasilnya :
f ( x) ( x 3)( x 3) polinom orde - 2 -bagi f(x) dengan (x-3) hasilnya : © Avid-06
15
Fisika Komputasi
f ( x) ( x 3) polinom orde - 1
Adapun langkah – langkah membagi suatu polinom orde-n dengan polinom orde -1 adalah : 1.Tuliskan koefisien-koefisien dari polinom semula an
an-1 an-2 … a0
2.Lakukan pembagian sintetik untuk suatu nilai x (x didapat dari Newton-Rapshon) sehingga didapatkan : bn
bn-1 bn-2 … b0
b0 = 0 = nilai fungsi 3.Ubah bn menjadi an-1, bn-1 menjadi an-2, … b1 menjadi a0 Setelah langkah ke-3 maka akan didapatkan koefisien-koefisien polinom baru yang ordenya 1 tingkat lebih rendah yaitu : an-1 an-2 … a0 Jika langkah tersebut dilanjutkan sampai didapatkan koefisien-koefisien polinom orde1, maka akar terakhir akan didapatkan dari polinom orde 1 ini. Contoh :
f ( x) x 3 11x 2 39 x 45 orde - 3 Langkah-langkah : 1. a3 1
a2 a1
a0
2. 5
-11 39 -45 1 -11
39
-45
5
-30
45
1 -6
9
0
x=5
3. f ( x) x 2 6 x 9 orde - 2 sehingga a2 = 1 a1 = -6 a0 = 9 Ulangi langkah-langkah tsb : 1
1 -6
2. 3
9
1 -6
9
3
-9
1 -3
0
x=3
3. f ( x) x 3 orde - 1 sehingga a1 = 1 © Avid-06
16
Fisika Komputasi
a0 = -3 Jadi akar-akar dari fungsi tsb adalah : x1 = 5 x2 = 3 x3 = 3
2.4 Program Komputer 1) Program komputer iterasi Newton-Rapshon orde 1, dimana persamaan yang dicari akarnya adalah :
f ( x) x 3 11x 2 39 x 45 //********************************************************************************** // Menghitung Akar dari Fungsi dengan Metode Newton Rapshon Orde –1 // Iterasi dibatasi hanya sampai 15 kali iterasi // Compiler : Visual C++ //********************************************************************************** #include #include <math.h>
//header untuk fungsi fabs()
#include //header untuk widht,setiosflags,... void main() { double y[100]; double x,fx,dx; int n; cout <<endl; cout << "Titik awal xo : "; cin >> y[0]; cout <<endl<<endl; //cout <<setprecision(2); for (n=1; n<=15; n++) //hanya 15 iterasi {
x=y[n-1]; //***************************** fx=x*(x*x)-11*(x*x)+39*x-45; //persamaan yg akan dicari akarnya dx=3*(x*x)-22*x+39;
// Nilai turunan pertama dari persamaan
//***************************** y[n]=x-fx/dx; //rumus Newton Rapshon orde-1 cout.width(15);
//tampilkan 15 digit
cout <<setiosflags(ios::left); //rata kiri © Avid-06
17
Fisika Komputasi
cout << n ; cout.width(15); cout <<setiosflags(ios::left);
cout <
Dalam program diatas iterasi dibatasi hanya sampai 15 iterasi. Sehingga program akan berhenti begitu 15 iterasi tercapai. Program berikut ini merupakan modifikasi dari program diatas dimana iterasi akan berhenti jika selisih antara 2 nilai iterasi yang berturutan lebih kecil atau sama dengan suatu konstanta ( epsilon). Adapun implementasi dari program tersebut adalah seperti dibawah ini : //*********************************************************************** // Menghitung akar dari fungsi dengan Newton Rapshon orde-1 // iterasi dibatasi oleh selisih antar 2 nilai iterasi yg berurutan // Compiler : Visual C++ //*********************************************************************** #include #include <math.h> void main() { double eps=1e-5; //konstanta (epsilon) double y[100]; double x,fx,dx,selisih; int n; cout <<endl; cout << "Titik awal xo : "; cin >> y[0]; cout <<endl<<endl; n=0; selisih = 10; while (selisih >= eps) { n=n+1; © Avid-06
18
Fisika Komputasi
x=y[n-1]; //bentuk berikut dpt berubah sesuai bentuk fungsinya //***************************** fx=x*(x*x)-11*(x*x)+39*x-45; //Fungsi yg dicari akarnya dx=3*(x*x)-22*x+39;
//Turunan pertama dari fungsi
//***************************** y[n]=x-fx/dx;
//pers. Metode Newton-Rapshon orde-1
selisih=fabs(y[n]-x);
//selisih 2 nilai iterasi yg berurutan
}
cout << "Jumlah iterasi = "<
= " << y[n] <<endl;
}
2) Program komputer iterasi Newton-Rapshon orde 2, dimana persamaan yang dicari akarnya adalah :
f ( x) x 3 11x 2 39 x 45
// ********************************************************************** // Menghitung akar dari fungsi dengan Newton Rapshon orde-2 // iterasi dibatasi oleh selisih antar 2 nilai yg berurutan // Compiler : Visual C++ // ********************************************************************** #include #include <math.h> #include void main() { double eps=1e-5; double y[100]; double x,fx,d1,d2,pen,selisih; int n; cout <<endl; cout << "Titik awal xo : "; cin >> y[0]; cout <<endl<<endl; n=0; selisih = 10; while (selisih >= eps) //iterasi sampai selisih 2 nilai berurutan <= eps © Avid-06
19
Fisika Komputasi
{ n=n+1; x=y[n-1]; //bentuk berikut dpt berubah sesuai bentuk fungsinya //***************************** fx=x*(x*x)-11*(x*x)+39*x-45; //Fungsi yg dicari akarnya d1=3*(x*x)-22*x+39;
//Turunan pertama dari fungsi
d2=6*x-22;
//Turunan kedua dari fungsi
//***************************** pen=d1-d2*fx/(2*d1); y[n]=x-fx/pen;
//pers. Iterasi Newton-Rapshon orde-2
cout.width(10); cout << setiosflags(ios::left); cout << n;
//tampilkan iterasi ke-n
cout.width(15); cout << setiosflags(ios::left); cout << y[n];
//tampilkan nilai iterasi ke-n
selisih=fabs(y[n]-x);
//selisih 2 nilai iterasi yg berurutan
cout.width(15); cout << setiosflags(ios::left); cout << selisih; cout << endl<<endl; } cout << "Jumlah iterasi = "<
= " << y[n] <<endl;
}
Latihan : 1) Buatlah program komputer untuk menghitung persoalan gas ideal. Dimana 1 mole gas ideal yang suhunya 3000 K volumenya diperbesar pada suhu tetap. Berapa liter volume gas tersebut ketika tekanannya 15 atmosfir (tetapan gas universal 0,0823 l.atm/mole0 K). Petunjuk : pV = nRT Dimana, n=1 mole, R=0.0823 l.atm/mole0K, dan T=3000K. Sehingga : p = 24.69/V untuk menentukan berapa volume gas pada saat tekanannya 15 atmosfir maka persamaan tersebut akan menjadi : p = 24.69/V -15 Atau jika p disubstitusi menjadi y dan V menjadi x maka : y = 24.69/x –15 © Avid-06
20
Fisika Komputasi
Dengan menggunakan iterasi Newton-Rapshon maka akan dapat ditentukan solusinya. Buatlah program komputer dari masalah ini.
2) Buatlah program komputer untuk masalah difraksi pada celah. Dimana suatu celah dengan lebar 0.4 mm disinari dengan cahaya dengan panjang gelombang 5900 A. Pola difraksi yang terjadi ditangkap oleh sebuah layar yang jaraknya 70 cm dari celah tsb. Hitung jarak dari tengah-tengah terang pusat sampai intensitas cahaya tinggal ½ nya. Petunjuk : Pola difraksi pada celah dinyatakan dengan :
I
I 0 Sin 2
2
dimana, I0 intensitas cahaya yang datang dan
1 kb sin . Jika kita mensubstitusikan 2
I/I0 menjadi y dan menjadi x maka persamaan tersebut akan menjadi
y
sin 2 x x2
untuk mendapatkan nilai dimana I/I0 = 0.5 maka persamaan tersebut diubah menjadi :
y
sin 2 x 0.5 x2
Dengan memakai iterasi Newton-Rapshon maka solusi persamaan ini akan dapat dicari. Buatlah programnya!.
Tugas selanjutnya adalah bagaimana kita membuat program komputer dari Pembagian Sintetik ini. Adapun implementasi dari program komputer untuk metode Pembagian Sintetik ini dinyatakan seperti berikut ini : //************************************************************* //Menentukan nilai fungsi dan turunan dengan Metode //pembagian Sintetik //Compiler : Visual C++ //************************************************************* #include #include void pemsin(double a[100],int n,double x,double b[100]); void main() { double a[100];
© Avid-06
21
Fisika Komputasi
double b[100]; double c[100]; double x; int i,orde; cout<<"Orde dari polinom = "; cin >> orde; for (i=orde; i>=0; i--) { cout<< "a"<>a[i]; //input koefisien dari masing-masing suku } cout <<endl; cout << "tentukan titik xo : "; cin >>x; cout <<endl<<endl; pemsin(a,orde,x,b); //panggil fungsi pemsin --->pembagian sintetik for (i=orde;i>=0;i--) { if (i!=0 ) { cout.width(15); //pesan tempat 15 digit cout<<setiosflags(ios::left); //tampilkan rata kiri cout<
22
Fisika Komputasi
for(i=orde;i>=0;i--) { if (i!=1) { cout<<setw(15); cout<<setiosflags(ios::left); cout<=0; j--) { b[j]=a[j]+x*b[j+1]; } }
Jika metode Pembagian Sintetik ini kita padukan dengan iterasi Newton-Rapshon dalam penentuan akar-akar persamaan Polinom maka bentuk program komputer dari Iterasi Newton-Rapshon akan mengalami perubahan seperti di bawah ini : //****************************************************************************** //Menentukan Akar dari Polinom dengan Newton-Rapshon dimana //nilai fungsi dan turunan ditentukan dgn metode Pembagian Sintetik //Compilar : Visual C++ © Avid-06
23
Fisika Komputasi
//****************************************************************************** #include #include #include <math.h> void pemsin(double a[100],int n,double x,double b[100]); //deklarasi fungsi pembagian sintetik void main() { double eps = 1e-5; double a[100]; double b[100]; double c[100]; double y[100]; double x,fx,dx; int i,n,orde; cout<<"Orde dari polinom = "; cin >> orde; for (i=orde; i>=0; i--) { cout<< "a"<>a[i]; //input koefisien dari masing-masing suku } cout <<endl; cout << "tentukan titik awal xo : "; cin >>y[0]; cout <<endl<<endl; n=0; while (fabs(y[n]-x)>=eps) { n=n+1; x=y[n-1]; pemsin(a,orde,x,b); fx=b[0]; pemsin(b,orde,x,c); dx=c[1]; y[n]=x-fx/dx; `
} cout << "Jumlah iterasi : "<
: "<
cout <<endl; } © Avid-06
24
Fisika Komputasi
void pemsin(double a[100],int n,double x,double b[100]) { int j; b[n]=a[n]; for (j=n-1; j>=0; j--) { b[j]=a[j]+x*b[j+1]; } }
Latihan : 1) Buatlah program komputer untuk menyelesaikan masalah gerak peluru. Dimana sebuah peluru ditembakkan dengan kecepatan awal V0 = 100 m/s dan sudut elevasi = 450. Percepatan grafitasi g = 9,8 m/dt 2 .Tentukan berapa jauh peluru akan melayang di udara. Petunjuk : Persamaan lintasan peluru dinyatakan dengan : y
g x 2 (tan ) x 2 2(v0 cos )
atau bisa ditulis menjadi polinom orde-2 yaitu :
y a2 x 2 a1 x dimana, a2
g 1.10 3 m 1 2(v0 cos ) 2
a1 tan 1 a0 0 Dengan memakai iterasi Newton-Rapshon yang dikombinasikan dengan metode pembagian Sintetik maka solusi permasalahan diatas akan dapat dicari. Buat program komputernya.
© Avid-06
25
Fisika Komputasi
III INTEGRAL DAN DIFERENSIAL 3.1 INTEGRAL Nilai Integral I dari suatu fungsi f(x) menyatakan luas bidang dibawah fungsi f(x) antara x=a dan x=b atau ditulis : b
I f ( x)dx …………………………………………………………...(3.1) a
Ada beberapa metode yang dapat dipakai untuk menghitung integral diantaranya : metode persegi panjang, metode trapesium, metode Simpson, metode Romberg dll. Berikut ini akan dijelaskan beberapa metode tersebut.
3.1.1 Metode Persegi Panjang Dalam metode ini luas bidang dibawah kurva f(x) antara x = a dan x = b dapat dicari dengan membagi bidang tersebut menjadi n buah pita yang berbentuk persegi panjang, yang panjangnya f(xi) dan lebarnya x , sehingga y
Pita ke-i Pita ke-1
f(x) Pita ke-n
x
a
b
x
Gambar 3.1 Metode Persegi Panjang luas masing-masing pita dapat dinyatakan :
Li f ( xi )x …………………..……………………………………………..(3.2) Karena antara x = a dan x = b terdapat n buah pita maka luas seluruhnya menjadi : n
L f ( xi )x ………………………………………………………………...(3.3) i 1
untuk Δx sangat sempit (Δx ≈ 0) maka persamaan (3.3) adalah persamaan integral seperti yang dirumuskan dalam persamaan (3.1). Metode untuk menghitung integral seperti yang disebutkan di atas dikenal dengan metode Persegi Panjang. © Avid-06
26
Fisika Komputasi
3.1.2 Metode Trapesium Dalam metode trapesium ini luas bidang di bawah fungsi f(x) antara x = a dan x = b dapat dicari dengan membagi bidang antara x = a dan x = b menjadi n buah pita yang berbentuk trapesium yang masing-masing lebarnya x seperti diperlihatkan dalam gambar 3.2 . y
f(xi)
f(xi+1)
Pita ke-i
f(x)
Pita ke-n Pita ke-1
x
a
xi xi+1
b
x
Gambar 3.2 Metode Trapesium Dari gambar 3.2 maka luas pita ke-i yang terletak antara xi dan xi+1 adalah :
Ai
1 x[ f ( xi ) f ( xi 1 )] 2
Sehingga untuk n buah pita maka luas seluruhnya menjadi :
A
1 1 1 x[ f ( x1 ) f ( x2 )] x[ f ( x2 ) f ( x3 )] ... x[ f ( xn ) f ( xn1 )] 2 2 2
karena f(x1) = f(a) dan f(xn+1) = f(b) maka persamaan tersebut menjadi :
A
1 x[ f (a) 2 f ( x2 ) 2 f ( x3 ) 2 f ( x4 ) ... 2 f ( xn ) f (b)] ……………(3.4) 2
contoh : Jika diketahui panas jenis suatu zat seperti tabel berikut ini :
© Avid-06
t ( 0C )
c (kkal/kg 0C)
-100
0,11904
-50
0,12486
0
0,13200
50
0,14046
100
0,15024
150
0,16134
200
0,17376
27
Fisika Komputasi
Tentukan panas yg diperlukan untuk memanaskan 1 kg zat dari –1000C - 200 0C. Jawab: a) Metode Persegi Panjang n
∆x
f(xn)
1
50
0,11904
2
50
0,12486
3
50
0,13200
4
50
0,14046
5
50
0,15024
6
50
0,16134
Luas(Ln)
Jumlah
a) Metode Trapesium n
∆x
f(xn)
f(xn+1)
Luas(Ln)
1
50
0,11904
0,12486
6,0975
2
50
0,12486
0,13200
6,4215
3
50
0,13200
0,14046
6,8115
4
50
0,14046
0,15024
7,2675
5
50
0,15024
0,16134
7,7895
6
50
0,16134
0,17376
8,3775
Jumlah
42,7650
Jadi panas yang diperlukan adalah 42,7650 kkal.
3.1.3 Program Komputer metode Trapesium 1) Membuat program komputer untuk menghitung integral dengan metode Trapesium. Dimana fungsi yang akan dicari integralnya adalah : f(x) = 6 – 6x5 Adapun program komputernya adalah seperti berikut ini : //******************************************************* //Menghitung integral dengan metode Trapesium //Compiler : Visual C++ //******************************************************* #include © Avid-06
28
Fisika Komputasi
#include <math.h> #include void main() { double eps=1e-5; double trap[100]; double x,x1,x2,delt,delx,pita,fx; int i; cout << " Batas bawah : "; cin >>x1; cout << " Batas Atas : "; cin >>x2; cout <<endl<<endl; i=0; delt=100;
//delt = selisih antara 2 hasil integral(iterasi) yg berurutan
while (fabs(delt)>eps)
//iterasi selama delt lebih besar dari eps
{ i=i+1; trap[i]=0; x=x1; pita=pow(2,i-1); //2^i-1 (pow = pangkat) delx=(x2-x1)/pita; //pita= jumlah pita while (x<x2)
//delx=lebar masing-masing pita;
{ //*************************** fx=6-6*pow(x,5);
//bentuk pers. 6-6x^5
//*************************** //fx=fabs(x); //contoh : pers. garis y=x //fx=1/x ; //contoh : pers. flux magnet pada kawat berarus if (x==x1 || x==x2) trap[i]=trap[i]+fx; else trap[i]=trap[i]+2*fx; x=x+delx; } trap[i]=trap[i]*delx/2; if (i==1) { cout<<setw(10); //10 lokasi digit cout<<setiosflags(ios::left); //rata kiri cout<< pita; © Avid-06
29
Fisika Komputasi
cout<<setw(15); cout<<setiosflags(ios::left); cout<< trap[i]; } else { delt=trap[i]-trap[i-1]; cout<<setw(10); cout<<setiosflags(ios::left); cout<
Coba jalankan program tersebut dengan batas-batas integral x=0 dan x=1. Bagaimanakah hasilnya.Jika programnya benar anda akan mendapatkan hasil integrasi = 4.99999761 Latihan : 1) Buatlah program komputer untuk menghitung integral dari fungsi : f(x) = 4 Eksekusi program tersebut dengan memasukkan batas-batas integrasi x = 0 dan x = 5. Jika programnya benar maka akan didapatkan hasil integrasi = 20. Program ini mengandung arti yaitu kita mencari luas segi empat dengan panjang 5 dan lebar 4. f(x) 4 f(x) = 4
x 0 5 2) Coba buat program komputer untuk menghitung integral dari fungsi : © Avid-06
30
Fisika Komputasi
f(x) = 1/x Coba eksekusi program tersebut dengan batas-batas integrasi adalah x=0.01 dan x=0.09. Jika programnya benar maka anda akan mendapatkan harga integrasi = 2.19723
3.1.4 Metode Simpson Kesalahan – kesalahan yang terjadi pada metode persegi panjang maupun metode trapesium adalah karena kurva antara x = a dan x = b didekati dengan potongan – potongan garis lurus. Kesalahan ini dapat diperkecil jika kita tidak menggunakan garis lurus tetapi kurva lain. Salah satu pendekatan yang dipakai adalah parabola atau polinom orde dua. Untuk menurunkan rumus integrasinya maka dimisalkan xi+1 ditempatkan di x = 0 sehingga xi berada di x = -∆x serta xi+2 berada di x = ∆x. Jika persamaan parabola yang dipakai adalah : f ( x) a2 x 2 a1 x a0 ………………………………………………………..(3.5)
maka luas pasangan pita di bawah parabola : x
Ai
(a
2
x 2 a1 x a0 )dx …………………………………………………....(3.6)
x
∆x
1 1 a 2 x 3 a1 x 2 a0 x 3 2 -∆x
2 a 2 x 3 2a0 x …………………………………………………………..(3.7) 3
Nilai fungsi pada ketiga titik tersebut adalah :
f ( xi ) a2 (x) 2 a1 (x) a0 ……………………………………………...(3.8)
f ( xi 1 ) a0 ……………………………………………………………………(3.9) f ( xi 2 ) a2 (x) 2 a1 (x) a0 ……………………………………………...(3.10)
Dari persamaan (3.8), (3.9), dan (3.10) maka didapatkan koefisien-koefisien : a2
f ( xi ) 2 f ( xi 1 ) f ( xi 2 ) ……………………………………………..(3.11) 2x 2
a1
f ( xi 2 ) f ( xi ) ……………………………………………………….…(3.12) 2x
a0 f ( xi 1 ) …………………………………………………………………..(3.13) © Avid-06
31
Fisika Komputasi
Sehingga persamaan (3.7) menjadi :
1 Ai x[ f ( xi 4 f ( xi 1 ) f ( xi 2 )] …………………………………………(3.14) 3 Maka luas semua pasangan pita adalah :
1 A x[ f ( x1 4 f ( x2 ) f ( x3 )] 3 1 x[ f ( x3 4 f ( x4 ) f ( x5 )] ... 3 1 x[ f ( xn 2 4 f ( xn1 ) f ( xn )] ……………………………………....(3.15) 3 atau
1 A x[ f ( x1 ) 4 f ( x2 ) 4 f ( x6 ) 4 f ( x8 ) ... 3 2 f ( x3 ) 2 f ( x5 ) 2 f ( x7 ) ... f ( xn )] ………………………………(3.16) contoh : 1)Jika diketahui panas jenis suatu zat seperti tabel berikut ini : t ( 0C )
c (kkal/kg 0C)
-100
0,11904
-50
0,12486
0
0,13200
50
0,14046
100
0,15024
150
0,16134
200
0,17376
Tentukan panas yang diperlukan untuk memanaskan 1 kg zat dari –1000C sampai 200 0C. Jawab: n
∆x
f(xn)
f(xn+1)
f(xn+2)
Luas(Ln)
1
50
0,11904
0,12486
0,13200
12,5080
2
50
0,13200
0,14046
0,15024
14,0680
3
50
0,15024
0,16134
0,17376
16,1560
Jumlah
42,7320
Jadi panas yang diperlukan adalah 42,7320 kkal 2)Suatu kawat yang sangat panjang dialiri arus listrik 30 A. Hitung besar flux magnet yang menembus suatu persegi panjang yang sisi panjangnya sejajar dan © Avid-06
32
Fisika Komputasi
sisi pendeknya tegak lurus kawat. Panjang persegi panjang 30 cm, lebar 8 cm. Jarak salah satu sisi panjang persegi panjang ke kawat adalah 1 cm ( μ0 = 4π . 10-7 weber/amp.m) 1 cm
r
8 cm
L 30 cm
dr
pita = dA = Ldr
Jawab : Besar flux magnet : d B.dA
Ambil suatu pita sejajar dengan kawat berarus panjang L = 30 cm, dan lebarnya dr. Maka luas pita dA = L.dr. Sehingga flux magnet yang menembus pita ini adalah :
B
0 i 2 r
Maka
d
0 i iL dr L dr 0 2 r 2 r
0 iL 0,091 dr 2 0,01r 3.1.5 Program komputer metode Simpson 1) Membuat program komputer untuk menghitung integral dengan metode Simpson. Dimana fungsi yang akan dicari integralnya adalah : f(x) = 6 – 6x5 Adapun program komputernya adalah seperti berikut ini : //******************************************************* //Menghitung integral dengan metode Simpson //Compiler : Visual C++ //******************************************************* #include #include <math.h> #include void main() { double eps=1e-5; © Avid-06
33
Fisika Komputasi
double simp[100]; double x,x1,x2,delt,delx,pita,fx; int i,j; cout << " Batas bawah : "; cin >>x1; cout << " Batas Atas : "; cin >>x2; cout <<endl<<endl; i=0; delt=100; //delt = selisih antara 2 hasil integral(iterasi) yg berurutan while (fabs(delt)>eps) //iterasi selama delt lebih besar dari eps { i=i+1; j=0; simp[i]=0; x=x1; pita=2*pow(2,i-1); //2^i-1 (pow = pangkat) delx=(x2-x1)/pita; //pita= jumlah pita while (x<x2)
//delx=lebar masing-masing pita;
{ //*************************** fx=6-6*pow(x,5);
//contoh : bentuk pers. 6-6x^5
//*************************** //fx=fabs(x); //contoh : pers. garis y=x //fx=1/x ; //contoh : pers. flux magnet pada kawat berarus j=j+1; if (x==x1 || x==x2) simp[i]=simp[i]+fx; else { if (x>x1 && x<x2 && j % 2 ==0) //suku genap simp[i]=simp[i]+4*fx; if (x>x1 && x<x2 && j % 2 ==1) //suku ganjil simp[i]=simp[i]+2*fx; } x=x+delx; } simp[i]=simp[i]*delx/3; if (i==1) { cout<<setw(10); //10 lokasi digit © Avid-06
34
Fisika Komputasi
cout<<setiosflags(ios::left); //rata kiri cout<< pita; cout<<setw(15); cout<<setiosflags(ios::left); cout<< simp[i]; } else { delt=simp[i]-simp[i-1]; cout<<setw(10); cout<<setiosflags(ios::left); cout<
Latihan : 1) Coba buat program yang sama seperti pada latihan untuk metode Trapesium. Bandingkan hasilnya.
3.2 DIFERENSIAL Diferensial atau sering disebut turunan dapat dihitung dengan memakai uraian deret Taylor. Jika uraian deret Taylor di sekitar x dinyatakan dengan f(x+h) dan f(x-h). Dimana masing-masing dinyatakan dengan persamaan :
1 1 f ( x h) f ( x) hf ' ( x) h 2 f '' ( x) h 3 f ''' ( x) ... 2 6
…………………..(3.17)
1 1 f ( x h) f ( x) hf ' ( x) h 2 f '' ( x) h 3 f ''' ( x) ... .....…………………..(3.18) 2 6 3.2.1 Turunan Pertama © Avid-06
35
Fisika Komputasi
Jika kita mengambil selisih antara kedua persamaan tersebut (3.17 dan 3.18) maka akan didapatkan :
f ' ( x)
f ( x h) f ( x h) 1 2 ''' h f ( x) ... .......…………………………...(3.19) 2h 6
Dari persamaan (3.19) jika kita mengambil nilai h yang sangat kecil maka suku-suku dengan h pangkat 2 atau lebih bisa diabaikan. Sehingga akan didapatkan persamaan Turunan Pertama dari suatu fungsi f(x) yaitu :
f ' ( x)
f ( x h) f ( x h) ………………………………………………..(3.20) 2h
3.2.2 Turunan Kedua Untuk mendapatkan Turunan Kedua dari fungsi f(x) maka kita menjumlahkan kedua persamaan f(x+h) dan f(x-h) dalam persamaan (3.17) dan 3.18). Sehingga didapatkan persamaan yaitu :
f '' ( x)
f ( x h) 2 f ( x ) f ( x h) 1 2 h f 12 h2
IV
( x) ...
...……………….(3.21)
Bila h sangat kecil maka suku dengan h pangkat 2 atau lebih bisa diabaikan. Sehingga akan didapatkan persamaan Turunan Kedua dari f(x) yaitu :
f '' ( x)
f ( x h) 2 f ( x ) f ( x h) ………………………………………(3.22) h2
Kedua metode untuk mencari turunan diatas disebut dengan metode Beda Sentral (central difference). Contoh : 1.Tentukan turunan pertama dari fungsi y = x2 saat x = 1 jawab :
f ' ( x)
f ( x h) f ( x h) 2h
a)Secara analitik
dy dx
2x 2 x 1
b)Secara Numerik (Metode Beda Sentral) © Avid-06
36
Fisika Komputasi
h=1
;
hn+1 = hn / 2
n
hn
f(x+h)
f(x-h)
f1(x) =f(x+h)-f(x-h)/2h
1
1,00
4,0000
0,0000
2,00
2
0,50
2,2500
0,2500
2,00
3
0,25
1,5625
0,5625
2,00
2.Tentukan turunan pertama dan kedua dari fungsi y = x2 + x saat x = 1 Jawab : Turunan Pertama : a)Secara analitik
dy dx
2x 1 3 x 1
b)Secara Numerik (Metode Beda Sentral) h=1
;
hn+1 = hn / 2
n
hn
f(x+h)
f(x-h)
f1(x) =f(x+h)-f(x-h)/2h
1
1,00
6,0000
0,0000
3,00
2
0,50
3,7500
0,7500
3,00
Turunan Kedua : a)Secara analitik
dy dx
2 x 1
b)Secara Numerik (Metode Beda Sentral) h=1
;
hn+1 = hn / 2
n
hn
f(x)
f(x+h)
f(x-h)
f11(x)
1
1,00
2
6,0000
0,0000
2,00
2
0,50
2
3,7500
0,7500
2,00
3.2.3 Program Komputer 1)Program komputer untuk menghitung Turunan Pertama dari fungsi. Dimana fungsi yang dicari turunannya adalah : © Avid-06
37
Fisika Komputasi
f(x) = x2– 5x Adapun implementasi program komputernya adalah seperti berikut ini : //*********************************************************************************** //Menghitung Turunan pertama suatu fungsi dengan metode beda sentral //compiler : Visual C++ //*********************************************************************************** #include #include <math.h> #include void main() { double eps=1e-5; double z[10]; double fx,del,dx,zz,x,y,h; int i,n; cout << " Masukkan nilai x : "; cin >>y; cout <<endl<<endl; h=1; n=0; zz=0; del=10; while (fabs(del)>=eps) //iterasi selama del lebih besar dari eps { n+=1; for (i=1;i<=2;i++) { if (i==1) x=y+h; else x=y-h; //fx=9.8*68.1*(1-exp(-12.5*x/68.1)); //fx=fx/12.5; //***************** fx=x*x-5*x; //fungsi yang dicari turunannya : x^2 - 5x //***************** z[i]=fx; } dx=(z[1]-z[2])/(2*h); del=zz-dx; cout.width(15); cout<
38
Fisika Komputasi
Coba jalankan (eksekusi) program tersebut. Jika anda masukkan nilai x =1 maka akan didapatkan hasil turunan pertama adalah –3. Cobalah untuk memasukkan nilai x yang lain dan lihat hasilnya. Bandingkan dengan perhitungan secara manual. 2) Membuat program komputer untuk menghitung Turunan kedua dari fungsi. Dimana fungsi yang dicari turunannya adalah : f(x) = x3– 5x Adapun implementasi program komputernya adalah seperti berikut ini : //****************************************************************** //Menghitung Turunan pertama suatu fungsi dengan metode beda sentral //compiler : Visual C++ //****************************************************************** #include #include <math.h> #include void main() { double eps=1e-5; double z[10]; double fx,del,dx,zz,x,y,h; int i,n; cout << " Masukkan nilai x : "; cin >>y; cout <<endl<<endl; h=1; n=0; zz=0; del=10; while (fabs(del)>=eps) //iterasi selama del lebih besar dari eps { n+=1; for (i=1;i<=3;i++) { switch(i) { case (1): { x=y+h; break; } case(2): { x=y; break; } case(3): { x=y-h; break; } } //fx=9.8*68.1*(1-exp(-12.5*x/68.1)); //fx=fx/12.5; //***************** fx=x*x*x-5*x; //fungsi yang dicari turunannya : x^3 - 5x //***************** © Avid-06
39
Fisika Komputasi
z[i]=fx; } dx=(z[1]-2*z[2]+z[3])/(h*h); del=zz-dx; cout.width(15); cout< yo cin >>x1>>y1; cout << " Masukkan titik akhir (x) : "; cin >>x2; cout <<endl<<endl; n=0; delt=100; while (fabs(delt)>eps) //iterasi selama delt lebih besar dari eps { n+=1; pita=pow(2,n); //2^n delx=(x2-x1)/pita; x=x1; y=y1; while (x<x2) { //************************* //fx=y; fx=4*x+5; //************************* y=y+fx*delx; //iterasi Euler orde-1 x=x+delx; } if (n==1) { cout.width(15); cout<
Coba jalankan program tersebut dan berikan masukan syarat awal (x0,y0) = (1 1), serta berikan masukan untuk nilai x = 2. Jika programnya berjalan dengan benar maka anda akan mendapatkan output berupa nilai fungsi di x = 2 adalah 11.999. 4.2.3 Metode Euler orde 2 © Avid-06
45
Fisika Komputasi
Dari uraian deret Taylor untuk f(x+h), jika kita menyertakan suku-suku dengan h pangkat 2 maka akan didapatkan persamaan :
f ( x h) f ( x) f ' ( x)h 1 f ' ' ( x)h 2 ...................................................(4.3) 2 dimana f ' ( x) f ( x, y)
f ' ' ( x) f ' ( x, y) f x ( x, y) f y ( x, y) f ( x, y) persamaan (4.3) dapat dinyatakan dalam bentuk iterasi :
yi 1 yi f ( x, y)h 1 f ' ( x, y)h 2 ...................................................................(4.4) 2 Persamaan ini disebut sebagai iterasi Euler orde 2. 4.2.4 Program komputer Metode Euler orde 2 Program komputer untuk mencari solusi persamaan diferensial biasa dengan iterasi Euler orde-2, dimana bentuk persamaan diferensial yang dicari solusinya adalah : f ’(x) = 4x+5 Adapun implementasi program komputernya adalah seperti berikut ini : //*************************************************************************** //Menghitung Persamaan Difrensial dengan metode Euler orde dua //compiler : Visual C++ //*************************************************************************** #include #include <math.h> #include void main() { double eps=1e-3; double x,x1,x2,y,y1,y2,fx,dx,delx,delt,pita; int n; cout << " Masukkan titik awal (xo,yo) : "; //masukkan nilai xo <spasi> yo cin >>x1>>y1; cout << " Masukkan titik akhir (x) : "; cin >>x2; cout <<endl<<endl; n=0; delt=100; while (fabs(delt)>eps) //iterasi selama delt lebih besar dari eps { n+=1; pita=pow(2,n); //2^n delx=(x2-x1)/pita; x=x1; y=y1; while (x<x2) { //************************ //fx=y; //bentuk fungsi //dx=y; //turunan fx=4*x+5; © Avid-06
46
Fisika Komputasi
dx=4; //************************* y=y+fx*delx+dx*pow(delx,2)/2; //iterasi Euler orde 2 x=x+delx; } if (n==1) { cout.width(15); cout<
adalah : "<
}
Coba jalankan program tersebut dan berikan masukan syarat awal (x0,y0) = (1 1), serta berikan masukan untuk nilai x = 2. Jika programnya berjalan dengan benar maka anda akan mendapatkan output berupa nilai fungsi di x = 2 adalah 11.9993.
4.3 Metode Runge -Kutta Persamaan umum iterasi Runge-Kutta adalah :
yn1 y n f ( xn , y n , h).......................................................................................(4.5) dimana
f ( xn , y n , h) a1k1 a2 k 2 ... am k m h xn1 xn dengan
k1 hf ( xn , y n ) k 2 hf ( xn p1h, y n q11k1 k3 hf ( xn p2 h, y n q21k1 q22k 2 )
...
...
...
k m hf ( xn pm1h, y n qm1,1k1 qm1, 2 k 2 ... qm1,m1k m1
4.3.1 Runge-Kutta orde 1 © Avid-06
47
Fisika Komputasi
Jika diambil m = 1 maka akan menghasilkan iterasi Runge-Kutta orde 1 yaitu :
yn1 yn a1k1 ..................................................................................................(4.6) yn1 y n a1hf ( xn , y n )....................................................................................(4.7) karena f ( xn , y n , h) a1k1 dan k1 hf ( xn , y n ) . Dengan mengambil a1 1 maka persamaan (4.7) sama dengan persamaan iterasi Euler orde-1 yaitu :
y n1 y n hf ( xn , y n ) 4.3.2 Runge-Kutta orde 2 Jika diambil m = 2 maka akan menghasilkan iterasi Runge-Kutta orde 2 yaitu :
yn1 yn a1k1 a2 k 2 .......................................................................................(4.8) dimana
k1 hf ( xn , y n ) k 2 hf ( xn p1h, y n q11k1 ) jadi dalam hal ini kita perlu menentukan nilai-nilai dari a1 , a2 , p1 , dan q11 . Untuk mendapatkan nilai dari konstanta-konstanta tersebut maka pertama-tama kita lihat uraian deret Taylor sampai orde 2 yaitu :
f ( x h) f ( x) f ' ( x)h 1 f ' ' ( x, y)h 2 ...................................................(4.8) 2 dimana f ' ( x) f ( x, y)
f ' ' ( x) f ' ( x, y) f x ( x, y) f y ( x, y) f ( x, y)
sehingga persamaan (4.7) menjadi :
f ( x h) f ( x) f ( x, y)h 1
f ( x, y) f y ( x, y) f ( x, y) h 2 ....................(4.9) 2 x Jika dinyatakan dengan rumus iterasi maka persamaan (4.8) menjadi :
1 f x ( xn , y n ) f y ( xn , y n ) f ( xn , y n ) h 2 ......................(4.10) 2 untuk k2 jika diuraikan dalam deret Taylor akan menjadi : y n1 y n f ( xn , y n )h
k 2 f ( xn , y n ) f x ( xn , y n ) p1h f y ( xn , y n )q11k1 h
sehingga rumus iterasi (persamaan (4.10)) menjadi : y n1 y n f ( xn , y n )a1h f ( xn , y n )a2 h f x ( xn , y n )a2 p1h 2
f y ( xn , y n ) f ( xn , y n )a2 q11h 2
y n f ( xn , y n )a1 f ( xn , y n )a2 h f x ( xn , y n )a2 p1 f y ( xn , y n ) f ( xn , y n )
a2 q11 h 2 .................................................................................................(4.11) © Avid-06
48
Fisika Komputasi
dari persamaan (4.10) dan (4.11) didapatkan :
a1 a2 1........................................................................................................(4.12) a2 p1 1 2 ..........................................................................................................(4.13) a2 q11 1 2 .........................................................................................................(4.14) a) Metode Heun Dalam metode ini diambil nilai dari konstanta :
a1 a2
1
2
p1 q11 1 Sehingga persamaan iterasi Runge-Kutta orde 2 menjadi :
yn1 y n 1 2 k1 1 2 k 2 .....................................................................................(4.15) dimana
k1 hf ( xn , y n ) k 2 hf ( xn h, y n k1 ) b) Metode Raltson Dalam metode ini diambil nilai dari konstanta :
a1
1
a2 2 3
3
p1 q11 3 4 Sehingga persamaan iterasi Runge-Kutta orde 2 menjadi :
yn1 y n 13 (k1 2k 2 )....................................................................................(4.16) dimana
k1 hf ( xn , y n ) k 2 hf ( xn 3 4 h, y n 3 4 k1 ) c) Metode Poligon Dalam metode ini diambil nilai dari konstanta :
a1 0
p1 q11
a2 1 1
2
Sehingga persamaan iterasi Runge-Kutta orde 2 menjadi :
yn1 y n k 2 ...................................................................................................(4.17) dimana
k1 hf ( xn , y n ) k 2 hf ( xn 1 2 h, y n 1 2 k1 ) Contoh : © Avid-06
49
Fisika Komputasi
Hitunglah solusi persamaan difrensial berikut pada x = 2 dengan metode Heun:
dy 4x 5 dx dimana syarat awal (x0,y0) = (1,1) Jawab : Jika kita membagi antara x=1 dan x=2 menjadi 10 pita maka lebar masing-masing pita (h) adalah : h
x 2 x1 2 1 0.1 pita 10
y n1 y n 1 2 k1 1 2 k 2 dimana
k1 hf ( xn , y n ) k 2 hf ( xn h, y n k1 ) sehingga :
y1 y0 12 k1 12 k 2 k1 0.1.(4.x0 5) 0.1(4.1 5) 0.9 k 2 0.1.(4.( x0 h) 5) 0.1(4.(1 0.1) 5) 0.1(9.4) 0.94 jadi
y1 y0 12 k1 12 k 2 1 12 (0.9) 12 (0.94)
y1 1 0.45 0.47 1.92 hasil selengkapnya dapat dilihat dalam tabel berikut ini. i
© Avid-06
xi
yi
k1
k2
yi+1
0
1
1
0.9
0.94
1.92
1
1.1
1.92
0.94
0.98
2.88
2
1.2
2.88
0.98
1.02
3.88
3
1.3
3.88
1.02
1.06
4.92
4
1.4
4.92
1.06
1.10
6.00
5
1.5
6.00
1.10
1.14
7.12
6
1.6
7.12
1.14
1.18
8.28
7
1.7
8.28
1.18
1.22
9.48
8
1.8
9.48
1.22
1.26
10.72
9
1.9
10.72
1.26
1.30
12.00
50
Fisika Komputasi
Jadi solusi persamaan difrensial pada x = 2 adalah 12.00 4.3.3 Program komputer Metode Runge Kutta orde 2 - Heun 1) Membuat program komputer untuk mencari solusi persamaan diferensial biasa dengan iterasi Runge-Kutta orde 2 berdasarkan metode Heun. Dimana bentuk persamaan diferensial yang dicari solusinya adalah : f ’(x) = 4x+5 Adapun implementasi program komputernya adalah seperti berikut ini : //******************************************************************************************** //Menghitung Persamaan Diferensial dengan metode Runge-kutta orde dua-Heun //compiler : Visual C++ //******************************************************************************************** #include #include <math.h> #include void main() { double eps=1e-3; double k[4]; double xx,x,x1,x2,yy,y,y1,y2,fx,delx,delt,pita; int i,n; cout << " Masukkan titik awal (xo,yo) : "; //masukkan nilai xo <spasi> yo cin >>x1>>y1; cout << " Masukkan titik akhir (x) : "; cin >>x2; cout <<endl<<endl; n=0; delt=100; while (fabs(delt)>eps) //iterasi selama delt lebih besar dari eps { n+=1; pita=pow(2,n); //2^n delx=(x2-x1)/pita; xx=x1; yy=y1; while (xx<x2) { x=xx; y=yy; for (i=1;i<=2;i++) { //************************ //fx=y; //bentuk fungsi fx=4*x+5; //fx=1-exp(-12.5*x/68.1); //fx=fx*9.8*68.1/12.5; //************************* k[i]=delx*fx; if (i==1) { x=xx+delx; y=yy+k[i]; } © Avid-06
51
Fisika Komputasi
} //for xx=xx+delx; yy=yy+(k[1]+k[2])/2; } //while if (n==1) { cout.width(15); cout<
}
Coba jalankan program tersebut dan berikan masukan syarat awal (x0,y0) = (1 1), serta berikan masukan untuk nilai x = 2. Jika programnya berjalan dengan benar maka anda akan mendapatkan output berupa nilai fungsi di x = 2 adalah 12.
2) Membuat program komputer untuk mencari solusi persamaan diferensial biasa dengan iterasi Runge-Kutta orde 2 berdasarkan metode Raltson. Dimana bentuk persamaan diferensial yang dicari solusinya adalah : f ’(x) = 4x+5
Adapun implementasi program komputernya adalah seperti berikut ini : //********************************************************************************************** //Menghitung Persamaan Diferensial dengan metode Runge-kutta orde dua-Raltson //compiler : Visual C++ //********************************************************************************************** #include #include <math.h> #include void main() { double eps=1e-3; double k[4]; double xx,x,x1,x2,yy,y,y1,y2,fx,delx,delt,pita; int i,n; cout << " Masukkan titik awal (xo,yo) : "; //masukkan nilai xo <spasi> yo cin >>x1>>y1; © Avid-06
52
Fisika Komputasi
cout << " Masukkan titik akhir (x) : "; cin >>x2; cout <<endl<<endl; n=0; delt=100; while (fabs(delt)>eps) //iterasi selama delt lebih besar dari eps { n+=1; pita=pow(2,n); //2^n delx=(x2-x1)/pita; xx=x1; yy=y1; while (xx<x2) { x=xx; y=yy; for (i=1;i<=2;i++) { //************************ //fx=y; //bentuk fungsi fx=4*x+5; //fx=1-exp(-12.5*x/68.1); //bentuk fungsi //fx=fx*9.8*68.1/12.5; //************************* k[i]=delx*fx; if (i==1) { x=xx+3*delx/4; y=yy+3*k[i]/4; } } //for xx=xx+delx; yy=yy+(k[1]+2*k[2])/3; //iterasi runge-kutta orde 2 - Raltson } //while if (n==1) { cout.width(15); cout<
} © Avid-06
53
Fisika Komputasi
Coba jalankan program tersebut dan berikan masukan syarat awal (x0,y0) = (1 1), serta berikan masukan untuk nilai x = 2. Jika programnya berjalan dengan benar maka anda akan mendapatkan output berupa nilai fungsi di x = 2 adalah 12.
3) Membuat program komputer untuk mencari solusi persamaan diferensial biasa dengan iterasi Runge-Kutta orde 2 berdasarkan metode Poligon. Dimana bentuk persamaan diferensial yang dicari solusinya adalah : f ’(x) = 4x+5
Adapun implementasi program komputernya adalah seperti berikut ini : //****************************************************************** //Menghitung Persamaan Diferensial dengan metode Runge-kutta orde dua-Poligon //compiler : Visual C++ //****************************************************************** #include #include <math.h> #include void main() { double eps=1e-3; double k[4]; double xx,x,x1,x2,yy,y,y1,y2,fx,delx,delt,pita; int i,n; cout << " Masukkan titik awal (xo,yo) : "; //masukkan nilai xo <spasi> yo cin >>x1>>y1; cout << " Masukkan titik akhir (x) : "; cin >>x2; cout <<endl<<endl; n=0; delt=100; while (fabs(delt)>eps) //iterasi selama delt lebih besar dari eps { n+=1; pita=pow(2,n); //2^n delx=(x2-x1)/pita; xx=x1; yy=y1; while (xx<x2) { x=xx; y=yy; for (i=1;i<=2;i++) { //************************ //fx=y; //bentuk fungsi fx=4*x+5; //fx=1-exp(-12.5*x/68.1); //bentuk fungsi //fx=fx*9.8*68.1/12.5; //************************* k[i]=delx*fx; if (i==1) © Avid-06
54
Fisika Komputasi
{ x=xx+delx/2; y=yy+k[i]/2; } } //for xx=xx+delx; yy=yy+k[2]; //iterasi runge-kutta orde 2 - Poligon } //while if (n==1) { cout.width(15); cout<
}
Coba jalankan program tersebut dan berikan masukan syarat awal (x0,y0) = (1 1), serta berikan masukan untuk nilai x = 2. Jika programnya berjalan dengan benar maka anda akan mendapatkan output berupa nilai fungsi di x = 2 adalah 12.
4.3.4 Runge-Kutta orde 3 Rumus iterasi dari metode ini didapatkan dengan menyertakan suku-suku orde-3 dalam uraian deret Taylor. Jika diambil m=3 dalam persamaan (4.1) maka didapat rumus iterasi Runge-Kutta orde 3 adalah :
yn1 y n 16 (k1 4k 2 k3 )............................................................................(4.18) dengan
k1 hf ( xn , y n ) k 2 hf ( xn 1 2 h, y n 1 2 k1 ) k3 hf ( xn h, y n k1 2k 2 ) 4.3.5 Program Komputer Runge -Kutta orde 3
© Avid-06
55
Fisika Komputasi
1) Membuat program komputer untuk mencari solusi persamaan diferensial biasa dengan iterasi Runge-Kutta orde 3. Dimana bentuk persamaan diferensial yang dicari solusinya adalah : f ’(x) = 4x+5 Adapun implementasi program komputernya adalah seperti berikut ini : //************************************************************************************** //Menghitung Persamaan Diferensial dengan metode Runge-Kutta orde Tiga //compiler : Visual C++ //************************************************************************************** #include #include <math.h> #include void main() { double eps=1e-3; double k[4]; double xx,x,x1,x2,yy,y,y1,y2,fx,delx,delt,pita; int i,n; cout << " Masukkan titik awal (xo,yo) : "; //masukkan nilai xo <spasi> yo cin >>x1>>y1; cout << " Masukkan titik akhir (x) : "; cin >>x2; cout <<endl<<endl; n=0; delt=100; while (fabs(delt)>eps) //iterasi selama delt lebih besar dari eps { n+=1; pita=pow(2,n); //2^n delx=(x2-x1)/pita; xx=x1; yy=y1; while (xx<x2) { x=xx; y=yy; for (i=1;i<=3;i++) { //************************ //fx=y; //bentuk fungsi fx=4*x+5; //fx=1-exp(-12.5*x/68.1); //bentuk fungsi f//x=fx*9.8*68.1/12.5; //************************* k[i]=delx*fx; switch(i) { case (1): { x=xx+delx/2; y=yy+k[i]/2; break; } case(2): © Avid-06
56
Fisika Komputasi
{ x=xx+delx/2; y=yy-k[1]+2*k[2]; break; } } } //for xx=xx+delx; yy=yy+(k[1]+4*k[2]+k[3])/6; //iterasi Runge-Kutta orde 3 } //while if (n==1) { cout.width(15); cout<
}
Coba jalankan program tersebut dan berikan masukan syarat awal (x0,y0) = (1 1), serta berikan masukan untuk nilai x = 2. Jika programnya berjalan dengan benar maka anda akan mendapatkan output berupa nilai fungsi di x = 2 adalah 11.993.
4.3.6 Runge-Kutta orde 4 Metode ini dinyatakan dengan menyertakan suku-suku orde-4 dalam uraian deret Taylor. Adapun rumus iterasi Runge-Kutta orde 4 adalah :
yn1 y n 16 (k1 2k 2 2k3 k 4 )..................................................................(4.19) dengan
k1 hf ( xn , y n ) k 2 hf ( xn 1 2 h, y n 1 2 k1 ) k3 hf ( xn 1 2 h, y n 1 2 k 2 ) k 4 hf ( xn h, y n k 3 ) © Avid-06
57
Fisika Komputasi
4.3.7 Program Komputer Runge-Kutta orde 4 1) Membuat program komputer untuk mencari solusi persamaan diferensial biasa dengan iterasi Runge-Kutta orde 4. Dimana bentuk persamaan diferensial yang dicari solusinya adalah : f ’(x) = 4x+5 Adapun implementasi program komputernya adalah seperti berikut ini : //****************************************************************** //Menghitung Persamaan Diferensial dengan metode Runge-Kutta orde Empat //compiler : Visual C++ //****************************************************************** #include #include <math.h> #include void main() { double eps=1e-3; double k[4]; double xx,x,x1,x2,yy,y,y1,y2,fx,delx,delt,pita; int i,n; cout << " Masukkan titik awal (xo,yo) : "; //masukkan nilai xo <spasi> yo cin >>x1>>y1; cout << " Masukkan titik akhir (x) : "; cin >>x2; cout <<endl<<endl; n=0; delt=100; while (fabs(delt)>eps) //iterasi selama delt lebih besar dari eps { n+=1; //cacah iterasi pita=pow(2,n); //2^n delx=(x2-x1)/pita; xx=x1; yy=y1; while (xx<x2) { x=xx; y=yy; for (i=1;i<=4;i++) { //************************ //fx=y; //bentuk fungsi fx=4*x+5; //fx=1-exp(-12.5*x/68.1); //bentuk fungsi //fx=fx*9.8*68.1/12.5; //************************* k[i]=delx*fx; if(i==1||i==2) { x=xx+delx/2; y=yy+k[i]/2; } else if (i==3) © Avid-06
58
Fisika Komputasi
{ x=xx+delx; y=yy+k[3]; } } //for xx=xx+delx; yy=yy+(k[1]+2*(k[2]+k[3])+k[4])/6; //iterasi Runge-Kutta orde 4 } //while if (n==1) { cout.width(15); cout<
}
Coba jalankan program tersebut dan berikan masukan syarat awal (x0,y0) = (1 1), serta berikan masukan untuk nilai x = 2. Jika programnya berjalan dengan benar maka anda akan mendapatkan output berupa nilai fungsi di x = 2 adalah 12.
© Avid-06
59
Fisika Komputasi
V MENYELESAIKAN SISTEM PERSAMAAN LINEAR 5.1 Sistem Persamaan Linear Persamaan Linear adalah suatu fungsi atau persamaan dimana pangkat tertinggi dari variabel-variabelnya adalah 1, misalnya :
ax by cz k Dimana x,y,z
= variabel
a,b,c
= koefisien
k
= konstanta
Sedangkan yang dimaksud dengan Sistem persamaan Linear adalah beberapa persamaan linear yang saling berhubungan dimana banyaknya koefisien, variabel, dan konstanta adalah sama. Agar sistem persamaan linear tersebut dapat dicari solusinya (dapat dihitung nilai variabelnya) maka banyaknya persamaan dalam sistem persamaan linear harus sama dengan banyaknya variabel. Contoh sistem persamaan linear : x 3y z 7 4 x 4 y 3z 15
5x 8 y 2 z 25 Atau persamaan diatas bisa ditulis :
x1 3x2 x3 7 4 x1 4 x2 3x3 15 5x1 8x2 2 x3 25 Tampak dalam sistem persamaan linear diatas banyaknya persamaan dan variabel adalah sama yaitu 3. Secara umum sistem persamaan linear dinyatakan dengan :
a11x1 a12 x2 a13 x3 ... a1n xn b1 a21x1 a22 x2 a23 x3 ... a2n xn b2 a31x1 a32 x2 a33 x3 ... a3n xn b3 .................................................................(5.1)
.......................................................... am1 x1 am 2 x2 am3 x3 ... amn xn bm Sistem persamaan linear tersebut terdiri dari m persamaan dan n variabel. Persamaan linear tersebut akan dapat diselesaikan jika jumlah m = n. Penyelesaian persaman linear © Avid-06
60
Fisika Komputasi
adalah mencari nilai xi (i = 1,2,...,n). Sistem persamaan linear (5.1) dapat dinyatakan dalam bentuk matrik yaitu :
a11 a 21 ... a m1
a12 a 22 ... am2
... a1n x1 b1 ... a 2 n x 2 b2 ..........................................................(5.2) ... ... ... ... ... a mn x n bm
Atau bisa dinyatakan dengan : Ax = B Dimana :
a11 a A = 21 ... a m1
a12 a 22 ... am2
... a1n x1 b1 b ... a 2 n x2 , x= , dan B = 2 ... ... ... ... ... a mn xn b m
Matrik A disebut matrik koefisien atau matrik Jacobian. Vektor x disebut sebagai vektor variabel dan vektor B disebut sebagai vektor Konstanta. Augmented Matrix (matrik perluasan) adalah perluasan dari matrik koefisien A dengan menambahkan vektor B pada kolom terakhir dari matrik A. Augmented (A) = [A B] Sehingga Augmented Matrik dari persamaan linear simultan (5.2) dapat dituliskan :
a11 a 21 ... a m1
a12 a 22 ... am2
... a1n ... a 2 n ... ... ... a mn
b1 b2 .........................................................................(5.3) ... bm
Contoh permasalahan yang dapat dinyatakan dalam bentuk sistem persamaan linear misalnya : Contoh 1 : Seorang pembuat boneka ingin membuat 2 macam boneka yaitu boneka A dan boneka B. Kedua boneka tersebut dibuat dengan menggunakan 2 macam bahan yaitu kain dan kancing. Andaikan boneka A memerlukan 10 meter kain dan 6 kancing sedangkan boneka B memerlukan 8 meter kain dan 8 kancing maka hitunglah berapa buah boneka A dan boneka B yang dapat dihasilkan jika kita memiliki 82 meter kain dan 62 kancing. Masalah ini dapat dimodelkan dengan menyatakan : x1 = jumlah boneka A © Avid-06
61
Fisika Komputasi
x2 = jumlah boneka B Sehingga untuk setiap bahan (kain dan kancing) dapat dinyatakan persamaan : 10 untuk boneka A + 8 untuk boneka B = 82 --- > kain 6 untuk boneka A + 8 untuk boneka B = 62 --- > kancing Persamaan diatas dapat dinyatakan dengan : 10 x1 + 8 x2 = 82 6 x1 + 8 x2 = 62 Penyelesaian dari persamaan tersebut adalah nilai x1 dan x2 yang memenuhi kedua persamaan diatas. Ada beberapa metode yang dapat dipakai untuk mencari solusi dari sistem persamaan linear. Metode – metode ini akan dijelaskan pada sub bab berikutnya. Teorema 5.1 Suatu sistem persamaan linear mempunyai penyelesaian tunggal bila memenuhi syarat-syarat sebagai berikut : - bentuk sistem persamaan linear bujursangkar, dimana jumlah persamaan linear sama dengan jumlah variabel bebas - sistem persamaan linear bersifat non-homogen, yaitu minimal ada satu nilai vektor konstanta B tidak nol atau ada bn ≠ 0 - determinan dari matrik koefisien A tidak sama dengan nol Untuk menyelesaikan permasalahan-permasalahan sistem persamaan linear dapat dilakukan dengan menggunakan metode-metode analitik seperti pemakaian metode grafis, aturan Crammer, atau invers matrik. Metode-metode tersebut dapat dilakukan dengan mudah apabila jumlah variabel dan jumlah persamaannya di bawah 4, tetapi bila ukurannya besar maka metode-metode diatas menjadi sulit dilakukan, sehingga pemakaian metode numerik menjadi suatu alternatif yang banyak digunakan. Metode numerik yang dapat digunakan untuk menyelesaikan permasalahan sistem persamaan linear antara lain : - metode Eliminasi Gauss - metode Eliminasi Gauss-Jordan - Metode Iterasi Gauss-Seidal 5.2 Metode Eliminasi Gauss Metode Eliminasi Gauss merupakan metode yang dikembangkan dari metode eliminasi, yaitu menghilangkan atau mengurangi jumlah variabel sehingga dapat diperoleh nilai dari suatu variabel. Untuk menggunakan metode Eliminasi Gauss ini, terlebih dahulu bentuk matrik diubah menjadi Augmented Matrik sebagai berikut : © Avid-06
62
Fisika Komputasi
a11 a 21 ... a m1
a12 a 22 ... am2
... a1n ... a 2 n ... ... ... a mn
b1 b2 .........................................................................(5.4) ... bm
Metode Eliminasi Gauss adalah suatu metode dimana Augmented matrik (5.4) pada bagian kiri diubah menjadi matrik segitiga atas atau segitiga bawah dengan menerapkan Operasi Baris Elementer (OBE).
a11 a12 a 21 a 22 a31 a32 ... ... a n1 a n 2
a13 a 23 a33 ... a n3
... ... ... ... ...
a1n a2n a3n ... a nn
b1 b2 b3 ... bn
c11 c12 0 c 22 0 0 ... ... 0 0
c13 c 23 c33 ... 0
... ... ... ... ...
c1n c2n c3n ... a nn
d1 d 2 d3 ... d n
Sehingga penyelesaiannya adalah :
xn
dn c nn
xn 1
1 cn1,n 1
(d n 1 cn 1,n xn )
................................. x2
1 (d 2 c23 x3 c24 x4 ... c2 n xn ) c22
x1
1 (d1 c12 x2 c13 x3 ... c1n xn ) c11
Operasi Baris Elementer adalah (OBE) adalah suatu operasi perubahan nilai elemen matrik berdasarkan barisnya, tanpa mengubah matriknya. OBE pada baris ke i+k dengan dasar baris ke – i dapat ditulis : aik, j aik, j c.a i,j
Dimana c adalah konstanta pengali yang diambil dari perbandingan antara nilai a ik,i dan
a i,i Contoh : Selesaikan sistem persamaan linear berikut ini :
x1 x2 x3 6 x1 2 x2 x3 2 © Avid-06
63
Fisika Komputasi
2 x1 x2 2 x3 10 Jawab : Augmented matrik dari sistem persamaan linear tersebut adalah :
1 1 1 6 1 2 1 2 2 1 2 10 Lakukan OBE sebagai berikut :
1 6 1 1 B2 B2 B1 0 1 2 4 B3 B3 2 B1 0 1 0 2
6 1 1 1 0 1 2 4 B3 B3 B2 0 0 2 6 Sehingga didapat penyelesaian yaitu :
x3
6 3 2
1 x2 (4 (2)(3)) 2 1 1 x1 (6 2 3) 1 1 5.2.1 Algoritma Metode Eliminasi Gauss Adapun algoritma dari metode Eliminasi Gauss adalah : 1.Masukkan matrik A dan vektor B beserta ukurannya 2.Buat Augmented Matrik [A|B], namakan dengan A 3.Untuk baris ke-i dimana i=1 sampai n, perhatikan apakah nilai ai,i sama dengan nol. Jika Ya : Tukarkan baris ke-i dengan baris i+k ≤n, dimana ai+k,i tidak sama dengan nol, bila tidak ada berarti perhitungan tidak bisa dilanjutkan dan proses dihentikan tanpa penyelesaian Jika tidak : lanjutkan 4.Untuk baris ke-j, dimana j = i+1 sampai n, Lakukan Operasi Baris Elementer OBE : - hitung c
a j, i a i, i
-untuk kolom k dimana k = 1 sampai n+1 hitung a j,k a j,k c . a i,k © Avid-06
64
Fisika Komputasi
5.Hitung nilai xi, untuk i = n sampai 1 (bergerak dari baris ke-n sampai baris ke-1)
xi
1 (b i a i, j1 x i 1 a i, j 2 x i 2 ... a i,n x n ) a i,i
dimana nilai i+k ≤ n
5.2.2 Program komputer Metode Eliminasi Gauss Berdasarkan algoritma yang dinyatakan dalam sub bab sebelumnya maka kita dapat menuliskan atau mengimplementasikan algoritma tersebut menjadi kode program (source code) dengan bahasa Visual C++ seperti berikut ini : //*****Metode Eliminasi Gauss - solusi sistem pers. Linear**** #include void main() { double A[10][10]; double B[10],X[10]; double c; int n,i,j,k; cout << " Ukuran Matrik : "; cin >>n; cout<<"Matrik A : "<<endl; for (i=1;i<=n;i++) { for (j=1;j<=n;j++) { cout<<"A["<>A[i][j]; } } cout<<endl; cout<<" Vektor B : "<<endl; for (i=1;i<=n;i++) { cout<<"B["<>B[i]; } cout <<endl; //Augmented Matrik for (i=1;i<=n;i++) { A[i][n+1]=B[i]; } for (i=1;i<=n;i++) { if (A[i][i]==0) { //tukar baris for (k=i+1;k<=n;k++) { } } else { © Avid-06
65
Fisika Komputasi
for (j=i+1;j<=n;j++) { //operasi baris elementer c= A[j][i]/A[i][i]; for (k=1;k<=n+1;k++) { A[j][k]=A[j][k]-c*A[i][k]; } } } } //tampilkan matrik hasil operasi baris elementer for (i=1;i<=n;i++) { for (j=1;j<=n+1;j++) { cout<<"A["<=1;i--) { if (i==n) { X[i]=A[i][n+1]/A[i][i]; } else { X[i]=A[i][n+1]; for (j=i+1;j<=n;j++) { X[i]=X[i]-A[i][j]*X[j]; } X[i]=X[i]/A[i][i]; } } //tampilkan X for (i=1;i<=n;i++) { cout<<"X["<
Cobalah jalankan program tersebut dan kemudian cobalah untuk menyelesaikan masalah sistem persamaan linear berikut ini :
x1 x2 x3 6 x1 2 x2 x3 2 2 x1 x2 2 x3 10 Berikut ini adalah hasil eksekusi program untuk masalah tersebut diatas : Ukuran Matrik : 3
© Avid-06
66
Fisika Komputasi
Matrik A : A[11] = 1 A[12] = 1 A[13] = 1 A[21] = 1 A[22] = 2 A[23] = -1 A[31] = 2 A[32] = 1 A[33] = 2 Vektor B : B[1] = 6 B[2] = 2 B[3] = 10 A[11] = 1 A[12] = 1 A[13] = 1 A[14] = 6 A[21] = 0 A[22] = 1 A[23] = -2 A[24] = -4 A[31] = 0 A[32] = 0 A[33] = -2 A[34] = -6 X[1] = 1 X[2] = 2 X[3] = 3
5.3 Metode Eliminasi Gauss-Jordan Metode Eliminasi Gauss – Jordan hampir mirip dengan metode Gauss, hanya saja pada metode ini, di bagian sebelah kiri dari Augmented matrik diubah menjadi matrik diagonal seperti diperlihatkan berikut ini :
a11 a12 a 21 a 22 a31 a32 ... ... a n1 a n 2
a13 a 23 a33 ... a n3
... ... ... ... ...
a1n a2n a3n ... a nn
b1 b2 b3 ... bn
1 0 0 ... 0
0 1 0 ... 0
0 0 1 ... 0
... ... ... ... ...
0 0 0 ... 1
d1 d 2 d3 ... d n
Penyelesaian dari sistem persamaan linear ini adalah nilai d1,d2,d3,...,dn dan atau : x1 = d1 , x2 = d2 , x3 = d3 , ..., xn = dn Cara yang digunakan untuk mendapatkan matrik diagonal dalam metode Gauss – Jordan adalah sama seperti metode Gauss yaitu dengan menerapkan Operasi Baris Elementer (OBE). Contoh : Selesaikan sistem persamaan linear berikut ini : © Avid-06
67
Fisika Komputasi
x1 x2 3 2 x1 4 x2 8 Jawab : Augmented matrik dari sistem persamaan linear tersebut adalah :
1 1 3 2 4 8 Lakukan OBE sebagai berikut :
1 1 3 B2 B2 2B1 0 2 2 1 1 3 B2 B2 2 0 1 1
B1 B1 B2 1 0 2 0 1 1 Sehingga didapat penyelesaian yaitu :
x1 2 x2 1 5.3.1 Algoritma Metode Eliminasi Gauss – Jordan Dari uraian diatas maka dapat dibuat algoritma dari metode Gauss – Jordan seperti berikut ini : 1.Masukkan matrik A dan vektor B beserta ukurannya 2.Buat Augmented Matrik [A|B], namakan dengan A 3.Untuk baris ke-i dimana i=1 sampai n, perhatikan apakah nilai ai,i sama dengan nol. Jika Ya : Tukarkan baris ke-i dengan baris i+k ≤n, dimana ai+k,i tidak sama dengan nol, bila tidak ada berarti perhitungan tidak bisa dilanjutkan dan proses dihentikan tanpa penyelesaian Jika tidak : lanjutkan Jadikan nilai diagonalnya menjadi satu, dengan cara untuk setiap kolom k dimana k = 1 sampai n+1, hitung ai ,k
ai ,k a i ,i
4.Untuk baris ke-j, dimana j = i+1 sampai n, Lakukan Operasi Baris Elementer OBE : - hitung c a j ,i -untuk kolom k dimana k = 1 sampai n+1 hitung
a j ,k a j ,k c.ai ,k © Avid-06
68
Fisika Komputasi
5.Hitung nilai xi, untuk i = 1 sampai n xi ai ,n1
5.3.2 Program komputer Metode Eliminasi Gauss – Jordan Dari algoritma Metode Gauss – Jordan yang sudah dijabarkan dalam sub bab sebelumnya maka selanjutnya kita dapat mengimplementasikan dalam bentuk program komputer ( Compiler Visual C++) : //Metode Gauss - Jordan : solusi sistem pers. linear #include void main() { double A[10][10]; double B[10],X[10]; double c; int n,i,j,k,brs,kol; cout << " Ukuran Matrik : "; cin >>n; cout<<"Matrik A : "<<endl; for (i=1;i<=n;i++) { for (j=1;j<=n;j++) { cout<<"A["<>A[i][j]; } } cout<<endl; cout<<" Vektor B : "<<endl; for (i=1;i<=n;i++) { cout<<"B["<>B[i]; } cout <<endl; //Augmented Matrik for (i=1;i<=n;i++) { A[i][n+1]=B[i]; } cout<<"Augmented Matrik"<<endl; for (i=1;i<=n;i++) { for (j=1;j<=n+1;j++) { cout<<"A["<
69
Fisika Komputasi
} } else { c=A[i][i]; for (kol=1;kol<=n+1;kol++) //elemen diagonal dibuat 1 { A[i][kol]=A[i][kol]/c; } for (brs=1;brs<=n;brs++) { //operasi baris elementer if (brs!=i) { c=A[brs][i]; for (kol=1;kol<=n+1;kol++) { A[brs][kol]=A[brs][kol]-c*A[i][kol]; } } } } } //tampilkan matrik hasil operasi baris elementer for (i=1;i<=n;i++) { for (j=1;j<=n+1;j++) { cout<<"A["<
Eksekusilah program tersebut dan kemudian cobalah untuk menyelesaikan masalah sistem persamaan linear berikut ini :
x1 x2 3 2 x1 4 x 8 Berikut ini adalah hasil eksekusi program untuk masalah tersebut diatas : Ukuran Matrik : 2 Matrik A : A[11] = 1 A[12] = 1 A[21] = 2
© Avid-06
70
Fisika Komputasi
A[22] = 4 Vektor B : B[1] = 3 B[2] = 8 Augmented Matrik A[11] = 1 A[12] = 1 A[13] = 3 A[21] = 2 A[22] = 4 A[23] = 8 A[11] = 1 A[12] = 0 A[13] = 2 A[21] = 0 A[22] = 1 A[23] = 1 X[1] = 2 X[2] = 1
5.4 Metode Iterasi Gauss – Seidal Berbeda dengan 2 metode sebelumnya yang menggunakan Operasi Baris Elementer (OBE) maka dalam Metode Gauss – Seidal kita menggunakan proses iterasi hingga diperoleh nilai yang konvergen (bertambah kecil sampai lebih kecil dari suatu konstanta ). Bila diberikan sistem persamaan linear :
a11x1 a12 x2 a13 x3 ... a1n xn b1 a21x1 a22 x2 a23 x3 ... a2n xn b2 a31x1 a32 x2 a33 x3 ... a3n xn b3
.......................................................... am1 x1 am 2 x2 am3 x3 ... amn xn bm Berilah nilai awal untuk setiap xi ( i = 1 sampai n) kemudian sistem persamaan linear diatas ditulis menjadi : x1
1 (b1 a12 x2 a13 x3 ... a1n xn ) a11
x2
1 (b2 a 21 x1 a 23 x3 ... a 2 n xn ) a 22
...............................................................
xn
1 (bn a n1 x1 a n 2 x2 ... a nn1 xn 1 ) a nn
Dengan menghitung nilai-nilai xi (i = 1 sampai n) dengan memakai persamaan diatas secara terus menerus hingga diperoleh suatu nilai xi yang sama (mendekati) dengan © Avid-06
71
Fisika Komputasi
nilai xi pada iterasi sebelumnya. Atau proses iterasi akan dihentikan jika sudah didapatkan selisih antara 2 nilai iterasi berurutan sudah lebih kecil dari suatu toleransi error yang ditentukan. Contoh : Selesaikan sistem persamaan linear berikut ini :
x1 x2 5 2 x1 4 x2 14 Jawab : Berikan nilai awal : x1 = 0 dan x2 = 0 Susun persamaan menjadi :
x1 5 x2 x2
1 (14 2 x1 ) 4
Iterasi 1 :
x1 5 0 5 x2
1 (14 2.5) 1 4
Iterasi 2 :
x1 5 1 4 x2
1 3 (14 2.4) 4 2
Iterasi 3 :
x1 5 x2
3 7 2 2
1 7 7 (14 2. ) 4 2 4
Iterasi 4 :
x1 5
x2
7 13 4 4
1 13 15 (14 2. ) 4 4 8
Iterasi 5 :
x1 5
© Avid-06
15 25 8 8
72
Fisika Komputasi
x2
1 25 31 (14 2. ) 4 8 16
Iterasi 6 :
x1 5 x2
31 49 16 16
1 49 63 (14 2. ) 4 16 32
Iterasi 7 :
x1 5
x2
63 97 32 32
1 97 127 (14 2. ) 4 32 64
Dapat dilihat nilai x1 dan x2 pada iterasi yang ke-7 sudah tidak berbeda jauh dengan nilai pada iterasi ke-6, sehingga proses iterasi dapat dihentikan dan didapatkan penyelesaian :
x1
97 32
x2
127 64
5.4.1 Algoritma Metode Iterasi Gauss – Seidal Algoritma dari metode Iterasi Gauss – Seidal dapat dinyatakan seperti berikut ini : 1.Masukkan matrik A dan vektor B beserta ukurannya n 2.Tentukan batas maksimum iterasi, max_iter 3.Tentukan toleransi error 4.Tentukan nilai awal dari xi, untuk i = 1 sampai n 5.Simpan xi dalam si, untuk i = 1 sampai n 6.Untuk i = 1 sampai n hitung :
xi
1 (bi ai , j x j ) a i ,i j i
ei xi s i 7.Iterasi <-- iterasi + 1 8.Jika iterasi lebih dari max_iter atau terdapat ei < untuk i = 1 sampai n maka proses dihentikan dan penyelesaiannya adalah xi (i =1 sampai n). Jika tidak , ulangi langkah 5 5.4.2 Program komputer Metode Iterasi Gauss – Seidal Dari algoritma Metode Iterasi Gauss – Seidal maka selanjutnya kita dapat mengimplementasikannnya dalam bentuk program komputer (bahasa Visual C++) : © Avid-06
73
Fisika Komputasi
//Metode Gauss - Seidal : solusi sistem pers. linear #include void main() { double A[10][10]; double B[10],X[10],S[10],err[10]; double eps,sum; bool loop; int n,i,j,max_iterasi,iterasi; cout << " Ukuran Matrik : "; cin >>n; cout<<"Matrik A : "<<endl; for (i=1;i<=n;i++) { for (j=1;j<=n;j++) { cout<<"A["<>A[i][j]; } } cout<<endl; cout<<" Vektor B : "<<endl; for (i=1;i<=n;i++) { cout<<"B["<>B[i]; } cout <<endl; max_iterasi = 10; eps=0.0001; iterasi=1; for (i=1;i<=n;i++) { X[i]=0; S[i]=X[i]; } loop=true; while ((loop==true) && (iterasi<=max_iterasi)) { for (i=1;i<=n;i++) { sum=0; for (j=1;j<=n;j++) { if (j!=i) { sum=sum+A[i][j]*X[j]; } } X[i]=(B[i]-sum)/A[i][i]; err[i]=X[i]-S[i]; if (err[i]<eps) loop==false; } iterasi=iterasi+1; } cout<<"jumlah iterasi = "<
