Praktikum Fisika Komputasi

Praktikum Fisika Komputasi PRAKTIKUM MINGGU 1

Deret Taylor Uraian deret Taylor disekitar xo dinyatakan dengan f(x) yaitu :

f ( x)  

f ( n ) ( x0 ) ( x  x0 ) n n!

f ( n ) ( x0 ) 

dn f dx n

dimana, x  x0

f ( 0 ) ( x0 )  f ( x0 )

jika x0 = 0 maka persamaan tersebut menjadi deret Mac Laurin yaitu :

f ( x)  

f ( n ) (0) n x n!

dimana,

f

(n)

dn f (0)  dx n

x 0

f ( 0) (0)  f (0)

Praktikum : 1)Deret MacLaurin dari fungsi sinus yaitu :

sin( x)  x 

x3 x5   ... 3! 5!

Program komputer dari fungsi sinus tersebut jika dinyatakan dalam bahasa Visual C++ adalah :

//********************************************************* // Menghitung fungsi sinus dengan deret MacLaurin // Compiler : Visual C++ //********************************************************* #include #include <math.h> void main() {

1

Praktikum Fisika Komputasi double eps=1e-5; double x,sinx,pem,pen,s,del; int i,j,m,k,tanda; 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; }

2

Praktikum Fisika Komputasi cout << "Jumlah suku = "<<m <<endl; cout << "sin(x)

= " << s <<endl;

}

Latihan : 1. Buatlah program komputer dari deret Mac Laurin untuk fungsi exp(x) :

exp( x)  1  x 

x 2 x3   ... 2! 3!

2. Buatlah program komputer dari deret Mac Laurin untuk fungsi cos(x) :

cos( x)  1 

x2 x4 x6    ... 2! 4! 6!

3

Praktikum Fisika Komputasi PRAKTIKUM MINGGU 2

Akar-akar dari Persamaan 2.1 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 yaitu:

x n 1  x n 

f ( xn ) f ' ( xn )

sedangkan jika kita melibatkan sampai suku dengan turunan kedua maka akan menghasilkan iterasi Newton Rapshon orde-2 :

x n 1  x n 

f ( xn ) f ( x n ) f '' ( x n ) ' f ( xn )  2 f ' ( xn )

Praktikum : 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;

4

Praktikum Fisika Komputasi 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 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

5

Praktikum Fisika Komputasi #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; 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++ // **********************************************************************

6

Praktikum Fisika Komputasi #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 { 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;

7

Praktikum Fisika Komputasi }

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/mole 0 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 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!. 8

Praktikum Fisika Komputasi 2.2 Pembagian Sintetik Penentuan nilai dan turunan dari suatu persamaan (Polinom) dapat ditentukan dengan cara Pembagian Sintetik. Misalkan polinom orde-n dinyatakan dalam bentuk persamaan seperti berikut ini : f ( x)  an x n  an1 x n1  an2 x n2  an2 x n3  ...

misal n=3 maka didapatkan polinom orde-3 : f ( x)  an x n  an1 x n1  an2 x n2  an2 x n3

Untuk mencari nilai fungsi dan turunannya pada suatu nilai x tertentu maka dapat ditentukan dengan Pembagian Sintetik. Adapun caranya adalah :

x

an

an1

an2

a n 3

(bn x) (bn1 x) (bn2 x) bn

bn1

bn2

+

bn3

dimana,

bn a n bn1  an1  bn x bn2  an2  bn1 x bn3  an3  bn2 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 , an1 , an2 , an3 diganti dengan nilai bn, bn1 , bn2 .

x

bn

bn1

bn2

(cn x) (cn1 x) cn

cn1

+

cn 2

dimana,

cn b n cn1  bn1  cn x cn2  bn2  cn1 x

9

Praktikum Fisika Komputasi Nilai dari c n  2 merupakan nilai turunan dari persamaan f(x) pada x yang bersangkutan. Sebagai salah satu contoh misalkan kita memiliki persamaan polinom orde 3 yaitu :

