METODE EULER, EULER-CROMER & VERLET
TUGAS
dibuat untuk memenuhi tugas mata kuliah fisika komputasi pada Program Studi Fisika Fakultas Matematika dan Ilmu Pengetahuan Alam Universitas Pdjadjaran
Muhammad Thariq Fathurrahman 140310140049
UNIVERSITAS PADJADJARAN FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM PROGRAM STUDI FISIKA JATINANGOR 2017
SOAL 1. Latihan membuat coding program dari contoh I-III dan analisa dari hasil plot yang diperoleh dan bandingkan dengan solusi analitis 2. Apa perbedaanya dengan metode/algoritma Euler dengan algoritma Euler-Cromer dan Verlet? 3. Buat program contoh I-III dengan menggunakan algoritma dari soal no. 2 dan bandingkan hasilnya!
Contoh I : Peluruhan Uranium Peluruhan 235U rata-rata lifetime-nya 105 tahun dNu −Nu Persamaan Differensial : dt = τ τ adalah konstanta peluruhan Solusi analitiknya : Nu = Nu(0)e−tτ Nu Metode Euler : Nu(t + Δt) = Nu(t) − τ Δt Parameter : Konstanta peluruhan (τ) : 105 Nilai awal Uranium (Nu) : 235 Nilai t awal (t0) = 0 Jumlah pias (n) = 5x105 Listing Program : %Soal 1 n=1e5; %jumlah pias dibatasi 5*10^5, bila lebih komputer rentan hang dt=1.5; kp=1e5; %konstanta peluruhan Nu=zeros(n,1); %inisialisasi array N uranium Nu t=zeros(n,1); %inisialisasi array waktu t Nu(1)=235; %inisialisasi nilai awal Nu t(1)=0; %inisialisasi nilai awal t %euler for i=1:n-1 Nu(i+1)=Nu(i)-(Nu(i)/kp)*dt; t(i+1)=t(i)+dt; end %perhitungan analitis tn=dt*n; ta=0:dt:tn; %t analitis Nua=Nu(1)*exp(-ta/tau); %solusi Nu analitis plot(t,Nu,'b',ta,Nua,'r--') title('Peluruhan Uranium[235]') xlabel('t (tahun)') ylabel('Nomor Massa U[235]')
Tampilan Program :
Analisa Dengan menggunakan nilai dt optimal yang semakin kecil, diperoleh pendekatan yang sangat mendekati solusi analitis, namun karena keterbatasan pada jumlah datanya, sehingga hanya dapat memodelkan pada rentang waktu yang kecil yang akibatnya plot grafik tak terlihat exponensial. Secara matematis, bentuk solusi analitisnya merupakan exponensial negative, sehingga ketika dt bernilai 1.5, dengan jumlah data yang sama diperoleh grafik eksponensial, jika dilihat dari grafik pada waktu 100 ribu tahun, massa atom terpadu nya menjadi sekitar 86,46 u, namun karena keterbatasan jumlah data maka nilai minimumnya tidak dapat terlihat dan diketahui. Pada kasus ini tak dapat menggunakan metode euler-cromer maupun verlet karena kasusnya merupakan kasus exponensial yang tak berosilasi.
Contoh II : Hambatan Udara Seorang pembalap sepeda bergerak dengan kecepatan v, daya P yang dikeluarkan atlit diperoleh secara eksperimen adalah 400 watt dalam waktu satu jam. Gesekan udara mempengaruhi kecepatan pembalap tersebut Persamaan Differensial : Metode Euler : v(t
dv
P = mv − CρA𝑒 dt m
2
+ Δt) = v(t) + (
Parameter : Daya (P) : 400 Watt Nilai kecepatan awal (v) : 5 m/s Nilai t awal (t0) : 0 Massa pembalap (m) : 75 kg
P mv(t)
−
2
CρA𝑒 (t) m
)Δt
Koefisien drag (C) : 0.5 Massa jenis udara (𝜌) : 1.2
Listing Program : %soal 2 n=5e5; dt=1e-4; p=400; %daya v=zeros(n,1); %inisialisasi array v t=zeros(n,1); %inisialisasi array waktu t v(1)=5; %inisialisasi nilai awal v t(1)=0; %inisialisasi nilai awal t C=0.5; %konstanta drag rho=1.2; %rho udara 1.2 m=75; %massa pembalap + sepeda (bebas) A=1.5; %Luas bidang terkena hambatan udara (bebas) %metode euler for i=1:n-1 v(i+1)=v(i)+((p/(m*v(i)))-((C*rho*A*v(i)*v(i))/m))*dt; t(i+1)=t(i)+dt; end %plot plot(t,v,'b') title('Kecepatan Pembalap Terhadap Waktu') xlabel('t (sekon)') ylabel('Kecepatan V (m/s)')
Tampilan Program:
Analisa Pada kasus pembalap sepeda ini, grafik berikut menggunakan nilai dt yang disesuaikan agar kurvanya dapat terlihat secara menyeluruh membentuk exponensial, dengan parameter standar, dan kecepatan awal 5m/s, pembalap dapat melaju hingga 7.5 m/s dalam waktu 20 detik, namun karena ada hambatan udara, kurvanya menjadi konstan, dengan kata lain hambatan udara menyebabkan kecepatan pembalap menjadi konstan pada titik tertentu. Bila massa diperbesar, waktu yang diperlukan pembalap untuk mencapai kecepatan 7.5 m/s bertambah, begitu pula sebaliknya bila massa diperkecil, waktu yang diperlukan akan semakin sedikit. Kemudian ketika luas permukaan A yang mana terkena hambatan udara diperbesar, maka kecepatan maksimum yang dapat dilakukan oleh pembalap tersebut berkurang, kemudian ketika massa jenis zat yang menghambat di modifikasi menjadi lebih kecil, menyebabkan v maksimum yang dapat dicapai pembalap bertambah. Pada kasus ini tak dapat diaplikasikan metode euler-cromer maupun verlet karena kasusnya merupakan kasus exponensial yang tak berosilasi.
Contoh III: Gerak Pendulum Persamaan Differensial :
dθ dt
= ω dan dω = − gl θ dt g
Solusi analitiknya : 0.2cos(√ l t) g
Metode Euler : ω(t + Δt) = ω(t) − l θ(t)Δt 𝜃(𝑡+Δ𝑡)=𝜃(𝑡)+𝜔(𝑡)Δ𝑡 g
Metode Euler - Cromer : ω(t + Δt) = ω(t) − l θ(t)Δt 𝜃(𝑡+Δ𝑡)=𝜃(𝑡)+𝜔(𝑡+1)Δ𝑡 1
Metode Verlet : θ(t + Δt) = 2 ∙ θ(t) − θ(t − Δt) + θ(t − Δt)(− m)Δt2 ω(t − Δt) = θ(t) − Parameter : Panjang tali (l) : 1 m Percepatan gravitasi (g) : 9.8 m/s2
θ(t−2Δt) 2Δt
Nilai kecepatan sudut awal (ω) : 0 Nilai posisi sudut awal (ɵ) : 0.2 Nilai t awal (t0) : 0
Listing Program : Metode Euler %soal 3 n=5e3; %jumlah pias dibatasi 5*10^3 %dt=1.5e-8; %dt optimal yang menghasilkan error terkecil dt=1e-2; %dt agar rentang dinamis grafik cukup panjang untuk diamati g=9.8; % percepatan gravitasi 9.8 m/s^2 l=1; % panjang tali omega=zeros(n,1); %inisialisasi array omega theta=zeros(n,1); %inisialisasi array theta t=zeros(n,1); %inisialisasi array waktu theta(1)=0.2; %inisialisasi nilai awal theta %metode euler for step=1:n-1 omega(step+1)=omega(step)-(g/l*theta(step))*dt; theta(step+1)=theta(step)+omega(step)*dt; t(step+1)=t(step)+dt; end time=0:1e-2:50; atheta=0.2*cos(sqrt(g/l)*t); %plot plot(t,theta,'b', t, atheta,'-r'); %t = ta jika verlet xlabel('waktu (s)'); ylabel('theta (radian)');
Metode Euler – Cromer %soal 3 n=5e3; %jumlah pias dibatasi 5*10^3 %dt=1.5e-8; %dt optimal yang menghasilkan error terkecil dt=1e-2; %dt agar rentang dinamis grafik cukup panjang untuk diamati g=9.8; % percepatan gravitasi 9.8 m/s^2 l=1; % panjang tali omega=zeros(n,1); %inisialisasi array omega theta=zeros(n,1); %inisialisasi array theta t=zeros(n,1); %inisialisasi array waktu theta(1)=0.2; %inisialisasi nilai awal theta %metode euler-cromer for step=1:n-1 omega(step+1)=omega(step)-(g/l*theta(step))*dt; theta(step+1)=theta(step)+omega(step+1)*dt; t(step+1)=t(step)+dt; end time=0:1e-2:50; atheta=0.2*cos(sqrt(g/l)*t); %plot
plot(t,theta,'b', t, atheta,'-r'); %t = ta jika verlet xlabel('waktu (s)'); ylabel('theta (radian)');
Metode Verlet %soal 3 n=5e3; %jumlah pias dibatasi 5*10^3 dt=1e-2; %dt agar rentang dinamis grafik cukup panjang untuk diamati g=9.8; % percepatan gravitasi 9.8 m/s^2 l=1; % panjang tali omega=zeros(n,1); %inisialisasi array omega theta=zeros(n,1); %inisialisasi array theta ta=zeros(n,1); %inisialisasi array waktu theta(1)=0.2; %inisialisasi nilai awal theta %metode verlet m=1; omega(1)=0.1; theta(2)=theta(1)+omega(1)*dt; ta=0:dt/pi:50; N=length (ta); for step=3:N; theta(step)=2*theta(step-1)-theta(step-2)+dt*dt*(1.0*1/m)*theta(step-1); omega(step-1)=(theta(step)-theta(step-2))/(2*dt); end omega(N)=(theta(N)-theta(N-1))/dt; time=0:dt/pi:50; atheta=0.2*cos(sqrt(g/l)*ta); %plot plot(ta,theta,'b', time, atheta,'-r'); %t = ta jika verlet xlabel('waktu (s)'); ylabel('theta (radian)');
Tampilan Program : Metode Euler
Metode Euler – Cromer
Metode Verlet
Analisa Dengan menggunakan metode euler, kasus gerak harmonic sederhana pendulum tersebut menunjukan kurva sinusoid yang semakin lama, posisinya semakin bertambah bahkan menjadi 10 kali lipat dari posisi simpangan awal 0.2, hal tersebut tak memenuhi kaidah gerak harmonic sederhana. Ketika menggunakan metode euler-cromer, dengan memodifikasi timestep persamaan kedua, dihasilkan grafik sinusoid yang amplitudonya konstan pada 0.2, ketika panjang tali ditambah, frekuensi ayunan akan melambat terlihat dari waktu yang dibutuhkan kurva dari puncak ke puncak yang semakin renggang, sebaliknya jika parameter tali diubah menjadi lebih pendek, frekuensi osilasi tali akan
semakin cepat. Untuk metode verlet hasil yang didapat mendekati terhadap nilai analitiknya yang diperoleh dengan menggunakan perhitungan analitik. Walaupun masih terdapat error yang terlihat dari tidak samanya antara nilai perhitungan analitik dan perhitungan numerik.
Perbedaan metoda/algoritma euler dengan euler-cromer dan verlet a. Metoda Euler Metode Euler merupakan metode numerik paling sederhana dan diturunkan dari deret Taylor, bentuk persamaan pada metoda euler adalah : 𝑦𝑛+1 − 𝑦𝑛 = 𝑓(𝑥𝑛 , 𝑦𝑛 ) 𝑥𝑛+1 − 𝑥𝑛 Kecepatan dalam perhitungan pada metoda ini cukup cepat, tetapi akurasi pada metoda euler limit atau terbatas dan solusi yang dihasilkan tidak stabil, atau kurang presisi. b. Metoda euler-cromer Metoda euler-cromer stabil untuk sistem yang berosilasi. Persamaan pada metoda euler cromer :
𝑣𝑛+1 = 𝑣𝑛 + 𝑎𝑛 ∆𝑡 𝑥𝑛+1 = 𝑥𝑛 + 𝑣𝑛+1 ∆𝑡 Solusi dari metoda euler cromer ini menghemat energi untuk masalah osilasi, tidak seperti metoda euler yang meningkat energi osilator dengan bertambahnya waktu. Sehingga metoda euler cromer sangat stabil pada sistemnya berosilasi. c. Metoda Verlet Metoda ini digunakan dalam simulasi partikel. Sehingga, dengan menggunkan metoda ini dalam mensimulasikan partikel dapat dikerjakan dengan cepat dengan presisi yang tinggi. Pendekatan matematika yang mendekati algoritma asli verlet :
1 𝑥𝑛+1 = 𝑥𝑛 + 𝑣𝑛 ∆𝑡 + 𝑎𝑛 (∆𝑡)2 2 1 𝑣𝑛+1 = 𝑣𝑛 + (𝑎𝑛+1 + 𝑎𝑛 )∆𝑡 2
Berdasarkan perbedaan dari ketiga metoda tadi, dapat dianalisakan dalam bentuk grafik berikut :
Dapat dilihat bahwa, dengan menggunakan metoda euler, perhitungan yang dilakukannya cukup tinggi tetapi tingkat kestabilannya tidak bagus, tetapi dengan metoda euler-cromer, hasilnya cukup stabil dibanding dengan menggunakan metoda euler. Sedangkan metoda verlet, terlihat sangat presisi. Sehingga dapat digunakan seefisien mungkin dibandingkan dengan menggunakan metode sebelumnya.
