INTEGRASI NUMERIK
TEORI DASAR 1.
Metode Gauss Quadrature b
Tinjauan Gauss dalam perhitungan integral I = ∫ f ( x)dx berdasarkan nilai F(x) a
dalam sub-interval yang tidak berjarak sama, melainkan simetris terhadap titik tengah b
interval. Jika I = ∫ f ( x)dx sebagai fungsi integran, maka dengan mengubah variabeL bebas a
x dengan t dalam hubungan x =
(a + b) + (b − a )t , batasan integrasi menjadi t = −1 dan 2
t = 1 . Nilai f (x) menjadi
(a + b) + (b − a )t f ( x) = f = f (t ) …(1) 2 (b − a ) dt , maka I menjadi 2
Karena dx =
(b − a ) f (t )dt 2 −∫1 1
I=
…(2) 1
Dengan menyatakan
∫ f (t )dt sebagai
−1 1
∫ f (t )dt = w
1
f (t1 ) + w2 f (t 2 ) + w3 f (t 3 ) + ... + wn f (t n )
…(3)
−1
Dengan t1 , t 2 , t 3 ,...t n adalah titik dalam interval t=-1 ke t=1 dengan hubungan variabel x dengan t : x1 = b
(a + b) + (b − a ) (a + b) + (b − a) t1 , x 2 = t 2 , dan seterusnya, maka 2 2
∫ f ( x)dx = a
(b − a ) [ w1 f (t1 ) + w2 f (t 2 ) + ... + wn f (t n )] 2
…(4)
Atau persamaan (4) diatas dapat ditulis sebagai
(b − a ) n ∑ wi f (ti ) …(5) 2 i =1
b
∫
f ( x)dx =
a
a) untuk orde 1 = 0.33332 b) untuk orde 2 =0.495447 c) untuk orde 3 = 7.819479 2.
Metode Monte Carlo
Dengan menggunakan metode ini, kita bisa mengatur kehalusan data yang dihasilkan dengan mengubah-ubah jumlah segmentasi yang diinginkan b
I = ∫ f ( x)dx = (b − a ) × f ( x), a
f ( x) = f ( x) f ( x) =
∑ f ( x) n
HASIL PRAKTIKUM Didapatkan hasil integrasi dari fungsi-fungsi berikut : Didapatkan hasil integrasi dari fungsi-fungsi berikut : Waktu Komputasi Waktu Komputasi Metoda GK (t1) Metoda MC (t2) 1 0.000000 0.016 dengan jumlah 2 x dx segmentasi 1 ∫
t1 : t2 0:0.016
0 1 2
∫ ∫ xy dxdy
0.000000
55.578 dengan jumlah segmentasi 1
Nilai Integral Metoda Gauss 2.3333333
Waktu Komputasi Metoda MC (t2) 2.4012
13.999999
13.881316
2
0 1
1
∫x
2
dx
0 1 2
2 ∫ ∫ xy dxdy 0 1
0:55
// Integrasi Numerik //Metode Gauss Quadrature //Integral lipat 1 untuk fungsi f(x)=x^2 terhadap dx #include <stdio.h> #include <stdlib.h> #include <math.h> #include
#include /* Daftar Variable a,b : batas integral n : orde */ int main() { clock_t awal, akhir; float w[10], t[10], a, b, f; int i,n; printf("Fungsi yang akan diintegrasikan adalah f(x)=x*x\n"); printf("Masukkan batas bawah : "); scanf("%f",&a); printf("Masukkan batas atas : "); scanf("%f",&b); printf("Masukan orde n yang akan dipilih 2 atau 3 \n"); scanf("%i",&n); if (n == 2){ w[1]= 1.; w[2]= 1.; t[1]= -0.577350269; t[2]= 0.577350269; f=0.; awal = clock(); for(i=1;i<=2;i++) { f=f+(b-a)/2*w[i]*pow(((a+b)+(b-a)*t[i])/2,2); } akhir = clock(); printf("\n\nHasil Integrasi fungsi tersebut adalah : %f",f); printf("\nWaktu : %f detik\n",(akhir-awal)*1./CLK_TCK); } else if (n == 3){ w[1]= 0.5555555555; w[2]= 0.8888888888; w[3]= 0.5555555555; t[1]= -0.774596669; t[2]= 0.; t[3]= 0.774596669;
f=0.; awal = clock(); for(i=1;i<=3;i++) { f=f+(b-a)/2*w[i]*pow(((a+b)+(b-a)*t[i])/2,2); } akhir = clock(); printf("\n\nHasil Integrasi fungsi tersebut adalah : %f",f); printf("\nWaktu : %f detik\n",(akhir-awal)*1./CLK_TCK); } else { printf(" Maaf saudara tidak memilih 2 ataupun 3 \n"); } getch(); return 0; } // Integrasi Numerik //Metode Gauss Quadrature //Integral lipat 2 untuk fungsi f(x)=x*y*y terhadap dxdy #include <stdio.h> #include <stdlib.h> #include <math.h> #include #include int main() { clock_t awal, akhir; float w[10], t[10], x[10], y[10], a, b, c, d, f; int i, j, n; printf("Fungsi yang akan diintegrasikan adalah f(x)=x*y*y\n"); printf("\nBatas bawah integrasi x: "); scanf("%f",&a); printf("Batas atas integrasi x: "); scanf("%f",&b); printf("Batas bawah integrasi y: "); scanf("%f",&c); printf("Batas atas integrasi y: "); scanf("%f",&d); printf("Masukan orde n yang akan dipilih 2 atau 3 \n"); scanf("%i",&n); if (n == 2){ w[1]= 1.; w[2]= 1.; t[1]= -0.577350269; t[2]= 0.577350269; f=0.;
awal = clock(); for(i=1;i<=2;i++) { x[i]=((a+b)+(b-a)*t[i])/2; for(j=1;j<=2;j++) { y[j]=((c+d)+(d-c)*t[j])/2; f=f+(b-a)*(d-c)/4*w[j]*w[i]*x[i]*pow(y[j],2); } } akhir = clock(); printf("\n\nHasil integrasi fungsi tersebut adalah : %f",f); printf("\nWaktu : %f detik\n",(akhir-awal)*1./CLK_TCK); } else if (n == 3){ w[1]= 0.5555555555; w[2]= 0.8888888888; w[3]= 0.5555555555; t[1]= -0.774596669; t[2]= 0.; t[3]= 0.774596669; f=0.; awal = clock(); for(i=1;i<=3;i++) { x[i]=((a+b)+(b-a)*t[i])/2; for(j=1;j<=3;j++) { y[j]=((c+d)+(d-c)*t[j])/2; f=f+(b-a)*(d-c)/4*w[j]*w[i]*x[i]*pow(y[j],2); } } akhir = clock(); printf("\n\nHasil integrasi fungsi tersebut adalah : %f",f); printf("\nWaktu : %f detik\n",(akhir-awal)*1./CLK_TCK); } else { printf("Anda tidak memilih 2 atau 3\n"); } getch(); return 0; }
//Modul 6 Integrasi Numerik Metode Monte Carlo //Integral lipat untuk fungsi f(x)=x*x terhadap dx #include <stdio.h> #include <stdlib.h> #include #include <math.h> #include int main () { clock_t awal, akhir; printf ("Diketahui Bentuk Fungsi: f(x)=x*x\n\n"); float a,b,Fo,F,dF; int i,n; printf ("Tentukan batas bawah :"); scanf ("%f",&a); printf ("Tentukan batas atas :"); scanf ("%f",&b); printf ("Tentukan jumlah segmentasi awal :"); scanf ("%i",&n); Fo = 0.; awal = clock(); do { F = 0; for (i=1;i<=n;i++) { F = F + pow(a+(rand()%1000/1000.)*(b-a),2.)*(b-a)/n; } dF = fabs(F-Fo); Fo = F; n = n+1; } while (dF>1e-3); akhir = clock(); printf ("\n\nHasil Integral = %f", F); printf("\nWaktu : %f detik\n",(akhir-awal)*1./CLK_TCK); getch(); return 0; } //Modul 6 Integrasi Numerik Metode Monte Carlo //Integral lipat 2 untuk fungsi f(x)=x*y*y terhadap dxdy #include <stdio.h> #include <stdlib.h> #include #include <math.h> #include
int main () { clock_t awal, akhir; printf ("Diketahui Bentuk Fungsi: f(x,y)=x*y*y\n\n"); double a,b,c,d,Fo,F,dF,x,y; int nx, ny, i, j; printf ("Tentukan batas bawah integrasi x :"); scanf ("%lf",&a); printf ("Tentukan Batas atas integrasi x :"); scanf ("%lf",&b); printf ("Tentukan jumlah segmentasi awal x:"); scanf ("%i",&nx); printf ("Tentukan batas bawah integrasi y :"); scanf ("%lf",&c); printf ("Tentukan batas atas integrasi y :"); scanf ("%lf",&d); printf ("Tentukan jumlah segmentasi awal y:"); scanf ("%i",&ny); Fo = 0.; awal = clock(); do { F = 0.; for (i=1;i<=nx;i++) { x = a+(rand()%1000/1000.)*(b-a); for (j=1;j<=ny;j++) { y = c+(rand()%1000/1000.)*(d-c); F = F + x*pow(y,2.)*(b-a)*(d-c)/(nx*ny); printf (" x = %lf y = %lf I = %lf\n", x, y, F); } } dF = fabs (F-Fo); Fo = F; nx = nx+1; ny = ny+1; } while (dF>1e-3); akhir = clock(); printf ("Hasil Integral = %lf", F); printf("Waktu : %f detik\n",(akhir-awal)*1./CLK_TCK); getch (); return 0; }
METODE DEKOMPOSISI MATRIKS LU Contoh aplikasi dari system persamaan linear yaitu sebagai berikut :
Gambar 1. rangkaian dengan 8R dan 4E
Dari hukum Kirchoff I dan II diperoleh : Loop I :
Loop II :
∑ E + ∑ iR = 0
∑ E + ∑ iR = 0
i1 R1 + i3 R3 + i 4 R4 − E1 = 0
i2 R2 + i5 R5 − i3 R3 − E 2 = 0
i1 R1 + i3 R3 + i 4 R4 = E1 ...[1]
i2 R2 + i5 R5 − i3 R3 = E 2 ...(3)
i1 = i2 + i3 ....[ 2]
i2 = i5 + i6 ...[ 4]
Loop III :
Loop IV :
∑ E + ∑ iR = 0
∑ E + ∑ iR = 0
− i4 R4 + i8 R8 + i7 R7 − E3 = 0
i5 R5 + i8 R8 − i6 R6 + E 4 = 0
− i4 R4 + i8 R8 + i7 R7 = E3 ...(5)
− i5 R5 − i8 R8 + i6 R6 = E 4 ...(7)
i1 = i 4 + i7 ...[6]
i7 = i6 + i8 ...(8)
Dari 8 persamaan tersebut, kita dapatkan 4 persamaan linear dengan 8 variabel. Untuk dapat membentuk sebuah matrik bujur sangkar, maka kita harus mereduksi 4 persamaan dengan 8 variabel tersebut menjadi 4 persamaan linear dengan 4 variabel. Sehingga persamaanpersamaan diatas menjadi : Loop I
: i1 R1 + i3 R3 + i 4 R4 = E1
…(9)
Loop II
: i2 R2 − 2i3 R5 + i4 R4 = E 2
…(10)
Loop III
: i1 ( R7 + R8 ) − i 2 R2 − i3 R8 − i4 ( R4 + R7 ) = E3
Loop IV
: − i1 R1 + i 2 ( R6 + R8 ) + i3 ( R5 + R6 + R8 ) + i 4 ( R5 + R6 ) = E 4 …(12)
…(11)
Persamaan-persamaan tersebut dapat dibuat matrik sebagai berikut : 0 R3 R4 R1 i1 E1 i E 0 R2 − 2 R5 R4 2 = 2 ( R7 + R8 ) − R2 − R8 − ( R 4 + R7 ) i 3 E 3 ( R6 + R8 ) ( R5 + R6 + R8 ) ( R5 + R6 ) i4 E 4 − R1
Dan matriks tersebut dapat diselesaikan menggunakan metode dekomposisi LU dengan rumusan sebagai berikut : RI = E R = LU UI = Y
…(13)
LY = E
Untuk menghitung besarnya masing-masing arus, digunakan program devC++ yaitu :
//Sistem Persamaan Linear //Rangkaian dengan 8 buah sumber tegangan (baterai) dan 4 buah tegangan #include<stdio.h> #include<stdlib.h> #include #include<math.h> int main() { float R[8], E[4], A[4][4], U[4][4], L[4][4], I[8], Y[4], D, X, F, G; int i, j, k; FILE *coba; if((coba = fopen("candra.txt", "w")) == NULL){ printf("File tidak dapat dibuka\n"); } printf ("C! \n"); for (i=1;i<=8;i++) { printf ("R%i = ",i); scanf ("%f",&R[i]); } for (i=1;i<=4;i++) { for (j=1;j<=4;j++) { A[i][j] = 0.; U[i][j] = 0.; L[i][j] = 0.; } Y[i]=0.; L[i][i] = 1.; } A[1][1] = R[1] ; A[1][2] = 2.*R[3] A[2][1] = R[2]+R[5] ; A[2][2] = R[3]
; A[1][3] = R[4] ; A[2][3] = 0.
;A[1][4] = 0. ; ;A[2][4] = R[5]
; A[3][1] = -1.*(R[6]+R[7]) ; A[3][2] = -1.*(R[7]+R[6]) ; A[3][3] = -1.*(R[4]+R[6]+R[7]) ;A[3][4] = R[6] ; A[4][1] = -2.*R[5] ; A[4][2] = -1.*R[5] ; A[4][3] = -1.*R[4] ;A[4][4] =R[8] ; for (j=1;j<=4;j++) { U[1][j] = A[1][j]; for (i=2;i<=j;i++) { F = 0.; for (k=1;k<=(i-1);k++)
{ F = F + L[i][k]*U[k][j]; } U[i][j] = A[i][j] - F; } for (i=j+1;i<=4;i++) { G = 0.; for (k=1;k<=(j-1);k++) { G = G + L[i][k]*U[k][j]; } L[i][j] = (A[i][j] - G)/U[j][j]; } } printf ("Silakan masukkan setiap nilai tegangan! \n"); for (i=1;i<=4;i++) { printf ("E%i = ",i); scanf ("%f",&E[i]); } printf("L11=%f\n",L[1][1]); printf("Y1=%f\n",Y[1]); for (i=1;i<=4;i++) { printf ("E%i = %f, ",i,E[i]); } printf("E1=%f\n",E[1]); Y[1]=E[1]; printf("Y1=%f\n",Y[1]); for (i=2;i<=4;i++) { D = 0.; for (j=1;j<=i-1;j++) { D = D + (L[i][j]*Y[j]); } Y[i] = (E[i] - D)/L[i][i]; printf("Y%i=%f",i,Y[i]); } for (i=1;i<=8;i++) { printf("I=%f\n",I[i]); I[i]=0.; printf("I=%f\n",I[i]); } I[4]=Y[4]/U[4][4];
printf("I=%f",I[4]); for (i=3;i>=1;i--) { X=0.; for (j=i+1;j<=4;j++) { X = X + (U[i][j]*I[j]); } I[i] = (Y[i] - X)/U[i][i]; } I[5] = I[1] - I[2]; I[6] = I[2] - I[3]; I[7] = I[4] - I[3]; I[8] = I[1] - I[4]; fprintf(coba,"Konfigurasi R\n"); for (i=1;i<=8;i++) { fprintf (coba,"R%i = %f, ",i,R[i]); } fprintf(coba,"\n\n"); fprintf(coba,"Konfigurasi E\n"); for (i=1;i<=4;i++) { fprintf (coba,"E%i = %f, ",i,E[i]); } fprintf(coba,"\n\n"); printf ("\nMaka Matriks A: \n"); fprintf (coba,"\nMaka Matriks A: \n"); for (i=1;i<=4;i++) { for (j=1;j<=4;j++) { printf ("%f",A[i][j]); fprintf (coba,"%f ",A[i][j]); } printf ("\n"); fprintf (coba,"\n"); } printf ("\nMaka Matriks L: \n"); fprintf (coba,"\nMaka Matriks L: \n"); for (i=1;i<=4;i++) { for (j=1;j<=4;j++) { printf ("%f ",L[i][j]); fprintf (coba,"%f ",L[i][j]);
} printf ("\n"); fprintf (coba,"\n"); } printf ("\nMaka Matriks U: \n"); fprintf (coba,"\nMaka Matriks U: \n"); for (i=1;i<=4;i++) { for (j=1;j<=4;j++) { printf ("%f ",U[i][j]); fprintf (coba,"%f ",U[i][j]); } printf ("\n"); fprintf (coba,"\n"); } printf ("\nSehingga Besar arus: \n"); fprintf (coba,"\nSehingga Besar arus: \n"); for (i=1;i<=8;i++) { printf ("I%i = %f \n",i,I[i]); fprintf (coba,"I%i = %f \n",i,I[i]); } getch (); return 0; } Sehingga didapatkan data-data sebagai berikut untuk besar hambatan yang berbeda-beda.
Sehingga didapatkan data-data sebagai berikut untuk besar hambatan yang berbeda-beda. R1
R2
R3
R4
R5
R6
R7
R8
E1
2
3
4
5
6
7
8
1.5 3
4.5 6
1
10
86
12
33
77
49
26
7.5 9
12
156
812 110
999
666
700
456 146
15
1200 900 2437 1710 1010 2510 729 1000 21
I1
I2
I3
I4
-0.210961
0.410900
-0.188610
-0.134765
0.029248
-0.008727 -0.016842
E2 E3
E4 15
12 15
9
18 2
25
I5
I6
I7
I8
0.648806 -0.621862
0.599510
0.837416 0.859767
0.217011
0.372108 -0.164013
-0.187763
0.155098 -0.506873
-0.01761
0.020256
0.040293 0.008884
-0.037867
0.020037 -0.049020
-0.012048
0.058439
0.078742 -0.004795
-0.070487
0.020302 -0.095584
PERSOALAN NILAI EIGEN Pada pengerjaan ini saya menggunakan nilai dari A adalah 2 sehingga vektornya menjadi:
0 0 A/ 2 − A/ 2 A 0 0 − A/ 2 A 0 0 0 0 1 0 0 −1 2 0 2 −1 0 0 0
0 0 0 1
0 0 0 A / 2
1 /ψ 2 /ψ = 3 /ψ 4 /ψ
1 /ψ 2 /ψ = E 3 /ψ 4 /ψ
E
1 /ψ 2 /ψ 3 /ψ 4 /ψ
1 /ψ 2 /ψ 3 /ψ 4 /ψ
Secara analitik kita mendapatkan nilai eigen dari persamaan di atas adalah dengan cara menghitung determinan matriks berikut: 1− E
0
0
0
0
−1− E
2
0
0
2
−1 − E
0
0
0
0
1− E
=0
Determinan dari matriks ini akan didapat:
Det = ( E − 1) 3 ( E + 3) ( E − 1) 3 ( E + 3) = 0
sehingga E1 = 1 E 2 = −3 Metode Shifted-Inverse Power Method ini hanya bisa menghitung nilai eigen yang berbeda. Jika nilai eigen ada yang sama maka kita harus menggunakan cara yang lain. Secara analitik maka kita akan dapatkan nilai vektor eigen adalah Untuk nilai eigen 1: Vektor eigen = 1 | + z >1 | + z > 2 +1 | − z >1 | + z > 2 +1 | + z >1 | − z > 2 +1 | − z >1 | − z > 2 X1 = (1) 2 + (1) 2 + (1) 2 + (1) 2
Untuk nilai eigen 3: Vektor eigen =
X1 =
1 | − z >1 | + z > 2 −1 | + z >1 | − z > 2 (1) 2 + (1) 2
Dengan menggunakan program dan nilai tebakan X mula-mula adalah [1 1 1 1]. Kita dapat memperoleh nilai eigen yang cocok dengan teori asalah kita menebak nilai alfa yang cukup dekat dengan nilai eigen tersebut (jangan tebak sama dengan nilai eigennya, hal ini akan membuat perhitungan error). Misalkan kita tebak alfa= 1.1 maka kita akan dapatkan nilai eigen 1 sedangkan kita tebak alfa =-3.1 maka kita akan dapatkan nilai eigen -3.
Vektor eigen yang didapat dari perhitungan program adalah sama dengan teori yaitu Untuk nilai eigen 1 X1 =
1 | + z >1 | + z > 2 +1 | − z >1 | + z > 2 +1 | + z >1 | − z > 2 +1 | − z >1 | − z > 2 (1) 2 + (1) 2 + (1) 2 + (1) 2
Untuk nilai eigen -3 X1 =
1 | − z >1 | + z > 2 −1 | + z >1 | − z > 2 (1) 2 + (1) 2
Sebenarnya perhitungan dari computer tidaklah memberikan nilai yang selalu tepat 1. terkadang nilainya adalah 0.99999, namun kita pahami bahwa dalam numeric kita melakukan pendekatan, jadi tidak eksak lagi 100%.
Dari semua hal ini dapat disimpulakan bahwa perhitungan nilai dan vektor eigen yang dilakukan dengan computer adalah benar karena cocok dengan hasil teori. #include<stdio.h> #include<math.h> #include<stdlib.h> #include int main()
{ int i,j,k,l; float A[5][5],X0[5],Y0[5],X1[5],Y1[5],alfa,B[5][5],L[5][5],U[5][5],Y[5],C1,C0,delta,lamda; for(i=1;i<=4;i++){ X0[i]=0; Y0[i]=0; X1[i]=0; Y1[i]=0; C1=0; Y[i]=0; C0=1000; for(j=1;j<=4;j++){ A[i][j]=0; U[i][j]=0; B[i][j]=0; L[i][j]=0; }} A[1][1]=1; A[2][2]=-1; A[3][3]=-1; A[4][4]=1; A[2][3]=2; A[3][2]=2; printf("Tebak nilai alfa: "); scanf("%f", &alfa); for(i=1;i<=4;i++){ printf("Tebak nilai X0 ke-%i: ", i); scanf("%f", &X0[i]); }
for(i=1;i<=4;i++){ L[i][i]=1; for(j=1;j<=4;j++){ if(j==i){ B[i][j]=A[i][j]-alfa; } else{ B[i][j]=A[i][j]; }}} /*for(i=1;i<=4;i++){ for(j=1;j<=4;j++){ printf("%f ", B[i][j]); } printf("\n"); } */ printf("\n"); do{ // gunakan metode dekomposisi LU for(j=1;j<=4;j++){ U[1][j]=B[1][j]; for(i=2;i<=j;i++){
U[i][j]= B[i][j]; for(k=1;k<=(i-1);k++){ U[i][j]= U[i][j] - L[i][k]*U[k][j]; } } for(i=(j+1);i<=4;i++){ L[i][j]=B[i][j]; for(k=1;k<=(j-1);k++){ L[i][j]=L[i][j]-L[i][k]*U[k][j]; } L[i][j]=L[i][j]/U[j][j]; } } /* for(i=1;i<=4;i++){ for(j=1;j<=4;j++){ printf("%f ", U[i][j]); } printf("\n"); } printf("\n"); for(i=1;i<=4;i++){ for(j=1;j<=4;j++){ printf("%f ", L[i][j]); } printf("\n"); } printf("\n");*/ for(i=1;i<=4;i++){ Y[i]=X0[i]; for(j=1;j<=(i-1);j++){ Y[i]=Y[i]-L[i][j]*Y[j]; } Y[i]=Y[i]/L[i][i]; } for(i=4;i>=1;i--){ Y0[i]=Y[i];
for(j=4;j>=(i+1);j--){ Y0[i]=Y0[i]-U[i][j]*Y0[j]; } Y0[i]=Y0[i]/U[i][i]; } C1=Y0[1]; for(i=2;i<=4;i++){ if(sqrt(C1*C1)<sqrt(Y0[i]*Y0[i])){ C1=Y0[i]; }} for(i=1;i<=4;i++){ X1[i]=Y0[i]/C1;} for(i=1;i<=4;i++){ X0[i]=X1[i];} delta=sqrt((C0-C1)*(C0-C1)); C0=C1; }while(delta>0.000001); lamda= 1./C1 + alfa; printf("Nilai eigen: %f\n", lamda); for(i=1;i<=4;i++){ printf("Vektor eigen yang ke-%i: %f\n", i,X1[i]); } getch(); return 0; }
Persamaan Differensial Biasa : Shooting Method
DATA Shooting Method merupakan metode numerik untuk menerapkan syarat betas pada solusi persamaan differensial biasa. Diberikan contoh soal sebagai berikut : Benda bermassa 1 kg diikatkan pada pegas tak bermassa secara horizontal dengan konstanta pegas 50 N/m, pegas tersebut diikatkan pada sebuah dinding kemudian system benda pegas disimpangkan sejauh 10 cm, jika tidak ada gesekan antara benda dengan lantai tentukan berapa sajakah kelajuan awal yang mungkin diberikan agar pada saat 10 detik dari simpangan awal tersebut benda berada pada titik kesetimbangan.
Metode ini memilih harga v1 awal dan v2 awal. Kemudian solusi y1 dan y2 dicari dari persamaan euler. Sehingga secara analitik dapat dianalisis sebagai berikut : y” = f (x,y,y’) dimana y (x1) = y1 dan y(x2) = y2 Persamaan gerak diturunkan 2 kali terhadap waktu t, yaitu : f1 (t , v, y ) =
dy = v Dan dt
f 2 (t , v, y ) =
d2y k =− y 2 dt m
Dengan k = 50 N/m dan m = 1 kg, maka didapatkan persamaan Euler (untuk orde pertama) yaitu, vi +1 = vi + f 2 (t , v, y )∆t = vi − 50 yi ∆t
dan, yi +1 = yi + f1 (t , v, y )∆t = yi + vi ∆t
Simpangan awal pada saat t = 0, adalah y = 0.1 dan saat t = 10, y = 0. kecepatan awal dalam persoalan kali ini ditebak jadi setelah mendapatkan nilainya digunakanlah persamaaan dibawah untuk mengetahui kecepatan awalnya :
v0 = v1 −
v2 − v1 y1 y 2 − y1
Setelah didapatkan kecepatan awalnya, dengan metode euler dihitung simpangannya untuk setiap dx tertentu dengan progam yang telah dibuat. Dari penghitungan komputasi yang dilakukan didapatkan bahwa kecepatan awal yang diperlukan agar pada saat t = 10 simpangan y = 0, yaitu v = Nilai v = 0.052989959717, v = 0.0201873, v = 0.013320922, v = 0.01264953 m/s. Nilai kecepatan awal ini selalu didapatkan untuk setiap tebakan dua kecepatan yang berbeda tanda karena program dibuat seperti itu, dengan nilai dt yang berbeda-beda didapatkan nilai kecepatan yang berbeda pula dan terdapat 4 nilai kecepatan yang memenuhi syarat , Dari penghitungan secara komputasi menggunakan didapatkan plot grafik simpangan terhadap waktu sebagai berikut:
Grafik t terhadap y dengan v = 0.0529899
Simpangan Pegas (m)
0.15 0.1 0.05 0 -0.05
0
2
4
6
8
10
12
-0.1 -0.15 waktu (s)
Pada percobaan kali ini, diberi suatu permasalahan nilai eigen yang digunakan untuk menyelesaikan persamaan diferensial. Dari grafik ini dapat dilihat bahwa pada saat waktu adalah 10 detik maka simpangannya adalah 0. Ini sesuai dengan syarat batas yang diberikan. Satu hal yang juga membuktikan bahwa metode komputasi yang digunakan
benar adalah dari hasilnya yaitu grafik simpanganya sinusoidal. Hasil ini sesuai dengan teori bahwa persamaan untuk benda adalah persamaan osilasi. // Shooting methode #include <stdio.h> #include <stdlib.h> #include #include <math.h> double f(double v) { double y,t,dt,k ; //masukan variable yang diketahui //dengan simpangan pegas = 10cm = 0.1 m y = 0.1; t = 0; dt = 0.00001; k = 50; do { v = v - k*y*dt; y = y + v*dt; t = t + dt; } while(t<9.999); return y; } int main() { double v1,v2,y1,y2; coba: printf ("Silakan Tebak dua nilai kecepatan awal yang berlainan tanda!\n"); printf ("Masukkan v1 = "); scanf ("%lf",&v1); printf ("Masukkan v2 = "); scanf ("%lf",&v2); y1 = f(v1); y2 = f(v2); if(y1*y2>0) { printf ("Maaf anda coba lagi!\n"); goto coba; } do { y1 = f(v1);
y2 = f(v2); printf ("y1=%lf\tv1=%lf\ty2=%lf\tv2=%lf\n",y1,v1,y2,v2); if(fabs(y2)>fabs(y1)) { v2 = (v1+v2)/2; } else { v1 = (v1+v2)/2; } } while(fabs(y2)>1e-6&&fabs(y1)>1e-6); printf ("Maka didapat hasil:\ny1=%5.12lf\t\tv1=%5.12lf\ny2=%5.12lf\t\tv2=%5.12lf",y1,v1,y2,v2); getch (); return 0; }
INTERPOLASI DATA Pada percobaan kali ini mengenai Interpolasi yang memiliki pengertian yang lebih luas merupakan upaya mendefinisikan suatu fungsi dekatan suatu fungsi analitik yang tidak diketahui atau pengganti fungsi rumit yang tak mungkin diperoleh persamaan analitiknya dengan masalah umum interpolasi ini yaitu menjabarkan data-data fungsi dalam fungsi dekatan. Berikut ini adalah tabel data dari hasil eksperimen : NO
Besaran Fisis Besaran Fisis X Y 1 1 2.5 2 3 2.5 3 4 3.0 4 6 8.0 5 9 12.0 6 12 16.0 7 20 18.0 8 25 17.5 9 30 15.5 10 40 8.5 11 50 1.5 Maka dari bantuan Microsoft excel didapatkan grafik sebagi berikut :
Besaran Fisis Y
Grafik Besaran Fisis X terhadap Y (Analitik) 20.0 18.0 16.0 14.0 12.0 10.0 8.0 6.0 4.0 2.0 0.0 0
5
10
15
20
25
30
35
40
45
50
Besaran Fisis X Grafik 1. Grafik data hasil analitik Setelah dibuat dari grafik tersebut, Langkah-langkah selanjutnya dalam melakukan interpolasi pada data adalah membagi data-data sesuai dengan karakteristiknya untuk didekati dengan interpolasi tertentu. Pada praktikum disini kita membagi data-data menjadi empat bagian yaitu untuk range x dari 1 sampai 3, dilakukan interpolasi linier. Untuk range x dari 4 sampai 9, dilakukan interpolasi Lagrange dengan polinomial tingkat 2. Begitu juga dengan data x dari 10 sampai 30 menggunakan interpolasi Lagrange polinomial tingkat 2. Kemudian, data x dari 30 sampai 50 menggunakan interpolasi Linier. Hal tersebut untuk mendefinisikan fungsi dekatan. Maka setelah langkah tersebut dilakukan pemograman menggunakan dev C++ yaitu sebagai berikut : #include <stdio.h> #include <stdlib.h> #include <math.h> #include int main() { int i ; float x1[100], x2[100], y1[100], y2[100], c0, c1, c2, L0, L1, L2;
55
FILE *hasil_candra; if((hasil_candra = fopen("hasil.txt", "w")) == NULL){ printf("file tidak dapat dibuka\n"); } //Dari hasil grafik kita menggunakan 4 range data yaitu sebagai berikut ini : //1.grafik linier dengan Range data (1 - 3) y1[1] = 2.5;
x1[1] = 1.;
y1[3] = 2.5;
x1[3] = 3.;
c0 = ((y1[1]*x1[3])-(y1[3]*x1[1]))/(x1[3]-x1[1]); c1 = (y1[3]-y1[1])/(x1[3]-x1[1]); for(i=1; i<=3; i++){ y1[i] = c0 + c1*i; printf("\nx=%i y=%f", i, y1[i]); fprintf(hasil_candra,"\n%f", y1[i]); } //Interpolasi Lagrange dengan n=2 dengan Range data (4 - 9) y1[4] = 3;
x1[4] = 4;
y1[6] = 8;
x1[6] = 6;
y1[9] = 12;
x1[9] = 9;
for(i=4; i<=8; i++){ L0 = ((i-x1[6])*(i-x1[9]))/((x1[4]-x1[6])*(x1[4]-x1[9])); L1 = ((i-x1[4])*(i-x1[9]))/((x1[6]-x1[4])*(x1[6]-x1[9])); L2 = ((i-x1[4])*(i-x1[6]))/((x1[9]-x1[4])*(x1[9]-x1[6])); y1[i] = y1[4]*L0 + y1[6]*L1 + y1[9]*L2; printf("\nx=%i y=%f", i, y1[i]); fprintf(hasil_candra,"\n%f",y1[i]); } //Interpolasi Lagrange dengan n=2 dengan Range data (9 - 30) y1[9] = 12;
x1[9] = 9;
y1[20] = 18;
x1[20] = 20;
y1[30] = 15.5;
x1[30] = 30;
for(i=9; i<=29; i++){ L0 = ((i-x1[20])*(i-x1[30]))/((x1[9]-x1[20])*(x1[9]-x1[30])); L1 = ((i-x1[9])*(i-x1[30]))/((x1[20]-x1[9])*(x1[20]-x1[30])); L2 = ((i-x1[9])*(i-x1[20]))/((x1[30]-x1[9])*(x1[30]-x1[20])); y1[i] = y1[9]*L0 + y1[20]*L1 + y1[30]*L2; printf("\nx=%i y=%f", i, y1[i]); fprintf(hasil_candra,"\n%f", y1[i]); }
//Grafik yang Linier dengan Range data (30 - 50) y1[30] = 15.5;
x1[30] = 30;
y1[50] = 1.5;
x1[50] = 50;
c0 = ((y1[30]*x1[50])-(y1[50]*x1[30]))/(x1[50]-x1[30]); c1 = (y1[50]-y1[30])/(x1[50]-x1[30]); for(i=30; i<=50; i++){ y1[i] = c0 + c1*i; printf("\nx=%i y=%f", i, y1[i]); fprintf(hasil_candra,"\n%f", y1[i]); } getch(); return 0; }
Didapat dari pemograman tersebut hasil sebagai berikut : Besaran Fisis X
Besaran Fisis Y
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
2.500000 2.500000 2.500000 3.000000 5.733334 8.000000 9.800000 11.133334 12.000000 12.924242 13.772728 14.545455 15.242424 15.863637 16.409090 16.878788 17.272728 17.590910 17.833334 18.000000 18.090910 18.106060 18.045454 17.909092 17.696970 17.409092 17.045456
16.606060 16.090910 15.500000 14.800000 14.100000 13.400001 12.700001 12.000000 11.300000 10.600000 9.900001 9.200001 8.500000 7.800001 7.100000 6.400001 5.700001 5.000000 4.300001 3.600001 2.900001 2.200001 1.500001
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
Maka kembali dengan bantuan Microsoft excel kita dapat membuat grafik yaitu :
Besaran Fisis Y
Grafik Besaran Fisis X terhadap Y (Hasil Program) 20.0 18.0 16.0 14.0 12.0 10.0 8.0 6.0 4.0 2.0 0.0 0
5
10
15
20
25
30
35
Besaran Fisis X
Grafik 2. Grafik data hasil pemograman
40
45
50
55
ANALISIS Dari grafik dapat dilihat bahwa grafik data diatas yang didapat dari hasil pemrograman mendekati atau hampir sama dengan hasil analitik. Dari praktikum ini dengan melakukan pendekatan interpolasi pada data-data diatas, maka kita dapat memperoleh nilai dari tiap x yang belum diketahui. Jadi kita dapat melakukan interpolasi pada data dengan menggunakan interpolasi linier, dan Interpolasi lagrange, Namun dalam pemograman, praktikan tidak menggunakan Interpolasi kuadratik disebabkan menghindari rumitnya program eliminasi gauss untuk mencari konstanta c0, c1, dan c2 secara berturut-turut, sehingga resolusi dari hasil pemograman mungkin tidak seakurat menggunakan Interpolasi kuadratik. Namun dalam hal ini praktikan mengetahui bahwasanya kita dapat menggunakan pendekatan dengan Interpolasi Linier dan Interpolasi Lagrange.
PERSAMAAN DIFFERENSIAL BIASA
DATA Metode Euler merupakan metode numerik yang paling sederhana untuk menyelesaikan persamaan differensial. Dengan persamaan rekursifnya : yi +1 = y i + f ( xi + y i )∆x f ( x, y ) =
dy dx
..........................(1)
Pada percobaan ini diberikan kasus sebagai berikut : Seorang penerjun payung memiliki massa efektif 80 kg, apabila gaya tahanan udara yang dialaminya sebanding dengan laju jatuhnya penerjun dengan koefisien gesekan udara b=0.1. Hitunglah berapa lama waktu yang diperlukanya untuk mencapai tanah jika ia melompat dari ketinggian 1000m.
PENGOLAHAN DATA Secara analitik :
∑ F = ma dy d2y =m 2 dt d t 2 d y dy m 2 + b − mg = 0 dt d t 2 d y b dy + − g = 0 ........................................................................................................(2) d 2 t m dt d2y dy + 10 −4 −g =0 2 dt d t
mg − b
Dalam bentuk lain menjadi : dv + 10 − 4 v − g = 0 dt
Untuk Menyelesaikan persoalan differensial diatas maka menggunakan metode Euler yang secara analitik sebagai berikut :
f 1 (x, y ) : f y =
dy d2y ………(3) dan f v = 2 ………(4) ,dan untuk turunan keduanya yaitu dt dt
f 2 ( x, y ) : ff y = f v =
d2y dt 2
dan ff v = −
b mg − bv …..(5) m m
Maka berdasarkan persamaan 1 dan 2 didapatkan nilai ketinggian, waktu dan kecepatan : Misalkan kita menggunakan nilai awal dengan i=0 Ketinggian
∆t 2 y i +1 = y i + f y ( xi , y i , vi )∆t + ff y ( xi , y i , v) 2 2 mg − bv ∆t y1 = y o + v∆t + ( ) m 2 0.2 ∆t 2 y1 = 1000 + v∆t + 9.8 − v 80 2
Dengan i=0,1,2,3,………………… Kecepatan ∆t 2 2 2 mg − bv − b mg − bv ∆t v1 = vo + ( )∆t + ) ( m m 2 m vi +1 = vi + f v ( xi , y i , vi )∆t + ff v ( xi , y i , v)
0.2 0.2 0.2 ∆t v1 = vo + 9.8 − v ∆t − v 9.8 − 80 80 80 2
2
Waktu t i +1 = t i + ∆t t1 = t 0 + ∆t t1 = 0 + ∆t t1 = ∆t dengan ∆t =
t akhir − t awal n
n adalah jumalah partisi Persamaan-persamaan diatas dapat kita masukkan kedalam sebuah program agar takhir dapat diketahui. Programnya adalah sebagai berikut: //Nama : Candra Mecca Sufyana //NIM : 10204039 //Praktikum Fiskom Modul07 Persamaan diferensial biasa #include<stdio.h> #include<math.h> #include<stdlib.h> #include int main() { float m,b,g,v0,v,y0,y,dt,t;
FILE *hasil_039; // biar umum maka h di soal = y aja di pemograman = ketinggian ;makasih mosi-mosi. v0=0; y0=1000; g=9.8; m=80; t=0; hasil_039=fopen("hasil_039.csv","w"); printf("Masukkan selang waktu yang anda inginkan!"); scanf("%f", &dt); printf("Masukkan koefisien gesekan kinetik udara!"); scanf("%f", &b); do{ fprintf(hasil_039,"%f,%f,%f\n", t,v0,y0); // ngitungnya disini ah... v = v0-(g+b*v0/m)*dt; y = y0 + (v0+v)*dt/2; t = t + dt; v0 = v; y0 = y; }while(y0>=0); fprintf(hasil_039,"%f,%f,%f\n", t,v0,y0); fclose(hasil_039); }
Dari pemograman tersebut maka didapatkan nilai kecepatan, ketinggian dan waktu yang dapat dilihat dari microsoft excel dengan dt = 0.01 dan b=0.1, tidak dicantumkan dalam laporan karena ga tau jadi ngeheng. Maka didapatkan grafik sebagai berikut :
Grafik waktu terhadap ketinggian 1200
ketinggian y
1000 800 600 400 200 0 -200 0
2
4
6
8
10
12
14
16
14
16
waktu t
Grafik waktu terhadap kecepatan 0 -20 0
2
4
6
8
10
12
kecepatan
-40 -60 -80 -100 -120 -140 -160 waktu (s)
Dari hasil grafik yang didapat dapat dianalisa bahwa hasil yang didapatkan dari pemograman sesuai dengan persoalan kasus gerak jatuh bebas yang diberikan dalam artian semakin lama waktu maka semakin rendah pula ketinggian, dan hubungan antara waktu dan kecepatan adalah linear walaupun seharusnya kecepatan bernilai positif, dalam hal ini kita tidak mengabsolutkan nilai kecepatan sehingga grafik kecepatan yang didapat negative.
AKAR PERSAMAAN Metode Newton Rapshon dan Iterasi Langkah – langkah : Hal-hal yang diketahui : Benda 1 • Massa : 1 kg • Mengalami gaya pegas dalam arah sumbu y dengan k =1 N/m • Kecepatan awal : v 1 = iˆ + 0.5 ˆj I.
Benda 2 • Massa : 1 kg • Mengalami gaya konstan sebesar 16 N dalam arah x • Kecepatan awal : v 2 = 4 iˆ + 2 ˆj
MENCARI PERSAMAAN GERAK BENDA I DAN BENDA II Benda 1 : d2 y dt 2 2 d y − ky = m 2 dt 2 d y + ω2 y = 0 dt 2
∑
Fy = m
Solusi dari persamaan diferensial tersebut adalah : y = A sin(ωt )
dan x = v 1, x t x =t
Dari hukum kekekalan energi didapat : 1 2 1 2 mv = ky 2 2
Karena m=k=1, maka y=v1,x Sehingga persamaan gerak untuk benda 1 adalah 1 v s 1 = tiˆ + sin( t ) ˆj 2
Benda 2 :
∑F
x
16 = m
=m
d 2x dt 2
2
d x
dt 2 d x 16 − = 0; m = 1 dt 2 m 2
d 2x dt 2
= 16
Maka solusi persamaan tersebut adalah x = 8t 2 + 4t
Dan
y = 2t
Sehingga persamaan gerak untuk benda 2 adalah
v s 2 = ( 8t 2 + 4 t )iˆ + 2tˆj
Dengan mengeliminasi t didapatkan 1 y 1 = sin( x ) 2 1 1 y2 = − + 1 + 2x 2 2
Karena benda tersebut bertemu disatu titik, berarti y1=y2, didapat f ( x ) = sin 2 ( x ) + 2 sin( x ) − 2x
II.
PROGRAM :
//MELIHAT GRAFIK DARI BENDA 1 DAN 2 dan MENENTUKAN TITIK TEMU KEDUA BENDA #include <stdio.h> #include <stdlib.h> #include #include <math.h> //MENENTUKAN TITIK int main () { float x, c, x0, y, y1, y2, dx, a, b, n, e, eps; int i; FILE *f; printf(" grafik lintasan kedua benda!\n"); printf (" Masukkan nilai x awal = "); scanf ("%f",&a); printf (" Masukkan nilai x akhir = "); scanf ("%f",&b); printf ("Silakan masukkan banyaknya titik yang akan digambar!\n"); printf (" n = "); scanf ("%f",&n); f = fopen ("plot39.xls","w"); for (i=0; i<=n; i++) { dx = (b-a)/n; x = a + i*dx; y1 = 0.5*sin(x); y2 = -0.5 + sqrt(0.5*(x+0.5));
printf (" x = %f y1 = %f y2 = %f \n",x,y1,y2); fprintf (f,"%f %f %f \n",x,y1,y2); }; fclose (f); //MENENTUKAN TITIK PERTEMUAN KEDUA BENDA printf ("\nGrafik dapat dilihat di file : grafik39.xls"); printf ("\n\n\n"); printf ("Saatnya menebak pertemuan kedua benda\n"); printf ("Tebak titik pertemuan kedua benda pada sumbu x!\n"); printf (" x = "); scanf ("%f",&c); e = 0.000001; do { x0 = c; c = c - (0.5*pow(sin(c),2.)+sin(c)-c)/(sin(c)*cos(c)+cos(c)-1.); eps = fabs(c-x0); } while (eps>e); y = 0.5*sin(c); printf ("\nmaka titik temu kedua benda tersebut adalah: \n"); printf (" x = %f y = %f \n\n",c,y); getch (); return 0; } III. MENGGAMBAR GRAFIK KEDUA BENDA Untuk melihat r grafik kedua benda, dibutuhkan titik-titik yang melalui persamaan garis diatas. Dan untuk menentukan titik-titknya digunakan komputasi dan dalam hal ini digunakan metode Newton Raphson. Dari program tersebut didapat 100 titik untuk membuat grafik. Kemudian 100 titik tersebut membentuk grafik seperti ini.
f(x)
Grafik 2 benda 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 0 -0.4 -0.6
Series1 Series2 1
2
3
x
Gambar1. Grafik 2 benda
4
IV.
MENCARI TITIK TEMU KEDUA BENDA Didapatkan titik temu kedua benda menurut metode Newton Rapshon adalah x = 1.494129, y = 0498531 Sedangkan berdasarkan didapat : x = 1.48 dan y = 0.49