f ( x)  x 3  11x 2  39 x  45 Jika kita ingin mencari nilai persamaan dan turunannya pada x=2 dengan Pembagian Sintetik maka hasilnya adalah : 2

1 -11 2

2

39 -45 -18 42

1

-9

21

1

-9

21

2

–14

-7

7

1

-3

Jadi didapatkan nilai f(x) pada x=2 adalah –3, sedangkan nilai turunannya adalah 7. 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]; 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

10

Praktikum Fisika Komputasi } 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<=0;i--) { if (i!=1) { cout<<setw(15); cout<<setiosflags(ios::left); cout<
11

Praktikum Fisika Komputasi } else { 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++ //****************************************************************************** #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];

12

Praktikum Fisika Komputasi 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; } 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 : 13

Praktikum Fisika Komputasi 1) Buatlah program komputer untuk menyelesaikan masalah gerak peluru. Dimana sebuah peluru ditembakkan dengan kecepatan awal V0 = 100 m/s dan sudut elevasi  = 450. Per cepatan 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  a 2 x 2  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.

14

Praktikum Fisika Komputasi PRAKTIKUM MINGGU 3

IITIN Integral daIn urTuEnG anRAL DAN TURUNAN 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 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 Karena antara x = a dan x = b terdapat n buah pita maka luas seluruhnya menjadi : n

L   f ( xi )x i 1

3.1.2 Metode Trapesium 15

Praktikum Fisika Komputasi 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 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 ( xn1 )] 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)] 2

PRAKTIKUM : 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++ //*******************************************************

16

Praktikum Fisika Komputasi #include #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

17

Praktikum Fisika Komputasi cout<< pita; 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

0

5

x

18

Praktikum Fisika Komputasi 2) Coba buat program komputer untuk menghitung integral dari fungsi : 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.3 Metode Simpson Jika pada metode Persegi panjang dan Trapesium kita menggunakan kurva yang berbentuk garis lurus maka pada metode Simpson kita akan memakai kurva yang berbentuk parabola atau polinom orde dua : f(x) = a2x2 + a1x + a0 Untuk menurunkan rumus integrasi maka dimisalkan xi+1 ditempatkan di x = 0, xi ditempatkan di x = -x dan xi+2 di x = x. Sehingga nilai integral di bawah kurva f(x) yaitu : x

Ai 1 

 (a

2

x 2  a1 x  a0 )dx

 x

selanjutnya dengan penurunan secara matematik akan didapatkan bahwa :

1 Ai 1  x[ f ( xi  4 f ( xi 1 )  f ( xi  2 )] 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 ( xn1 )  f ( xn )] 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 )] PRAKTIKUM : 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 : 19

Praktikum Fisika Komputasi //******************************************************* //Menghitung integral dengan metode Simpson //Compiler : Visual C++ //******************************************************* #include #include <math.h> #include void main() { double eps=1e-5; 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

20

Praktikum Fisika Komputasi 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 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.

21

Praktikum Fisika Komputasi PRAKTIKUM MINGGU 4

Integral dan Turunan (lanjutan) 3.2 TURUNAN 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 1 1 f ( x  h)  f ( x)  hf ' ( x)  h 2 f '' ( x)  h 3 f ''' ( x)  ... 2 6 Jika kita mengambil selisih antara kedua persamaan tersebut maka akan didapatkan :

f ' ( x) 

