Uygulama Notları: 9¶
FİZ219 - Bilgisayar Programlama I | 23/12/2019¶
Fonksiyonlar II: Eğik Atış Problemi
- Problemin Tanımı ve Uçuş Süresinin Hesabı
- Analitik çözüm
- Sayısal çözüm: yarılama metodu (bisection)
- Hazır gelen kök bulma fonksiyonu: roots()
- Konumların hesaplanması
- Fonksiyonlaştırmak
- Bonus: Güzel grafikler, cici grafikler (Ek bilgi, ilgilenirseniz diye)
- "Meydan okumalar" (Challenges)
- İleri düzey fonksiyonlaştırmak
- Korsan gemisini vuran top
Emre S. Tasci emre.tasci@hacettepe.edu.tr
Eğik Atış Problemi¶
Problemin Tanımı ve Uçuş Süresinin Hesabı¶
$v_0$ ilk hızı ve yatayla $\theta$ açısı yapacak şekilde atılan $m$ kütleli bir cismin üzerine düşey yönde, -y yönünde $F_g=mg$ yerçekimi kuvveti etki ederken, yatay yönde etki eden herhangi bir kuvvet bulunmamaktadır.
Bu durumda, $\vec{r(t)} = \vec{r_0} + \vec{v_0} t + \frac{1}{2} \vec{a} t^2$ genel hareket denklemi, x ve y eksenleri yönlerinde bileşenlerine ayrıldığında: \begin{equation*} v_{0,x} = v_0 \cos{\theta}\\ v_{0,y} = v_0 \sin{\theta}\\ x_0 = y_0 = 0\\ a_x = 0\\ a_y = -g\\ \end{equation*}
\begin{align} x(t) &= x_0 + v_{0,x} t + {\tfrac{1}{2}} a_x t^2\\ \rightarrow x(t) &= v_{0,x} t = v_0 \cos{\theta} t \end{align}\begin{align} y(t) &= y_0 + v_{0,y} t + {\tfrac{1}{2}} a_y t^2\\ \rightarrow y(t) &= v_{0,y} t + {\tfrac{1}{2}} (-g) t^2 = v_0 \sin{\theta} t - {\tfrac{1}{2}} g t^2 \end{align}olarak bulunur.
Cismin toplam uçuş zamanı $t_d$, atıldıktan sonra y değerinin 0 olduğu andır. Bu anı bulmak için:$$y(t_d) = v_0 \sin{\theta} t_d - {\tfrac{1}{2}} g t_d^2 = 0$$ denklemini çözeriz.
Problemimizde parametrelerimizi: $v_0 = 30\, \text{m/s}$, $\theta = 60^{o}$, $g = -9.81\, \text{m/s}^2$ olarak alalım.
Analitik çözüm:¶
İkinci dereceden $ax^2 + b x + c = 0$ denkleminin çözümü: $$ x_{1,2} = \frac {-b \pm \sqrt{b^2 - 4 a c}}{2a} $$ şeklinde bulunur. Katsayıları elimizdeki $t_d$ için uygularsak:$$a= - {\tfrac{1}{2}} g,\;b= v_0 \sin{\theta},\;c=0\\ t_{d_{1,2}} = \frac{-v_0 \sin{\theta} \pm \sqrt{(v_0 \sin{\theta})^2 - 4 (- {\tfrac{1}{2}} g) (0)}}{2 (- {\tfrac{1}{2}} g)} = \frac{-v_0 \sin{\theta} \pm v_0 \sin{\theta} }{-g}={0, \frac{2 v_0 \sin{\theta}}{g}}$$
İlk kök ($t_d = 0$) zaten atış anına geldiğinden bariz; bizi ilgilendiren ikinci, $t_d = \frac{2 v_0 \sin{\theta}}{g}$ kökü. $(x,y)$'nin zamana göre değerlerini hesaplarken aralığımızı sıfırdan bu ana kadar alacağız.
Problemimize özel değerleri yerine koyduğumuzda: $$t_d = \frac{2 v_0 \sin{\theta}}{g} = \frac{2\cdot 30 \sin{(60)}}{9.81}\,\text{s}$$
v0 = 30;
theta = 60;
g = 9.81;
t_d = 2*v0*sind(theta)/g
saniye olarak bulunur. Hatta daha da hassas yazdıralım:
printf("Toplam ucus zamani: %9.6f s\n",t_d)
Sayısal çözüm: yarılama metodu (bisection)¶
Sürekli bir f(x) fonksiyonun verilen belirli bir [a,b] aralığında işaret değiştirdiğini biliyorsak, o halde fonksiyonun o aralıkta en az bir kere yatay ekseni kestiğinin, yani y'nin en az bir kere sıfır değerini aldığının, yani en az bir tane kökü olduğunun garantisi vardır. O halde yapabileceğimiz bir şey, a ile b'nin tam ortasındaki c değerinde fonksiyonun işaretine bakmak, eğer f(a) ile f(c) aynı işaretse bu sefer kökü aramamıza c ile b arasında devam etmek (bu sefer c'ye a deyip, b'yi aynı tutup, yeni c'yi bu yeni a (eski c) ile b'nin orta noktası almak); f(a) ile f(c) zıt işaretliyse, arayışımızı a ile c arasına odaklamak suretiyle gerçek köke istediğimiz hassasiyette yaklaşabiliriz.
Önce, y yönündeki hareket denklemini fonksiyon olarak tanımlayalım:
function f=atis_y(v0,theta,t,g=9.81)
f = v0*sind(theta)*t - 0.5*g*t^2;
endfunction
v0 = 30;
theta = 60;
% atıldıktan az sonraki (t=0.01s) yüksekliğe bakalım:
a = 0.01;
printf("atildiktan %6.3f saniye sonraki yukseklik: %10.5f m.\n",a,atis_y(v0,theta,a))
% bir de 10 saniye sonraki yüksekliğe bakalım:
b = 10;
printf("atildiktan %6.3f saniye sonraki yukseklik: %10.5f m.\n",b,atis_y(v0,theta,b))
Görüldüğü üzere, 10 saniye sonra yükseklik -230.7 olarak verildiğinden, 0 ile 10. saniyeler arasında bir yerde cismimiz yatay eksenden geçmiş. Peki ama nerede?
% 0.01. ile 10 saniyenin tam ortasındaki anda ne yükseklikte idi?
c = (a + b) / 2;
printf("atildiktan %6.3f saniye sonraki yukseklik: %10.5f m.\n",c,atis_y(v0,theta,c))
Atıldıktan 5.005 saniye sonra yüksekliği pozitif; demek ki cisim hala yer seviyesinin üzerindeymiş. Bu durumda "yere çarpış", 5.005 ile 10. saniye arasında bir anda yer alıyor. Bu durumda biraz önce a=0.01s, b=10s için yaptığımız hesabı şimdi a=5.005s ile b=10s için yapacağız.. Bunu sistematiğe bağlarsak:
clear;
function f=atis_y(v0,theta,t,g=9.81)
f = v0*sind(theta)*t - 0.5*g*t^2;
endfunction
v0 = 30;
theta = 60;
istenilen_hassasiyet = 0.00001;
a = 0.01;
f_a = atis_y(v0,theta,a);
b = 10;
f_b = atis_y(v0,theta,b);
adim = 1;
while(abs(a-b)>istenilen_hassasiyet)
% arama araliginin boyu |a-b|, istenilen hassasiyetten
% buyuk oldugu surece, yarilamaya devam ediyoruz
c = (a+b)/2;
f_c = atis_y(v0,theta,c);
printf("%2d: %9.6f %10.5f | %9.6f %10.5f | %9.6f %10.5f\n",adim,a,f_a,c,f_c,b,f_b);
if(sign(f_c) == sign(f_a))
% c ile a'nın yükseklikleri aynı işaretli, o zaman
% arada kök olmasının garantisi yok, arayışa c ile b
% arasında devam edelim
a = c;
f_a = f_c;
else
% c ile b'nin yükseklikleri aynı işaretli, o zaman
% arada kök olmasının garantisi yok, arayışa a ile c
% arasında devam edelim
b = c;
f_b = f_c;
endif
adim = adim + 1;
endwhile
Görüldüğü üzere, bu şekilde, gerçek (analitik) değere verilen hassasiyet mertebesinde yaklaşmış olduk.
Hazır gelen kök bulma fonksiyonu: roots()¶
Bir polinom denkleminin köklerini roots() komutu ile hesaplayabiliriz. Fonksiyonun parametreleri en yüksek dereceli terimin katsayısından başlayarak, sabit katsayıya kadar yazılmış katsayılar vektörüdür. Örneğin: $3x^3 -4x +7 = 0$ denkleminin kökünü bulmak için:
roots([3 0 -4 1])
yazarız (burada $x^2$ olmadığı için, onun katsayısını 0 aldığımıza dikkat edin).
Problemimizdeki ikinci dereceden parabol denklemimiz: $$y(t_d) = v_0 \sin{\theta} t_d - {\tfrac{1}{2}} g t_d^2 = 0$$ idi. Katsayıları roots() fonksiyonunda yerine koyacak olursak:
v0 = 30;
theta = 60;
g = 9.81;
roots([-0.5*g v0*sind(theta) 0])
şeklinde doğrudan kökü bulmuş oluruz. (Bir de, polinom olsun olmasın, en uygun metotla kök bulan fzero() fonksiyonu var ama şimdilik onu bir kenarda bırakıyoruz)
Konumların hesaplanması¶
Artık uçuş süresini bildiğimize göre, bu uçuş süresi boyunca cismimizin koordinatlarını hareket denklemlerinden hesaplayabilir, ardından grafiğini de çizdirebiliriz.
Zaman aralığını $[0,t_d]$ arasında, başlangıç için 5 nokta alıp inceleyelim, bunları yaparken, sıfırdan başlayıp, toparlayalım:
clear;
function f=atis_y(v0,theta,t,g=9.81)
f = v0*sind(theta)*t - 0.5*g*t^2;
endfunction
v0 = 30;
theta = 60;
g = 9.81;
t_d = 2*v0*sind(theta)/g
N = 5;
t_ler = linspace(0,t_d,N)
y_konumlari = [];
for t = t_ler
y_konumlari = [y_konumlari, atis_y(v0,theta,t)];
endfor
y_konumlari
plot(t_ler,y_konumlari,"-*")
Pek fena görünmüyor, fakat iki şeyi fark ettiyseniz ne mutlu size: programı yazarken bir önceki uygulama notlarında fonksiyonları yazarken yapmamanızı tavsiye ettiğim iki şeyi bizzat ben yapmış durumdayım. Nedir o iki şey?..
- Fonksiyonu ayrı, kendi dosyasında değil, doğrudan programın içinde tanımladım:
Bunu yapmak biraz mecburiyetten çünkü içeriğini görün istedim. - Oradaki for... döngüsü:
Elimde linspace ile çıkardığım t değerlerini içeren diziyi/listeyi/vektörü doğrudan fonksiyona beslemek varken, gayet masraflı bir şekilde for döngüsüyle tek tek fonksiyona besledim -- bu affedilecek bir şey değil. Ama fonksiyona bakacak olursanız bu haliyle çalışamaz (nedenini görebildiniz mi?.. oradaki t^2 işlemi eleman bazlı çalışmalı, bu yüzden t.^2 olarak yazılmalıydı).
Haydi kodu bir daha, bu sefer daha verimli yazalım, hazır yazarken hesaplanan nokta sayısını da 100 yapalım (hatta oldu olacak, x'in konumlarını da hesaplayalım):
clear;
function f=atis_y(v0,theta,t,g=9.81)
f = v0*sind(theta)*t - 0.5*g*t.^2;
endfunction
function f=atis_x(v0,theta,t,g=9.81)
f = v0*cosd(theta)*t;
endfunction
v0 = 30;
theta = 60;
g = 9.81;
t_d = 2*v0*sind(theta)/g;
N = 100;
t_ler = linspace(0,t_d,N);
y_konumlari = atis_y(v0,theta,t_ler);
x_konumlari = atis_x(v0,theta,t_ler);
plot(t_ler,y_konumlari,"-*")
xlabel("t (s)")
ylabel("y (m)")
title("Zamana Gore Yukseklik Grafigi")
plot(x_konumlari,y_konumlari,"r-*")
xlabel("x (m)")
ylabel("y (m)")
title("Parcacigin Yorunge Grafigi")
Fonksiyonlaştırmak¶
Şimdi öyle bir fonksiyon yazalım ki, ona ilk hızımızı ve atış açımızı (opsiyonel olarak da yerçekimsel ivmeyi) verdiğimizde, bize cismimizin toplamda kaç saniye uçtuğunu, ulaştığı maksimum yüksekliği ve yatay yönde kat ettiği menzili söylesin.
clear;
function [t_d,h_max,R] = ucus_analiz(v0,theta,g=9.81)
v0_y = v0*sind(theta);
t_d = 2*v0_y/g;
t_d_bolu_2 = t_d /2;
h_max = v0_y*t_d_bolu_2 - 0.5*g*t_d_bolu_2^2;
R = v0*cosd(theta) * t_d;
endfunction
ilk_hiz = 30;
aci = 60;
[ucus_zamani,max_yukseklik,menzil] = ucus_analiz(ilk_hiz,aci);
printf("Ucus zamani: %7.4fs | Maksimum Yukseklik: %6.3fm | Menzil: %6.3fm\n",ucus_zamani,max_yukseklik,menzil)
ara vermeden, aynı ilk hız için hangi açıda en uzun menzile ulaşacak bulalım:
ucus_zamanlari = [];
max_yukseklikler = [];
menziller = [];
acilar = 1:90;
for aci = acilar
[ucus_zamani,max_yukseklik,menzil] = ucus_analiz(ilk_hiz,aci);
ucus_zamanlari = [ucus_zamanlari, ucus_zamani];
max_yukseklikler = [max_yukseklikler, max_yukseklik];
menziller = [menziller, menzil];
endfor
sonuclar = [acilar' ucus_zamanlari' max_yukseklikler' menziller']
Görüldüğü üzere, en uzun uçuş süresine ve en yükseğe tam tepeye atıldığında, $\theta=90^{o}$ olduğunda ulaşıyor ($t_d = 6.11621\text{s}$ ve $h_{max} = 45.87156\text{m}$ ile) ama en uzağa 91.74312m ile $\theta=45^{o}$ ile atıldığında ulaşıyor ki, Fizik I dersini düşündüğümüzde bu bir sürpriz olmamalı! ;)
Bir de açılara bağlı olarak grafiklerini çizdirip, tamamlayalım:
plot(acilar,ucus_zamanlari,"-*");
xlabel("Aci (derece)");
ylabel("Ucus zamani (s)");
title("Acilara gore ucus zamanlari");
plot(acilar,max_yukseklikler,"r-o");
xlabel("Aci (derece)");
ylabel("Maksimum yukseklik (m)");
title("Acilara gore maksimum yukseklikler");
plot(acilar,menziller,"k-^");
xlabel("Aci (derece)");
ylabel("Menzil (m)");
title("Acilara gore menziller");
Bonus: Güzel grafikler, cici grafikler (Ek bilgi, ilgilenirseniz diye)¶
Şimdiye kadar grafiklerimizi hep bir grafik bir panelin içinde ya da iki veya daha fazla grafiği üst üste bindirip, yine tek bir panel içinde olacak şekilde çizegeldik.
Fakat, elimizde birden fazla grafik olduğunda bunları derli toplu yerleştirmek için subplot() komutunu kullanıyoruz. Bu komut, panelimizi belirttiğimiz satır ve sütuna bölüyor, sonra soldan sağa, yukarıdan aşağıya giderek bunları numaralandırıyor. Örneğin, subplot(2,3,5) ile kast ettiğimiz panelimizin 2 satır, 3 sütuna bölünüp, alttaki satırın orta sütununu aktif hale getirip, bir sonraki plot() komutumuzun sonucunun oraya çizdirilmesi. 5.'nin orası olduğunu nasıl hesapladım? En üst, en soldakine 1 deyip, soldan sağa, yukarıdan aşağı ilerleyip numaralandırarak:
+-----+-----+-----+ | 1 | 2 | 3 | +-----+-----+-----+ | 4 | 5 | 6 | +-----+-----+-----+
subplot(2,2,1);
plot(acilar,ucus_zamanlari,"-*");
xlabel("Aci (derece)");
ylabel("Ucus zamani (s)");
title("Acilara gore ucus zamanlari");
subplot(2,2,2);
plot(acilar,max_yukseklikler,"r-o");
xlabel("Aci (derece)");
ylabel("Maksimum yukseklik (m)");
title("Acilara gore maksimum yukseklikler");
subplot(2,2,3);
plot(acilar,menziller,"k-^");
xlabel("Aci (derece)");
ylabel("Menzil (m)");
title("Acilara gore menziller");
subplot(2,2,4);
plot(menziller/2,max_yukseklikler,"ms");
xlabel("? 8)");
ylabel("? ;)");
title("Bakalim bunu yorumlayabilecek misiniz?");
Meydan okumalar (Challenges)¶
Kendini zorlamak isteyenlere iki tane meydan okuma. Ödev değildir, kesinlikle sınavınızdan önce ilgilenmenizi tavsiye etmiyorum, örneğin finallerden sonraki ara tatilde canınız sıkılırsa aklınızda olsun. Çözebilirseniz bana e-posta ile gönderebilirsiniz ama ancak bu dönemin notlarını verdikten sonra (ve o zamandan sonra da hatırlatırsanız!) incelerim. (Özetle demek istediğim, ilginizi çekerse yapmaya çalışın fakat yapın ya da yapmayın, ders notunuza doğrudan bir katkısı olmayacak, o nedenle eğer ilk iki denemenizde olmuyorsa, finallerin sonrasına bırakmanızı tavsiye ederim -- her hâlükârda notlar teslim edilene kadar kontrol etmeyeceğim 8)
İleri düzey fonksiyonlaştırmak¶
Yukarıda, "Fonksiyonlaştırmak" başlıklı kısımda, cismi en uzağa taşıyacak açıyı bulmaya çalışırken, açılar üzerinden for... döngüsü ile ilerledik. Halbuki bize yakışan:
[ucus_zamanlari,max_yukseklikler,menziller] = ucus_analiz(ilk_hiz,acilar);
şeklinde, çat! diye, tek bir çağırışta, hiçbir for... döngüsü kullanmadan sonuca ulaşabilmek olmalıydı. Fonksiyonu bu şekilde (ve tabii ki fonksiyonun içinde de for... kullanmadan yazabilir misiniz?
Korsan gemisini vuran top¶
Halliday & Resnick'in "Fiziğin Temelleri" kitabında değişik bir örnek vardır: Sahildeki bir top $v_0 = 82\text{m/s}$ ilk hızla gülle atabiliyor. $R=560\text{m}$ ötedeki bir korsan gemisini vurması için gerekli açıyı sorar. Soru çok ilginç değil ama çözümü sonucunda bir değil, iki doğru açı çıkar.
Bu açıları, analitik yoldan değil de, nümerik yoldan tespit edip, iki atışın yörüngesini (yani x-y grafiklerini) aynı grafik üzerinde çizdirebilir misiniz?