f ( x  h)  f ( x  h) 1 2 '''  h f ( x)  ... 2h 6

jika kita mengambil 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) 2h

Untuk mendapatkan Turunan Kedua dari fungsi f(x) maka kita menjumlahkan kedua persamaan f(x+h) dan f(x-h). Sehingga didapatkan persamaan yaitu :

f '' ( x) 

f ( x  h)  2 f ( x)  f ( x  h) 1 2 iv  h f ( x)  ... 12 h2

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) h2

Kedua metode untuk mencari turunan diatas disebut dengan metode Beda Sentral. PRAKTIKUM : 1) Membuat program komputer untuk menghitung Turunan Pertama dari fungsi. Dimana fungsi yang dicari turunannya adalah : 22

Praktikum Fisika Komputasi 2

f(x) = x – 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<
23

Praktikum 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 //*****************

24

Praktikum 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;

27

Praktikum Fisika Komputasi 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. 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, y)h  1 f ' ( x, y)h 2 2 persamaan ini dapat dinyatakan dalam bentuk iterasi :

yi 1  yi  f ( x, y)h  1 f ' ( x, y)h 2 2 Persamaan ini disebut sebagai iterasi Euler orde 2. PRAKTIKUM :

28

Praktikum Fisika Komputasi 1) Membuat 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; dx=4; //************************* y=y+fx*delx+pow(delx,2)/2; //iterasi Euler orde 2 x=x+delx; } if (n==1) { cout.width(15); cout<
29

Praktikum Fisika Komputasi } cout<<endl; y2=y; } cout<<"Nilai Fungsi di x= "<<x2<<" cout<<endl;

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.

30

Praktikum Fisika Komputasi PRAKTIKUM MINGGU 6

Persamaan Diferensial Biasa (lanjutan) 4.3 Metode Runge-Kutta Persamaan umum iterasi Runge-Kutta adalah :

y n1  y n  f ( xn , y n , h) dimana

f ( xn , y n , h)  a1k1  a2 k 2  ...  am k m h  xn1  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  pm1h, y n  qm1,1k1  qm1, 2 k 2  ...  qm1,m1k m1

4.3.1 Runge-Kutta orde 1 Jika diambil m = 1 maka akan menghasilkan iterasi Runge-Kutta orde 1 yaitu :

y n1  y n  a1k1 yn1  yn  a1hf ( xn , y n )

1

karena f ( xn , y n , h)  a1k1 dan k1  hf ( xn , y n ) . Dengan mengambil a1  1 maka persamaan ini sama dengan persamaan iterasi Euler orde-1. 4.3.1 Runge-Kutta orde 2 Jika diambil m = 2 maka akan memberikan iterasi Runge-Kutta orde 2 yaitu :

yn1  yn  a1k1  a2 k 2 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 . Dengan melakukan perhitungan matematis maka didapatkan nilai dari konstanta tersebut seperti dijelaskan di bawah ini. a) Metode Heun

31

Praktikum Fisika Komputasi Dalam metode ini diambil nilai dari konstanta :

a1  a2 

1

2

p1  q11  1 Sehingga persamaan iterasi Runge-Kutta orde 2 menjadi :

y n1  y n  1 2 k1  1 2 k 2 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 :

y n1  y n  13 (k1  2k 2 ) 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 :

y n1  y n  k 2 dimana

k1  hf ( xn , y n ) k 2  hf ( xn  1 2 h, y n  1 2 k1 )

PRAKTIKUM : 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 32

Praktikum Fisika Komputasi 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]; } } //for xx=xx+delx; yy=yy+(k[1]+k[2])/2; } //while if (n==1) { cout.width(15); cout<
33

Praktikum Fisika Komputasi else { delt=yy-y2; 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; 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;

34

Praktikum Fisika Komputasi 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<
}

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.

35

Praktikum Fisika Komputasi 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) { x=xx+delx/2; y=yy+k[i]/2; } } //for xx=xx+delx; yy=yy+k[2]; //iterasi runge-kutta orde 2 - Poligon } //while

36

Praktikum Fisika Komputasi 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.2 Runge-Kutta orde 3 Rumus iterasi Runge-Kutta orde 3 adalah :

y n1  y n  16 (k1  4k 2  k 3 ) 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 ) PRAKTIKUM : 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++

37

Praktikum Fisika Komputasi //************************************************************************************** #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): { 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)

38

Praktikum Fisika Komputasi { 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.2 Runge-Kutta orde 4 Rumus iterasi Runge-Kutta orde 4 adalah :

y n1  y n  16 (k1  2k 2  2k 3  k 4 ) 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 ) PRAKTIKUM : 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

39

Praktikum Fisika Komputasi //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) { 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<
40

Praktikum Fisika Komputasi 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.

41