Arnol'd'un serbest düşme problemi ve çözümü

Fiziksel teoriler bizim manavdan aldığımız yemeklik kuru soğana benzer. Katman katmandırlar. Bir sahada teori kurabilecek kadar bilgi birikiminiz varsa ilkin en basit ve sade bir model kurmakla işe başlar ardından bu modeli çözer ve deneysel verilere karşı sınarsınız. Modelin kendi içindeki tutarsızlıkları veyahut deneysel verilerle uyuşmazlıkları zuhur ettikçe bu modeli ya da teoriyi kademe kademe genişletir ve soğanın kabuğundan en orta katmanına doğru bir yolculuğa koyulursunuz.

Rahmetli Rus matematikçi Vladimir Igorevich Arnol'd'un kuşaklar boyu fizik öğrencilerinin üzerine iz bırakan klasik mekanik kitabında buna dair bir örnek var. Arnol'd önce Galileo'nun serbest düşmeyi sabit ivmeli bir vakıa olarak telakki etmesini örnek veriyor ve bu problemin kinematikle çözülebileceğini vurguluyor. Galileo'nun yaklaşımını lise düzeyindeki fizik derslerinden beri biliyoruz. Yerçekimi ivmesinin $g$ olduğu bir ortamda $h$ yüksekliğinden bırakılan bir taş $T:=\sqrt{2h/g}$ sürede yere düşer. Burada kütlenin gözükmemesi serbest düşmenin kütleden bağımsız olduğunu söyler. Atmosfer direncinin olmadığı bir ortamda bir kuş tüyü ile bir top güllesi aynı yükseklikten bırakıldıklarında aynı sürede yere düşer. Bu bilgi kuru soğanın kabuğudur. Arnol'd soğanın kabuğunu açıp içine yol almak için Newton'ın genel çekim formülünden dem vuruyor ve diyor ki yerçekimi ivmesini sabit kabul edemeyiz. Zira Newton formülüne göre çekim kuvveti taşın yüksekliği ile değişiyor. Bu gözlemi bir örnek olarak kitabında vermiş ama serbest düşmenin ne dinamiğini çözmüş ne de bunu okura bir alıştırma olarak bırakmış. Biz bu postada Hoca'nın verdiği örneğin ayrıntılarına gireceğiz.

Soru: Atmosfer sürtünmesini tamamen ihmal edebileceğimiz, kütlesi $M$, yarıçapı $r$ olan ve eylemsiz küresel bir gezegenin yüzeyindeki çekim ivmesi $g$ olsun. Gezegen yüzeyinden $h$ yüksekliğinde kütlesi $m$ olan ve test parçacığı olarak kabul edebileceğimiz bir taş serbest düşmeye bırakılıyor. Taş ne kadar sürede yere düşer?

Çözüm: Gezegen yüzeyindeki çekim ivmesi Newton teorisine göre $g := GM/r^{2}$ formülüyle verilir. Burada $G$ evrensel genel çekim sabitidir. Taşın herhangi bir anda gezegen yüzeyine olan yüksekliği $z(t)$ fonksiyonu ile verilsin. Sorunun kurulumundan $z(0)=h$ ve $\dot{z}(0)=0$ olması gerektiği barizdir. Toplam mekanik enerji ve onun korunumuyla başlayalım. \begin{eqnarray}\nonumber E(t) &=& \frac{1}{2}m \dot{z}^{2}(t) - \frac{GMm}{r+z(t)} \\ \nonumber &=& \frac{1}{2}m \dot{z}^{2}(t) - \frac{mgr^{2}}{r+z(t)} \\ \nonumber &=& E(0) \\ \nonumber &=& -\frac{mgr^{2}}{r+h} \end{eqnarray} Bu denklemin ikinci ve dördüncü satırlarını birbirine eşit kıldığımızda taşın kütlesi $m$ sadeleşecek ve Galileo formülünde olduğu gibi serbest düşmenin kinematiği kütleden bağımsız hale gelecektir. \[ -\frac{gr^{2}}{r+h} = \frac{1}{2} \dot{z}^{2}(t) - \frac{gr^{2}}{r+z(t)} \] Bu denklem ve problemin kurulumu bize konum ve zaman değişkenlerinin birimlerini öneriyor. Konumu $r$ biriminde ölçersek o zaman birimsiz konum değişkeni $x := z/r$ ile, zamanı $ \sqrt{r/g}$ ile ölçersek birimsiz zaman değişkeni de $\tau := t \sqrt{g/r}$ ile verilmelidir. Son olarak hız değişkeni birimsiz niceliklerde $\dot{z} = \frac{dz}{dt} = \frac{r}{\sqrt{r/g}}\frac{dx}{d\tau} = \sqrt{gr} \frac{dx}{d\tau}$ olur. $x_{\circ} := h/r$ tanımıyla eneri korunumu birimsiz niceliklerde parametre kalabalığından kurtularak sadece başlangıç şartını ihtiva eden bir denklem halini alır. \[ \frac{1}{2} \left( \frac{dx}{d\tau} \right)^{2} = \frac{1}{1+x}-\frac{1}{1+x_{\circ}}=\frac{x_{\circ}-x}{(1+x_{\circ})(1+x)} \ \ \ \text{ya da} \ \ \ \sqrt{\frac{1+x}{x_{\circ}-x}}\frac{dx}{d\tau} = -\sqrt{\frac{2}{1+x_{\circ}}} . \] Karekökü alırken $(-)$ işaretini tercih ettik çünkü taş $-x$ yönünde yere düşüyor ve hızı da bu yönde olmalı.

Formel olarak bir boyutlu bu problemin çözümü aşağıdaki gibidir. \[ -\sqrt{\frac{2}{1+x_{\circ}}} \tau = \int \limits_{x_{\circ}}^{x(\tau)} \sqrt{\frac{1+x}{x_{\circ}-x}}dx = \int \limits_{x_{\circ}}^{x(\tau)} \sqrt{\frac{1+x_{\circ}}{x_{\circ}-x} - 1}dx \] Bu integrali hesalayıp hesaplayamadığınıza bakılmaksızın matematikçiler problemin çözümünü bu noktada bitmiş addederler. Lakin biz bu integrali önce cebirsel ve daha sonra trigonometrik fonksiyon dönüşümlerini kullanarak hesaplayacağız. MIT'nin Integration Bee yarışmalarında eleme sorusu olacak derecede çetrefilli bir integral jimnastiği bizi bekliyor.

İlkin $u := \frac{1+x_{\circ}}{x_{\circ}-x} \iff x = x_{\circ} - \frac{1+x_{\circ}}{u}$ ile tanımladığımız koordinat dönüşümünde integralin sınırları $x = x_{\circ} \iff u=\infty$ ve $x=x(\tau) \iff u(\tau) = \frac{1+x_{\circ}}{x_{\circ}-x(\tau)}$ olur. $dx = \frac{1+x_{\circ}}{u^{2}}du$ ilişkisini integralde yerine koyduğumuzda çalışmamız gerekli integral aşağıdaki gibi olur. \[ \frac{\sqrt{2}\tau}{(1+x_{\circ})^{3/2}} = \int\limits_{u(\tau)}^{\infty} \frac{\sqrt{u-1}}{u^{2}}du \] $x \in [0,x_{\circ}]$ olduğundan $u \in \Big[ \frac{1+x_{\circ}}{x_{\circ}},\infty \Big)$ olduğu görülür.

İkinci olarak $u:=1/\cos^{2}(\varphi)$ dönüşümünü tanımlıyoruz ki integranddaki köklü ifadeyi temizleyelim. Buradan $\varphi(\tau) = \arccos(u^{-1/2}(\tau))$ olur. İntegralin sınırları için $u=\infty \iff \varphi = \pi/2$ ve $u=u(\tau) \iff \varphi = \varphi(\tau) := \arccos (u^{-1/2}(\tau))$ konulmalıdır. $du = \frac{2 \sin (\varphi)}{\cos^{3}(\varphi)} d\varphi$ ve $\frac{\sqrt{u-1}}{u^{2}} = \sin(\varphi)\cos^{3}(\varphi)$ eşitliklerini kullandığımızda integral aşağıdaki basit forma gelir. \[ \frac{\sqrt{2}\tau}{(1+x_{\circ})^{3/2}} = \int\limits_{\varphi(\tau)}^{\pi/2} 2\sin^{2}(\varphi) d\varphi = \int\limits_{\varphi(\tau)}^{\pi/2} (1-\cos(2\varphi)) d \varphi = \frac{\pi}{2} - \varphi(\tau) + \sin(\varphi(\tau)) \cos(\varphi(\tau)) \]

İntegral hesabını bitirdiğimize göre değişken dönüşümünü ta en başa sarıp çözümü $x(\tau)$ cinsinden ifade etmemiz gerekiyor. $\cos \varphi = 1/\sqrt{u}$ olacak şekilde bir dik üçgen çizdiğimizde $\omega := \frac{\pi}{2} - \varphi = \arctan \frac{1}{\sqrt{u-1}}$ tümleyen açısını tanımlayabiliriz. Aynı zamanda $\sin \varphi \cos \varphi = \frac{\sqrt{u-1}}{u}$ olur. Böylece çözüm $u(\tau)$ cinsinden ifade edildiğinde \[ \frac{\sqrt{2}\tau}{(1+x_{\circ})^{3/2}} = \frac{\sqrt{u(\tau)-1}}{u(\tau)} + \arctan \frac{1}{\sqrt{u(\tau)-1}} \] formuna gelir. Nihayet $u(\tau) = \frac{1+x_{\circ}}{x_{\circ}-x(\tau)}$ koyup zaman değişkenini sol tarafta tek başına bırakırsak aşağıdaki aşkın denkleme ulaşırız. \[ \boxed{ \tau = \frac{(1+x_{\circ})^{3/2}}{\sqrt{2}} \arctan \left( \sqrt{\frac{x_{\circ}-x(\tau)}{1+x(\tau)}} \right) + \frac{\sqrt{1+x_{\circ}}}{\sqrt{2}} \sqrt{(1+x(\tau))(x_{\circ}-x(\tau))} } \] Bu aşkın denklemi çözüp $x(\tau)$ fonksiyonunu açık bir biçimde yazabileceğimizi sanmıyorum. Öyle yapabilseydik çok daha güzel bir çözüm olurdu. Ama koyunun olmadığı yerde keçi Abdurrahman Çelebiymiş, mecbur kutu içindeki bu formülle idare edeceğiz.

Şimdi taşın yere düşme anı için $x(\tau_{*}) = 0$ tanımını kullanırsak o zaman birimsiz niceliklerde çakılma zamanını aşağıdaki formülle verebiliriz. \[ \tau_{*} = \frac{(1+x_{\circ})^{3/2}}{\sqrt{2}} \arctan (\sqrt{x_{\circ}}) + \sqrt{\frac{x_{\circ}(1+x_{\circ})}{2}} \] Bu sonucu daha primitif olan Galileo formülüyle karşılaştıralım. $T = \sqrt{2h/g}$ birimsiz niceliklerde $\sqrt{2x_{\circ}}$ olur. $\tau_{*}$ formülünde $x_{\circ} \ll 1$ olduğunu kullanıp Taylor serilerinden faydalandığımızda \[ \tau_{*} = \sqrt{2x_{\circ}} + \frac{5\sqrt{2}}{6} x_{\circ}^{3/2} + \mathcal{O}(x_{\circ}^{5/2}) \] hem Galileo formülünü elde ediyor hem de ona gelen düzeltmeleri görüyoruz. İyi bir teorik çalışmada ileri katmanda bulunan sonuçlar gerekli varsayımlar altında daha iptidai katmanlarda bulunan sonuçları göstermelidir.

Yandaki grafikte serbest düşme zamanı yüksekliğe karşı hem Galileo hem de Newton teorisi kapsamlarında karşılaştırmalı olarak gösteriliyor. Galileo teorisi yükseklikle yerçekimi ivmesinin azaldığını öngörmediği için taş bu teori kapsamında daha erken yere düşüyor. Her iki teori de yüksekliğin çok küçük olduğu limitlerde örtüşüyor. Somut bir örnek verelim. Uluslararası Uzay İstasyonu'nun yerden ortalama yüksekliği kabaca 300 km olsun. Dünyanın yarıçapını da yine kabaca 6000 km alalım. O zaman $x_{\circ} \approx 0.05$ olacaktır. İstasyonun dünyaya serbest düştüğü senaryoda Galileo ve Newton teorilerinin sonuçları grafiğe göre bu bölgede neredeyse ayırdedilemez derecede yakındır.

Şimdi günün en önemli sorusu tabii ki Arnol'd'un bu problemi neden çözmeden bıraktığıdır. Çok sık yaptığı gibi "Bu basit problemin çözümünü okura bırakıyorum." dahi dememiş. Bu sorunun cevabını bilmiyorum. Belki integralle uğraşmak istememiştir...

Carl Bender'ın derlediği iki kar küreme problemi

Aralık ayının sonlarındayız. Calculus ve adi diferansiyel denklemler finallerinin eli kulağında. Herkese final sınavlarında başarılar ve yıl başında bol kar yağışı diliyorum. Bu günlerin anlam ve önemini havi, aşağıdaki soruları Bender ve Orszag'ın yazdığı İleri Matematiksel Metodlar kitabının birinci faslında gördüm. Çözümüyle beraber Yerölçüsü Blogu'nda paylaşmak istedim. Erbabı, Bender'a asimptotik analiz ve bilhassa perturbasyon teorisi literatüründen, rahmetli Orszag'a ise spektral yöntemlerdeki çalışmalarından aşinadır.

R. P. Agnew'ın Kar Küreme Problemi: Günlerden bir gün, yoğun ama sabit hızda bir kar yağışı başlamış. Öğlen 12:00'de çalışmaya başlayan bir kar küreme ekibi çalışmalarının ilk saatinde 2 km, ikinci saatinde ise 1 km yol temizleyebilmiş. Buna göre kar yağışı ne zaman başlamıştır? (İpucu: Kar küreme makinesinin hızı, kar yüksekliği ile ters orantılıdır.)

Çözüm: Zamanı saat biriminde ölçelim. Zamanın sıfır noktası öğlen saat 12:00 olarak tespit edilsin. Kar yağışı $\tau \lt 0$ anında başlamış olsun. Kar sabit hızda yağdığına göre, kar yüksekliği ile zaman arasında anlık olarak $h(t) = \alpha (t-\tau)$ şeklinde bir ilişki olmalıdır. Burada $\alpha \gt 0$ olmak üzere bir orantı sabitidir. Kar küreme makinesinin hızı kar yüksekliği ile ters orantılı olduğuna göre $t \gt 0$ için aşağıdaki adi diferansiyel denklem geçerli olur. \[ \frac{dx(t)}{dt} = \frac{\beta}{\alpha} \frac{1}{t - \tau} = \frac{\kappa}{t-\tau} \] Burada $\beta \gt 0$ kar yüksekliği ile kar küreme makinesinin hızı arasında bir (ters) orantı sabiti olup, $0 \lt \kappa := \beta / \alpha$ olarak tanımlanmıştır. $x(0)=0$ başlangıç şartına tabi $x(t)$ fonksiyonu ise kar küreme makinesinin konumunu veriyor. Zannedersem tarif edilebilecek en basit adi diferansiyel denklem şablonuna ait olan bu denklemin çözümü Calculus'un temel teoremiyle kolayca yapılabilir. Sonuç aşağıdadır. \[ x(t) - x(0) = x(t) = \kappa \left( \log(t-\tau) - \log(-\tau) \right) \]

Şimdi soruda verilenleri kullanalım. Çalışmanın ilk saatinde 2 km, ikinci saatinde 1 km yol alınmış. Yani $x(1)-x(0) = 2$ km ve $x(2)-x(1)=1$ km alınacak. (Anlamlı sayılar jandarması, lütfen otur oturduğun yerde!) Veyahut \begin{eqnarray}\nonumber 2 &=& \kappa \left( \log(1-\tau) - \log(-\tau) \right), \\ \nonumber 1 &=& \kappa \left( \log(2-\tau) - \log(1-\tau)\right), \end{eqnarray} olur. $\kappa$ orantı sabitini elemek için bu iki denklem taraf tarafa bölünür ve yeniden düzenlenirse aşağıdaki aşkın (transcendental) denkleme ulaşılır. \[ 2\log(2-\tau) - 3\log(1-\tau) + \log(-\tau) = 0 \] Bu denklem aşkın gibi gözüküyor ancak her iki tarafın üstelini aldığımızda onu cebirsel hale getirmek mümkün. \[ 1 = \frac{(2-\tau)^{2}(-\tau)}{(1-\tau)^{3}} = \frac{-\tau^{3}+4\tau^{2}-4\tau}{1-3\tau+3\tau^{2}-\tau^{3}} \] Gerekli sadeleştirmelerden sonra denklemin kübik değil kuadratik olduğu görülür. \[ \tau^{2} - \tau - 1 = 0 \] Bu denklemin negatif kökünü aldığımızda \[ \tau = \frac{1-\sqrt{5}}{2} \approx -0,6180 \ \text{saat} \ \ \ \text{ya da } \ \ \tau\approx -37,08 \ \text{dakika} \] bulunur. Buna göre kar yağışı saat 11:23 sularında başlamıştır.

---

Yukarıdaki orijinal problemin çözümü nisbeten basitti. Riyazi ve fenni bilimlerdeki her fikri genelleştirmek ve bir kadame daha zorlaştırmak geleneği uyarınca, Bender ve Orszag bu problemin daha karmaşık bir versiyonuna da yer vermiş.

M. S. Klamkin'in Büyük Kar Küreme Konvoyu: (ZOR) Yevmen min-el-eyyam, yine yoğun ama sabit hızda bir kar yağışı başlamıştı. Özdeş üç kar küreme makinesini kullanan üç ekip, aynı noktadan başlayarak, aynı yolu, aynı istikamette temizlemek amacıyla saat 12:00, 13:00 ve 14:00'te harekete geçti. Bu konvoydaki üç ekip de tek bir noktada tekrar karşılaştıklarına göre kar yağışı saat kaçta başlamıştır?

Çözüm: Bir önceki sorunun çözümünde olduğu gibi zamanı yine saat biriminde ölçecek ve zamanın sıfır noktasını saat 12:00 olarak alacağız. Kar yağışı $\tau \lt 0$ anında başlasın. Sırasıyla birinci, ikinci ve üçüncü ekibin kat ettiği yollar $x(t)$, $y(t)$ ve $z(t)$ fonksiyonları ile verilsin. Bu fonksiyonların başlangıç şartları $x(0)=y(1)=z(2)=0$ şeklindedir. Çözümün nihayetinde bizden $T \gt 2$ olmak üzere $x(T)=y(T)=z(T)$ denklemini çözüp $\tau$ değerini bulmamız bekleniyor. Birimsiz niceliklere geçerek sağda solda $\kappa$ niceliğini yazmak yerine $\xi := x / \kappa$, $\eta := y/\kappa$ ve $\zeta := z / \kappa$ ile notasyonu bir kademe daha sadeleştirebiliriz.

İlk ekibin yol fonksiyonu için hareket denklemini ve çözümünü aslında yukarıda verdik. Burada sonucu birimsiz niceliklerde yeniden üretmekle yetineceğiz. \[ \frac{d \xi (t)}{dt} = \frac{1}{t-\tau} \ \ \ \Rightarrow \ \ \ \xi(t) = \log \left( \frac{t-\tau}{-\tau} \right) \ \ \ \text{ve} \ \ \ t = t_{1}(\xi) := -\tau \left( e^{\xi} - 1\right). \]

İkinci ekibin kürediği yol için benzer bir hareket denklemi geçerli ancak ilk ekip yolu onlardan önce temizlediği için, $\eta$ noktasında yolda yeniden biriken kar toplamda $t-\tau$ değil de $t - t_{1}(\eta)$ kadar sürede birikmiş olmalıdır. Bu da çözülmesi gereken hareket denklemini aşağıdaki forma kor. \[ \frac{d\eta (t)}{dt} = \frac{1}{t-t_{1}(\eta)} = \frac{1}{t-\tau+\tau e^{\eta}} \] $(t-\tau)/\tau =: s$ tanımıyla yeni değişkenlerde bu denklem Lambert-$W$ fonksiyonunun türevine dönüşür. Ama daha kolayı var: $t=t_{2}(\eta)$ tanımıyla bu denklem lineer hale gelerek Leibniz'in integral faktörü tekniğiyle kolayca çözülür. \[ \frac{dt}{d\eta} = t - \tau + \tau e^{\eta} \ \ \ \Rightarrow \ \ \ \frac{dt}{d\eta} - t = - \tau + \tau e^{\eta} \ \ \ \Rightarrow \ \ \ \frac{d}{d\eta} \left( t e^{-\eta} \right) = -\tau e^{-\eta} + \tau \] $t_{2}(0)=1$ başlangıç şartı kullanılarak bu toplam diferansiyelin integrali rahatça alınır. Sonuç aşağıdaki gibidir. \[ t = t_{2}(\eta) = \tau + \left( 1 - \tau + \tau \eta \right) e^{\eta} \]

Ara Yorum: Problemin kurucu denklemleri (constitutive equations) göz önüne alındığında birinci ekiple ikinci ekip arasında nitel olarak net bir fark oluşuyor. İlk ekibin önünde hep kar var ve hatta kar kalınlığı giderek artıyor. Bu da ilk ekibin giderek yavaşlamasına neden oluyor. Ancak ikinci ekip, ilk ekibin daha önceden temizlediği yolu tekrar temizlediği için ilk ekibe yetişene kadar kar kalınlığı yol boyunca giderek azalıyor. Dolayısıyla ikinci ekip hızlanarak ilerliyor. Hatta ilk ekibe yetiştiği anda önündeki kar kalınlığı teorik olarak sıfırlandığı için $\dot{\eta}$ denkleminde bir tekillik (singularite, sıfır payda) oluşuyor. Bu singulariteye sonlu zamanda erişiliyor. Öte yandan bizim çözdüğümüz diferansiyel denklemde singularite $dt_{2}/d \eta = 0$ anlamına gelmektedir. Tekilliğin $t^{\star} = t_{2}(\eta^{\star})$ zamanında yakalandığını varsaydığımızda yaptığımız analizin $t_{2} \in (1,t^{\star})$ için geçerli olduğu hatırdan çıkarılmamalıdır.

Şimdi singularitenin tanımından $t^{\star} = -\tau e^{\eta^{\star}} + \tau$ olduğu görülür. Ayrıca bu değer $t_{2}(\eta^{\star})$ değerine eşitlenirse o zaman \[ t^{\star} = -\tau e^{-1/\tau} + \tau \ \ \ \text{ve} \ \ \ \eta^{\star} = -\frac{1}{\tau} \] ara sonuçlarına ulaşılır. Bu aşamada cingözlük edip $t^{\star} = t_{1}(\eta^{\star})$ denklemini çözerek $\tau$ değerine ulaşabileceğinizi sanıyorsanız boşuna denemeyin. $a=a$ şablonunda bir totolojiyle karşılaşıyorsunuz. Çare yok. Mecbur üçüncü ekip için de benzer bir analiz yapmak durumundayız.

Üçüncü ekip $\zeta(2) = 0$ başlangıç şartına tabi olarak yola koyulacak ve $\zeta$ konumuna $t$ zamanında eriştiğinde önündeki kar kalınlığının akümülasyonu için geçmesi gereken zaman $t-t_{2}(\zeta)$ kadar olacak. O zaman çözmemiz gereken hareket denklemi aşağıdaki gibidir. \[ \frac{d \zeta (t)}{dt} = \frac{1}{t-t_{2}(\zeta)} = \frac{1}{t-\tau-(1-\tau+\tau \zeta)e^{\zeta}} \] Bu denklemi çözmek için kendimizi boş yere parelememizin bir manası yok zira $t=t_{3}(\zeta)$ fonksiyonu için geçerli olan diferansiyel denklem de lineer ve integral faktörünü çok kolayca bularak bir toplam diferansiyele ulaşabiliyoruz. \[ \frac{dt}{d\zeta} = t-\tau-(1-\tau+\tau \zeta)e^{\zeta} \ \ \ \Rightarrow \ \ \ \frac{dt}{d\zeta} - t = -\tau-(1-\tau+\tau \zeta)e^{\zeta} \ \ \ \Rightarrow \ \ \ \frac{d}{d\zeta} \left( t e^{-\zeta} \right) = -\tau e^{-\zeta} - 1 + \tau - \tau \zeta \] Calculus'un temel teoremini kullanarak integral alırken dikkat etmemiz gereken yegane püf nokta $t(0)=t_{3}(0)=2$ olduğudur. İntegralin sonucu aşağıdaki gibidir. \[ t = t_{3}(\zeta) = \tau + \left( 2 - \tau - (1-\tau)\zeta - \frac{\tau}{2}\zeta^{2} \right) e^{\zeta} \] Allah'tan başka bir ekip yok da işler büsbütün zıvanadan çıkmıyor.

Nihayet çözüme ulaşacağız! Daha önceden yazdığımız singularite argümanı üçüncü ekip ikinci ekibi yakaladığı anda da geçerli. Ancak bizden istenen üçüncü ekibin ikinci ekibi yakalamasıyla ikinci ekibin ilk ekibi yakalamasının aynı konum ve zamanda gerçekleşmesi. Yani $t^{\star} = t_{3}(\eta^{\star})$ denklemini çözmemiz gerekiyor. $t^{\star}$ ve $\eta^{\star}$ için yukarıda verdiğimiz ifadeleri yerine koyduğumuzda neredeyse her şey birbirini götürüyor ve şaşırtıcı derece kolay birinci dereceden bir cebirsel denklemle $\tau = -1/2$ sonucuna ulaşıyoruz. (Lütfen deneyiniz.) Dakika biriminde bu değer -30 dakika olduğundan kar yağışı tam olarak saat 11:30'da başlamıştır. Ekipler $t^{\star} = \tfrac{1}{2} ( e^{2}-1) \approx 3,1945$ saat sonra ortak tekilliğe ulaşırlar. Geleneksel ifadesiyle bu 15:12 sularına tekabül eder.

Lineer hava sürtünmesinin olduğu bir ortamda yukarı atılan bir top asla atıldığı hızla yere çakılmaz

Slater'ın Mekanik kitabındaki ilk fasılda verilen alıştırma problemlerine göz gezdirirken aşağıdaki soruyu gördüm. Muhtemelen yazar Lambert-$W$ fonksiyonlarını kullanmadan bu problemi çözmüştür.

Soru: $t=0$ anında başlangıç hızı $v > 0$ olacak şekilde bir top yerden yukarı doğru atılsın. Hava sürtünmesi ise hızın tersi yönünde ve hızla orantılı olsun. Top yere $v$ hızından daha düşük bir hızla çakılır.

Çözüm: pdf indir

Yarışan ve tersinmez iki otokatalitik tepkime trişkadan olmayan bir dengede buluşabilir

A+B→2A ve 2A+B→3B otokataliz tepkimelerini ele alalım. Bu tepkimelerin ilkinde A ikinci dereceden, ikincisinde ise B üçüncü dereceden kendi oluşumlarını hızlandırmakta olduklarından bu tip tepkimelere otokatalitik tepkimeler denir. Tersinmez tepkimelerin hız sabitleri $\alpha > 0$ ve $\beta > 0$ olsun. A ve B maddelerinin derişimlerini sırasıyla $a(t) \geq 0$ ve $b(t) \geq 0$ fonksiyonlarıyla temsil edelim. Kütle aksiyon kanunu kabul edersek aşağıdaki hareket denklemlerini hemen yazabiliriz. \begin{eqnarray} \nonumber \dot{a}(t) &=& \alpha a(t)b(t) - 2\beta a^{2}(t) b(t) \\ \nonumber \dot{b}(t) &=& -\alpha a(t)b(t) + 2\beta a^{2}(t) b(t) \end{eqnarray} Burada $\dot{f}(t) := df(t)/dt$ ilgili fonksiyonun zamana göre türevini temsil etmektedir. Dikkat edilirse $\dot{a} + \dot{b} = 0$ olduğu görülür ki bu bize kütlenin korunumu kanunundan başka bir şey vermez. \begin{equation*} a(t) + b(t) = a(0) + b(0) =: m \end{equation*} $m \geq 0$ toplam kütle olup, hız sabitleriyle birlikte problemi tanımlayan parametrelerden biridir.

Klasik mekanikte olduğu gibi her korunum kanunu hareket denklemlerinden birisini elemekte kullanılabilir. $b(t) = m - a(t)$ yazıldığında aslında çözülmesi gerekli sadece bir adet diferansiyel denklem olduğu görülür. \begin{equation*} \dot{a} = \alpha ab - 2\beta a^{2}b = ab(\alpha - 2\beta a) = a(m-a)(\alpha - 2\beta a) \end{equation*} Problemi parametre kalabalığından kurtarmak için sırasıyla birimsiz konsantrasyon ve birimsiz zaman niceliklerini $a =: xm$ ve $\tau := \alpha mt$ tanımlayalım. $x \in [0,1]$ ve $\tau \in [0,\infty)$ olduğu barizdir. Bu tanımlar yerlerine konulduğunda aşağıdaki nispeten daha sade diferansiyle denkleme ulaşıyoruz. \begin{equation*} x^{\prime} = x(1-x)(1 - \kappa x) \end{equation*} Burada $0 \leq \kappa := 2\beta m / \alpha $ bütün parametrelerin öbeklendiği birimsiz bir parametredir ve dinamiğin karakteri bu niceliğin aldığı değerlere göre tayin edilir. (Unutmadan $f^{\prime}(\tau) := df(\tau)/d\tau$ olduğunu belirtelim.)

Çalıştığımız dinamik sistemin denge noktaları $x^{\prime}=0$ denkleminin kökleriyle verilir. Bu kökler toplamda üç adettir: $x_{*} \in \{ 0 , 1 , \kappa^{-1} \}$. $x_{*}=0$ noktasında hiç A maddesi olmadığından, $x_{*} = 1$ noktasında ise hiç B maddesi olmadığından mekanizmadaki reaksiyonlar ilerleyemez ve sistem dengededir. Biz bu tip bariz dengelere trişkadan dengeler diyoruz. Şimdi genelliği biraz kaybedelim ve $\kappa > 1$ olduğunu varsayalım. $\kappa$ değeri hız sabitlerinin dışında toplam madde miktarına da bağlı olduğu için böylesi bir durumu her zaman ayarlayabiliriz. Mekanizmamızdaki her iki tepkime de tersinmez ve böylesi tepkimlerde genellikle sistem reaktantlardan en az biri tükenene değin ilerler ve bütün derişimlerin pozitif olduğu bir denge noktası pek yoktur. Ama burada otokataliz mekanizması böylesi bir dengeyi mümkün kılmaktadır.

$\kappa > 1$ için $x \in (0,\kappa^{-1})$ olsun. O zaman $x^{\prime} > 0$ olacak ve sistem $\tau \to \infty$ limitinde $x \to \kappa^{-1} $ dengesine varacaktır. Benzer şekilde $x \in (\kappa^{-1},1)$ durumunda da $x^{\prime} \lt 0 $ olacağından asimptotik olarak yine aynı dengeye varılır. Toparladığımızda $x_{*} = \kappa^{-1}$ dengesinin kararlı olduğunu da görüyoruz.

İşaret: Dengelerin varlığı ile bu dengelere nasıl varıldığı tamamen ayrı sorunlar. Bu postada bahsettiğimiz dengeye varış dinamiğini çözmek için Lambert'in üç terimli denklemini çözmek gerekiyor. Bu denklem $z^{n} + q = z$ formundadır ve $n$ rasyonel olmak zorunda değildir.

Bütün birinci dereceden denge tepkimeleri hızla dengeye gelir

Kapalı ortamda ve sabit sıcaklıkta gerçekleşen en basit denge tepkimesini ele alalım. \begin{equation*} {\rm A} \ \rightleftharpoons \ {\rm B} \end{equation*} Gösterilen tepkimede ileri hız sabiti $\alpha > 0$, geri hız sabiti de $\beta > 0$ olsun. (Her iki hız sabitinin birimi de $1/{\rm s}$. Her nedense kimyasal kinetikte bu birim için Hz denmiyor...) A ve B maddelerinin zamana göre derişimleri ise $a(t)$ ve $b(t)$ ile verilsin. Sürekli karıştırmanın reaktör homojenliğini temin ettiğini varsayalım. Sistemin tanımını bitirmek için kütle aksiyon kanununu kullanarak hareket denklemlerini yazacağız. \begin{eqnarray} \nonumber \dot{a}(t) &=& -\alpha a(t) + \beta b(t) \\ \nonumber \dot{b}(t) &=& \alpha a(t) - \beta b(t) \end{eqnarray} Burada $\dot{x}(t)$ ile $x(t)$ fonksiyonunun zamana göre türevi gösterilmiştir. Başlangıç şartlarının en genel haliyle $a_{\rm o} := a(0) \geq 0$ ve $b_{\rm o} := b(0) \geq 0$ olduğunu belirtmekle yetineceğiz.

Dikkatli okur bu aşamada $\dot{a}(t) + \dot{b}(t) = \tfrac{d}{dt} (a(t)+b(t)) = 0$ olduğunu farketmiş olmalı. Bir niceliğin zamana göre türevinin sıfır olması, o niceliğin korunduğu veya sabit kaldığı manasına gelir. O zaman analizin temel teoremi uyarınca \begin{equation*} a(t) + b(t) = a(0) + b(0) = a_{\rm o} + b_{\rm o} =: m \geq 0 \end{equation*} yazabiliriz. Bu denklemdeki $m$ niceliğine mesela toplam kütle, ilgili denkleme de kütle dengesi diyebiliriz. Bir dinamik sistem problemini çalışırken eğer bir korunum kanunu bulunursa, bu eşitlik eldeki diferansiyel denklemlerden birisini elemekte kullanılabilir. Örneğin $b(t) = m - a(t)$ kullanılarak $a(t)$ için verilen adi diferansiyel denklem \begin{equation*} \dot{a}(t) = \beta m - ( \alpha + \beta ) a(t) \end{equation*} haline getirilebilir. $\rho := a/m$, $\kappa := \alpha / \beta$ ve $\tau := \beta t$ ile sırasıyla indirgenmiş derişim, denge sabiti ve birimsiz zaman niceliklerini tanımlayalım. Bu tanımlarla ve kütle korunumu ve pozitifliği kullanarak $\rho \in [0,1]$ olması gerektiği rahatça görülür. Bu tanımlarla $\rho$ için çözmemiz gereken hareket denklemi aşağıdaki gibi olur. \begin{equation*} \rho ^{\prime} (\tau) = 1 - (\kappa + 1) \rho (\tau) \end{equation*} Burada $\prime$ ile gösterilen türev zincir kuralı kullanılarak $\tau$ değişkenine göre alınmıştır. Problemin başlangıç şartı ise $\rho(0) = a(0)/m = a_{\rm o}/m$ ile verilecektir. İndirgenmiş niceliklerde problemin parametre sayısının $\kappa$ ve $\rho_{\rm o}$ olmak üzere ikiye düştüğünü gözleyiniz.

$\rho$ için yazdığımız adi diferansiyel denklemi çözmek için normalde Leibniz'in icat ettiği integral alma faktörü ile problemi bir tam diferansiyel haline getirmemiz gerekiyor. Ancak elimizdeki problem Leibniz tekniğinin bütün ayrıntılarını vermeden de çözülebilir. İntegral alma faktörü özü itibariyle bir diferansiyel denklemde çarpımın türevine ait $(fg)^{\prime} = f^{\prime}g+fg^{\prime}$ ifadesini ya bulmak ya da üretmekten ibarettir. Şimdi \begin{eqnarray} 1 &=& \rho^{\prime} + (\kappa+1) \rho(\tau) \\ \nonumber &=& \frac{1}{\exp((\kappa+1)\tau)} \left( \rho^{\prime} \exp((\kappa+1)\tau) + \rho(\tau) (\kappa+1)\exp((\kappa+1)\tau) \right) \\ \nonumber &=& \frac{1}{\exp((\kappa+1)\tau)} \frac{d}{d\tau} \left( \exp((\kappa+1)\tau) \rho(\tau) \right) \end{eqnarray} olduğundan $\rho$ fonksiyonu bir tam diferansiyel içine alınmış olur. Artık basitçe her iki tarafın integralini alarak matematiksel çözümü bitireceğiz. \begin{equation*} \int\limits_{0}^{\tau} \exp((\kappa+1)\sigma) d \sigma = \int\limits_{0}^{\tau} \frac{d}{d\sigma} \left( \exp((\kappa+1)\sigma) \rho(\sigma) \right) d \sigma \end{equation*} Denklemin sol tarafı için üstel fonksiyonun integralini, sağ tarafı için de analizin temel teoremini uygulayacağız. \begin{equation*} \frac{1}{\kappa + 1} \left( \exp((\kappa+1)\tau) - 1 \right) = \exp((\kappa+1)\tau) \rho(\tau) - \rho_{\rm o} \end{equation*} $\rho(\tau)$ fonksiyonunu yalnız bırakacak şekilde bu denklemi yeniden düzenleyerek matematiksel manipülasyonu noktalayacağız. \begin{equation*} \rho(\tau) = \frac{1}{\kappa + 1} + \left( \rho_{\rm o} - \frac{1}{\kappa+1} \right) \exp(-(\kappa+1)\tau) \end{equation*}

Kimyasal kinetikte bir problemin denge noktası nasıl bulunur? Cevap: üç yolla.

1. Hareket denklemlerini sıfırlayan derişimler kütle denkliği şartına tabi olacak şekilde çözülür. Diğer bir deyişle $\dot{a} = \dot{b} = -\alpha a + \beta b = 0$ ile $a + b = m$ denklemlerinin ortak çözümü bulunur.
2. Elimizde analitik çözümün bulunması halinde sistemin $\tau \to \infty$ limitinde dengeye geldiği varsayılarak, ki bu durumda gerçekten de öyledir, denge derişimi bulunur. \begin{equation*} \rho_{\rm d} := \lim_{\tau \to \infty} \rho (\tau) = \frac{1}{\kappa + 1} \end{equation*}
3. Son olarak problem sanki bir termodinamik problemiymiş gibi muamele edilir ve her ikisi de birinci dereceden olan \begin{equation*} \frac{\alpha}{\beta} = \kappa = \frac{b_{\rm d}}{a_{\rm d}} \ \ \ {\rm ve} \ \ \ a_{\rm d} + b_{\rm d} = m \end{equation*} iki bilinmeyenli iki denklemin ortak çözümü bulunur.

Biz bu postada en pahalı olan ikinci yolu tercih ettik zira postanın başlığında yer alan iddiayı ancak böyle kanıtlayabilirdik.

Çalıştığımız dinamik sistemin dengeye yaklaşma hızını nitel olarak ölçmek için aşağıdaki manipülasyonu takip edin. \begin{equation*} |\rho(\tau) - \rho_{\rm d}| = |\rho_{\rm o} - \rho_{\rm d}| \exp(-(\kappa + 1) \tau) \leq \exp(-(\kappa + 1) \tau) \end{equation*} $\rho$ niceliğinin tanımı gereği $|\rho_{\rm o} - \rho_{\rm d}| \leq 1$ olduğunu gözleyiniz. Bu eşitsizlik bize sistemin denge noktasına üstel hızda yakınsadığını söylemektedir ki üstel hız mevcut analitik fonksiyonlar ile elde edebileceğimiz en hızlı davranışlardan birisini temin eder. $|\rho(\tau) - \rho_{\rm d}| \leq 1 / \log (\kappa \tau)$ gibi bir davranış bulsaydık, o zaman reaksiyonun %95 oranında dengeye gelmesi için $1 / \log (\kappa \tau) = 0,05$ ya da $\kappa \tau = \exp (1/0,05) = 485.165.195,4 \approx 5 \times 10^{8}$ olması gerekirdi.

Öte yandan dengeye gelme zamanı ise birimsiz niceliklerde $\tau _{\rm d} := (\kappa + 1)^{-1}$ ile birimli niceliklerde ise $t_{\rm d}:=(\alpha + \beta)^{-1}$ ile verilir. Örneğin başlangıçtan $3 \tau_{\rm d}$ süre sonra sistem %95 itibariyle dengeye gelmiştir. ($\exp(-3) = 0,049787$) Her iki hız sabitinin toplamının sistemin dengeye gelme zamanını belirlediğini gözleyiniz.

Döngüsel bir kimyasal reaksiyon şebekesinde, konsantrasyonlar da şebeke gibi her zaman periyodik olur mu?

Hayır, olmaz.

Bu postada başlıktaki soruya verdiğimiz olumsuz yanıtı bir karşı örnekle gerekçelendireceğiz. Biyokimyada, bazı moleküllerin farklı izomerlerinin birbirlerine dönüşme kinetiği, oyun şebekesi (play network) dediğimiz yandaki döngüsel sisteme benzer. Şebekede yer alan üç reaksiyonun hız sabitleri sırasıyla $\alpha,\beta,\gamma$ ile ilgili reaksiyon okunun üzerinde gösteriliyor. Kimyasal kinetiğin en temel varsayımından, kütle aksiyon kanunundan başlayarak bu şebekede gerçekleşen reaksiyonlarda yer alan maddelerin konsantrasyonları için üç adet hareket denklemi yazacağız. \begin{eqnarray} \nonumber &&\dot{a}(t) = -\alpha a(t) + \gamma c(t) \\ \nonumber &&\dot{b}(t) = \alpha a(t) - \beta b(t) \\ \nonumber &&\dot{c}(t) = \beta b(t) - \gamma c(t) \end{eqnarray} Burada mesela $a(t)$ ile A maddesinin zamana bağlı derişimini, $\dot{a}(t)$ ile bu derişimin zamana göre türevini temsil ediyoruz. Hız sabitlerinin pozitif, başlangıç konsantrasyonlarının pozitif veya sıfır olması haricinde bu problem için başka da talep ettiğimiz bir şart yok. Şimdi hız denklemlerini taraf tarafa toplarsak $\dot{a} + \dot{b} + \dot{c} = 0$ elde ediyoruz. Dinamikte bir şeyin zamana göre türevi her zaman sıfırsa, bu, o niceliğin korunduğu anlamına gelir. Çalıştığımız problemde söz konusu olan basit bir toplam kütlenin korunumundan başka bir şey de değildir: \begin{equation*} a(t) + b(t) + c(t) = a(0) + b(0) + c(0) =: m \end{equation*} Dikkat edilirse kütle korunumunu kullanarak $c(t)$ maddesinin derişimi $a(t)$ ve $b(t)$ cinsinden yazılabilmektedir: $c(t) = m - a(t) - b(t)$. Bu ilişkiyi ilk iki hareket denklemine koyarak, çalışılması gerekli diferansiyel denklem sayısını üçten ikiye düşürebiliriz. Böylesi bir indirgeme dinamik sistemlerde çok tipiktir. Her bağımsız korunum kanunuyla, diferansiyel denklem sisteminden bir denklem elenebilir. Hatta hareket denklemi kadar korunum kanunu varsa, o zaman diferansiyel denklemleri çözmeye bile gerek kalmaz!

Adet olduğu üzere, bu sistem için denge durumunu soruşturmakla işe başlayacağız. Dinamik bir sistemi oluşturan durum değişkenlerinin hepsinin zamana göre türevinin sıfır olduğu noktalara denge noktaları denir. Bu tanımı kullandığımızda $\dot{a} = \dot{b} = \dot{c} = 0$ denklemlerinin ortak çözümü, bize denge konsantrasyonlarını verecektir. Basit bir alıştırma ile bu denge konsantrasyonlarını hesaplayabiliriz. \begin{eqnarray} \nonumber a_{\rm d} &=& \frac{\beta \gamma}{\alpha \beta + \alpha \gamma + \beta \gamma} m \\ \nonumber b_{\rm d} &=& \frac{\alpha \gamma}{\alpha \beta + \alpha \gamma + \beta \gamma} m \\ \nonumber c_{\rm d} &=& \frac{\alpha \beta}{\alpha \beta + \alpha \gamma + \beta \gamma} m \end{eqnarray} Şebekedeki bütün reaksiyonlar tek yönlü, tersinmez (irreversible) gibi görünüyor ama sistemde yer alan hiçbir maddenin denge konsantrasyonu sıfır değil! Bu genellikle kimyada görmeye alışık olmadığımız ve şebekenin döngüselliğinin temin ettiği bir durum. Ayrıca bu örnekten ilham alarak tersinirlik kavramını biraz daha genişletiyor ve aşağıdaki tanımı yapıyoruz.

Tanım: (Tersinirlik ve zayıf tersinirlik) Bir kimyasal reaksiyon şebekesinde yer alan bütün tepkimeler tersinir (reversible) ise, o zaman o şebekeye tersinir şebeke denir. Bir kimyasal reaksiyon şebekesinde yer alan reaksiyon oklarının her iki tarafında yer alan reaktant ve ürünlere topluca kompleks denir. Bir kimyasal reaksiyon şebekesinde yer alan herhangi bir ${\mathcal C}_{1} \to {\mathcal C}_{2}$ reaksiyonu için, ${\mathcal C}_{2}$ kompleksi ile başlayıp ${\mathcal C}_{1}$ kompleksi ile biten bir yol (yani reaksiyon zinciri) bulunabiliyorsa, o zaman o şebekeye zayıf tersinir şebeke (weakly reversible) denir.

Bütün tersinir şebekelerin aynı zamanda zayıf tersinir olduğu çok barizdir. Burada çalıştığımız oyun şebekesi tersinir değil. Çünkü, örneğin ${\rm A} \to {\rm B}$ reaksiyonu var ama ${\rm B} \to {\rm A}$ yok. Öte yandan ${\rm B} \to {\rm C} \to {\rm A}$ kanalıyla B maddesinden A maddesine ulaşmak mümkün olduğundan çalıştığımız oyun şebekesi zayıf tersinirdir.

Şebekenin dinamiğini çözmeden önce birimsiz niceliklere geçeceğiz. $t =: \tau/\alpha$, $A := a/m$, $B := b/m$, $C :=c/m$, $g := \gamma/\alpha$ ve $h := \beta/\alpha$ tanımlayalım. $A+B+C=1$ olduğunu gözleyiniz. Dahası kütlenin korunumu gereğince $C = 1 - A - B$ yazabiliriz. Bu, $\dot{c}$ için yazılan diferansiyel denklemi fuzuli (redundant) kılar. İndirgenmiş birimlerde çalışmamız gereken dinamik sistem aşağıdaki gibidir. \begin{eqnarray} \nonumber A^{\prime}(\tau) &=& -(1+g)A(\tau) - gB(\tau) + g \ \ \ (1) \\ \nonumber B^{\prime}(\tau) &=& A(\tau) - h B(\tau) \ \ \ (2) \end{eqnarray} $\tau$ değişkenine göre türevi $\prime$ ile temsil ettik. Şimdi (1) nolu denklemin bir daha ($\tau$ değişkenine göre) türevini alır ve (2) nolu denklemi de kullanırsak, aşağıdaki ara sonucu elde ediyoruz. \begin{equation*} A^{\prime \prime}(\tau) = -(1+g)A^{\prime}(\tau) - gA(\tau) + ghB(\tau) \ \ \ (3) \end{equation*} Amacımız sadece $A$ değişkenine bağlı bir adi diferansiyel denklem elde etmek. Bu yüzden (1) nolu denklemden $B(\tau)$ ifadesini çekip, (3) nolu denklemde kullanınca çalışmamız gereken denklem aşağıdaki forma girmektedir. \begin{equation*} A^{\prime\prime}(\tau) + (1+g+h)A^{\prime}(\tau) + (g+h+gh)A(\tau) = gh \ \ \ (4) \end{equation*} Bu denklemin sadece $A$ değişkenine bağlı olduğunu gözleyiniz. Sistemin toplamda $g$ ve $h$ olmak üzere iki adet parametresi vardır.

(4) nolu denklem ikinci dereceden, sabit katsayılı, homojen olmayan, lineer bir adi diferansiyel denklemdir ve adi diferansiyel denklemlerin teorisinde çözüm yöntemi vardır. İlkin tekabül eden karakteristik denklemi çözeceğiz. \begin{equation*} \lambda^{2} + (1+g+h)\lambda + (g+h+gh) = 0 \ \ \ (5) \end{equation*} İkinci dereceden denklemin diskriminantı aşağıdaki gibidir. \begin{equation*} \Delta := (1+g+h)^{2}-4(g+h+gh) = (g-h)^{2} - 2(g+h) + 1 \end{equation*} Bu diskriminant hem pozitif hem de negatif olabilir. Örneğin $g=h=1$ için $\Delta = -1$ ama $g=5$, $h=1$ için $\Delta = 5$ olmaktadır. Karakteristik değerler aşağıdaki gibi verilir. \begin{equation*} \lambda_{1,2} :=\frac{-(1+g+h)\pm\sqrt{\Delta}}{2} \ \ \ (6) \end{equation*} Viete-Girard formüllerinden $\lambda_{1}+\lambda_{2} < 0$ ve $\lambda_{1}\lambda_{2} > 0$ olduğundan, her durumda köklerin gerçel kısımları negatiftir.

$A_{\rm p} := gh/(g+h+gh)$ ifadesinin (4) nolu denklemi sağladığını gözleyiniz. (Basitçe yerine koymanız yeterlidir.) İndirgenmiş birimlerde A malzemesinin denge konsantrasyonunu da veren bu ifadeye kısmi çözüm diyeceğiz. Kısmi çözüm aynı zamanda $A_{\rm p} = a_{\rm d}/m$ denklemini de sağladığından, aslında indirgenmiş birimlerde A maddesinin denge konsantrasyonudur. Artık aradığımız çözümü nihayet verebiliriz. \begin{equation*} A(\tau) = A_{\rm p} + c_{1}e^{\lambda_{1}\tau} + c_{2}e^{\lambda_{2}\tau} \ \ \ (6a) \end{equation*} $c_{1,2}$ katsayıları başlangıç şartlarından temin edilmelidirler. $\tau=0$ koyduğumuzda \begin{equation*} c_{1}+c_{2} = A(0)-A_{\rm p} \ \ \ (7) \end{equation*} denklemini elde ediyoruz. İkinci bir denkleme daha ihtiyacımız var. Bu amaçla (1) nolu denklemde de $\tau=0$ koyacağız. \begin{equation*} \lambda_{1}c_{1} + \lambda_{2}c_{2} = -(1+g)A(0) - gB(0) + g \ \ \ (8) \end{equation*} (7) ve (8) nolu denklemler beraber çözüldüklerinde aşağıdaki sonuçları elde ediyoruz. \begin{eqnarray}\nonumber c_{1} &=& \frac{-\lambda_{2}(A(0)-A_{\rm p})-(1+g)A(0)-gB(0)+g}{\sqrt{\Delta}} \ \ \ (9) \\ \nonumber c_{2} &=& \frac{\lambda_{1}(A(0)-A_{\rm p})+(1+g)A(0)+gB(0)-g}{\sqrt{\Delta}} \ \ \ (10) \end{eqnarray}

Durum I. $\Delta > 0$ için (6a) nolu denklemi kullanarak $\lim_{\tau \to \infty} = A_{\rm p}$ olduğunu gösterebiliriz. Dahası yine aynı denklemin türevini aldıktan sonra $\lim_{\tau \to \infty}A^{\prime}(\tau) = 0$ olduğunu da gösterebiliriz. Bu bize (1) nolu denklem kanalıyla $\lim_{\tau \to \infty}B(\tau) = g/(g+h+gh)$ sonucunu verir. İndirgenmiş birimlerde bu, B maddesinin denge konsantrasyonundan başka bir şey değildir. Kütle korunumu ile C maddesinin de $\tau \to \infty$ limitinde, dengeye geldiği gösterilir. Sistemin dengeye varma süresi, indirgenmiş birimlerde $|\lambda_{2}|^{-1}$ ile kestirilebilir. (Neden $|\lambda_{1}|^{-1}$ değil?) Birimli niceliklerde bu değer $T \sim (2/\alpha) / (1+g+h-\sqrt{\Delta})$ kadardır.

Durum II. $\Delta = 0$ için $\lambda_{1} = \lambda_{2} =: \Lambda = -(1+g+h)/2$ olacaktır. Bu durumda çözümü \begin{equation*} A(\tau) = A_{\rm p} + d_{1}e^{\Lambda \tau} + d_{2}\tau e^{\Lambda \tau} \ \ \ (11) \end{equation*} şeklinde ifade etmeliyiz. ((11) nolu denklemi (4) nolu denkleme koyarak bu çözümün doğruluğunu gösteriniz.) Burada $d_{1}$ ve $d_{2}$ başlangıç şartlarından tayin edilmelidirler. Basit bir alıştırmayla bu katsayıları çözebiliriz. \begin{eqnarray}\nonumber d_{1} &=& A(0) - A_{\rm p} \\ \nonumber d_{2} &=& -\Lambda (A(0)-A_{\rm p}) - (1+g)A(0) - gB(0) + g \end{eqnarray} $\Lambda < 0$ olduğundan bu sistem de üstel hızda dengeye varır ve dengeye varma zamanı birimli niceliklerde $T \sim (2/\alpha)/(1+g+h)$ ile kestirilebilir.

Durum III. $\Delta < 0$ için özdeğerlerin gerçel olmadığını gözleyiniz. Bu durumu çalışmak için öncelikle $\delta := \sqrt{|\Delta|}$ tanımıyla işe başlıyoruz. O zaman özdeğerler $\lambda_{1,2} = \Lambda \pm \tfrac{i\delta}{2}$ ile verilecektir. Burada $\Lambda$ bir önceki paragrafta tanımlandığı gibidir. Bu tanımlarla ve trigonometrik fonksiyonların ($2i\sin(\theta) = e^{i\theta}-e^{-i\theta}$ ve $2\cos(\theta) = e^{i\theta}+e^{-i\theta}$ gibi) bazı özelliklerini kullanarak A maddesinin konsantrasyonu aşağıdaki gibi verilir. \begin{equation*} A(\tau) = A_{\rm p} + e^{\Lambda \tau}(A(0)-A_{\rm p}) \cos\left( \frac{\delta \tau}{2} \right) - e^{\Lambda \tau}\frac{2}{\delta}\left((1+g)A(0)+gB(0)-g+\Lambda(A(0)-A_{\rm p})\right) \sin\left( \frac{\delta \tau}{2} \right) \end{equation*} Yukarıdaki denklem bazı tanımlarla sadeleştirilebilir. \begin{eqnarray} \nonumber &&D:=\sqrt{(A(0)-A_{\rm p})^{2}+\frac{4}{\delta^{2}}\left((1+g)A(0)+gB(0)-g+\Lambda(A(0)-A_{\rm p})\right)^{2}} \\ \nonumber &&\cos \varphi :=\frac{A(0)-A_{\rm p}}{D} \\ \nonumber &&\sin \varphi := \frac{2}{\delta}\frac{(1+g)A(0)+gB(0)-g+\Lambda(A(0)-A_{\rm p})}{D} \end{eqnarray} Nihayet $\cos(x+y) = \cos x \cos y - \sin x \sin y$ özdeşliğini kullandığımızda A malzemesinin konsantrasyonunu daha sade bir formda sunabiliyoruz. \begin{equation*} A(\tau) = A_{\rm p} + De^{\Lambda \tau} \cos \left( \frac{\delta \tau}{2}+\varphi \right) \ \ \ (12) \end{equation*}

İşaret. $\Lambda < 0$ olduğundan, burada da $\lim_{\tau \to \infty}A(\tau)=A_{\rm p}$ olur, yani tepkime dengeye üstel hızda gelir. Ama (12) nolu denklem şimdiye kadar hiç karşılaşmadığımız bir davranışa, salınımlara (oscillations) sahiptir. Biyokimyada ve kimya mühendisliğinde salınım yapan tepkimeler önemli bir yer tutar. Ne yazık ki üstel terimin hızla sıfıra gitmesinden ötürü, oyun şebekesinde salınımları uzun süre gözlemek mümkün değildir. Bu meyanda Deficiency-0 teoreminin pek çok kimyasal sistemde salınımları gözlemenin imkansız olduğunu söylediğini kaydedelim.

İşaret. Birimsiz niceliklerde tepkime zamanı $|\Lambda|^{-1}=2/(1+g+h)$ ile kestirilebilir. Öte yandan salınımların periyotu tam olarak \begin{equation*} \frac{4\pi}{\delta} = \frac{4\pi}{\sqrt{ 2(g+h) - (g-h)^{2} - 1}} \end{equation*} kadardır. En az bir tam salınım gözlemek için $|\Lambda|^{-1} \ge 4\pi \delta^{-1}$ şartını kullanmamız gerekiyor. Bu şart aşağıdaki eşitsizliği gerektirmektedir. \begin{equation*} 2(1-4\pi^{2})(g+h+gh) - (1+4\pi^{2})(g^{2}+h^{2}) \ge 4\pi^{2} - 1 \end{equation*} Ne var ki bu eşitsizlik absurddur. Zira negatif bir niceliğin pozitif bir nicelikten büyük olduğunu söyler. Diğer bir ifadeyle oyun şebekesinde bir tam salınım dahi gözlenmeden şebeke dengeye gelir.

Üçüncü dereceden bir analog gürültü filtresi tasarımı

Arka plan ve motivasyon

Hepimiz spektroskopi öğrendik ama ne dört başı mamur bir gürültü devresi tasarımı yaptık ne de paraboloid aynaların mükemmel kolimasyon ve odaklama özelliklerini ispatladık. Daha önce yerölçüsünde paraboloid aynalarla ilgili çalışmalar yaptığımızdan bu postada gürültü devresi tasarımına bakacağız. Neden? Zira spektroskopik cihazlar, diğer tüm elektronik cihazlar gibi gürültü üretir. Gürültünün muhtelif fiziksel nedenleri var, ancak bu postada gürültünün muayenesinden ziyade bertarafına odaklanacağız. Kabaca söylemek gerekirse, bir elektrik devresindeki yinelenebilir çıktıya -ne kadar zayıf olursa olsun- sinyal ve yinelenemeyen kısma -ne kadar güçlü olursa olsun- gürültü denir. Daha teknik bir söylemle sinyal deterministik, gürültü stokastiktir. Bir çalışmanın odağı, gürültünün kendisinin istatistiksel bir analizi değilse, o zaman gürültünün varlığı genellikle istenmez. Gürültüyü temizlemenin bir yolu, onun tanımından faydalanır. Aynı deneyi defalarca yineler ve daha sonra da çıktılarını eklersiniz. Tesadüfi gürültü genlikleri genellikle toplamda çok küçük sayılara baliğ olurken, deterministik sinyal kesin bir şekilde bu süreçte güçlenir. Her ne kadar kelimenin gerçek anlamıyla mükemmel olsa da, bu yöntem çok pahalı ve zaman alabilecek bir deneyi defalarca yineleme imkanımızın olduğunu varsaymaktadır ki bu genellikle hiçbir zaman mümkün olmaz. Bu noktada elektronik alternatif bir çözüm sunar. Bu postada analog ve muhtemelen cihaza gömülü lineer devre unsurlarından oluşmuş bir gürültü devresi tasarımını, işin kararlılık boyutunu da ele alarak irdeleyeceğiz.

Gerek optik gerekse manyetik spektroskopide alışık olduğumuz laboratuvar tecrübesi, gürültü ve sinyalin hemen hemen her zaman, frekans uzayında ayrıştıklarını, gürültünün yüksek frekanslı sinyalin ise nisbeten daha düşük frekanslı olduğunu söylemektedir. Buna geçici varsayım (working hypothesis) diyeceğiz ve pek çok sisteme uygulaması olduğundan böylesi bir probleme maliyeti düşük bir çözüm önereceğiz. Geçici varsayım uyarınca, günlük konuşma dilindeki anlamıyla bir sinyalin Fourier analizini yapan ve onun düşük ve yüksek frekanslı bileşenlerini ayrıştırıp farklı bir biçimde muamele eden bir devre olsa, bu kuşkusuz spektroskopik cihazlarda çok büyük fayda temin edecektir. Okurun anaokulu düzeyinde RC devresinin analizine ve frekans davranışına aşina olduğunu varsayıyoruz. Bu postada inceleyeceğimiz üç ızgaralı analog devre aşağıdaki panoda gösteriliyor. Devrenin bir varyantını Omar Wing'in ders kitabında gördüm. Bu postada amacım söz konusu devreyi en genel haliyle inceleyip, devre parametrelerini düşük frekanslı çıktıyı en ideal bir şekilde iletirken, yüksek frekanslı olanları bastıracak şekilde optimize etmektir. Bu işi yaparken bir kararlılık analizini de sunmayı hedefliyorum.

Devrenin dinamiğine ilişkin cebirsel ve diferansiyel denklemlerin türetilmesi

Bu devrede değerleri üzerinde gösterilen iki direnç iki de sığaç mevcuttur. Söylemeye lüzum yok, $\alpha$ ve $\beta$ birimsiz ölçek çarpanları olarak bu devrede rol alıyorlar. Kirchhoff'un akım ve voltaj kanunlarını bu devreye uyguladığımızda aşağıdaki karışık cebirsel ve diferansiyel denklem sistemini elde ediyoruz. Direnç, sığaç, bobin gibi pasif elektronik devre unsurlarını bizim için özel kılan şey, bunlarla kurulan devrelerin dinamiğinin lineer cebirsel ve lineer adi diferansiyel denklemlere tabi olmasıdır. Lineer denklemlere ilişkin güçlü ve zengin bir yöntem mühimmatını bu problemlere uygulayabiliriz. Önce denklemleri yazalım. \begin{eqnarray} &{\rm A1:}&\ i(t) = i_{1}(t) + i_{2}(t) \nonumber \\ &{\rm A2:}&\ i_{1}(t) = i_{3}(t) + i_{4}(t) \nonumber \\ &{\rm V1:}&\ Ri(t) + \frac{1}{C}q_{2}(t) -V_{I}(t)=0 \nonumber \\ &{\rm V2:}&\ L\frac{{\rm d}i_{1}(t)}{{\rm d}t} + \frac{1}{\alpha C} q_{4}(t) - \frac{1}{C}q_{2}(t) =0 \nonumber \\ &{\rm V3:}&\ \beta R i_{3}(t) - \frac{1}{\alpha C} q_{4}(t) =0 \end{eqnarray} Burada voltaj denklemlerinde yer alan $q_{2}(t)$ ve $q_{4}(t)$, sırasıyla, $C$ ve $\alpha C$ sığaçlarında depolanan anlık yükleri temsil ediyorlar. Öte yandan elektronik laboratuvarlarında adet yük yerine akımları ölçmek ve bir durum değişkeni olarak yükleri topyekun elemektir. Buna göre, ihtiyaç halinde, $\dot{q}_{2}(t) = i_{2}(t)$ vb. eşitlikler sayesinde denklemlerimizden yükleri eleyebiliriz. Her zaman olduğu gibi bir değişkenin tepesine konulan nokta onun zamana göre türevini simgelemektedir. Burada göstermeyeceğiz ama voltaj denklemlerinde yer alan yükleri elemek amacıyla söz konusu denklemlerin her iki tarafının da zamana göre türevini almamız yeterli.

Temel devre teorisinden biliyoruz ki $RC$ çarpımının birimi zaman. Bu gözlem bize $t =: \tau RC$ kanalıyla yeni bir birimsiz zaman değişkeni tanımlama imkanı sunuyor. Değişkenin tepesine koyduğumuz bir üssü (prime) işareti ile $\tau$'ya göre türev almayı temsil edersek, o zaman adi türevlerdeki zincir kuralı uyarınca zamana bağlı herhangi bir fonksiyon için $\dot{f}(t) = (RC)^{-1} f^{\prime}(\tau)$ eşitliğini not edelim. Bundan böyle Kirchoff denklemlerini birimsiz zaman değişkenleriyle ifade edilmiş bir biçimde kullanacağız. Az sonra nedenleri daha belirgin bir halde tebarüz edeceği gibi, $i_{3}$ değişkenini problemin merkezi değişkeni addedip, geri kalan diğer durum değişkenlerini $i_{3}$ ve onun zamana göre türevleri cinsinden ifade edeceğiz. Şimdi (V3) denkleminin önce her iki tarafının zamana göre türevini alıp ($0=\beta R \dot{i}_{3} - i_{4}/\alpha C$) ardından birimsiz zaman niceliğine geçtiğimizde aşağıdaki sonucu elde ediyoruz. \begin{equation*} {\rm V3a}: \ i_{4}(\tau) = \alpha \beta i_{3}^{\prime}(\tau) \end{equation*} Ayrıca (V3) denklemini yeniden düzenlediğimizde basitçe aşağıdaki denklemi de elde ederiz. \begin{equation*} {\rm V3b}: \ q_{4}(\tau) = \alpha \beta RC i_{3}(\tau) \end{equation*} (V3a) denkleminden (A2) denkleminde faydalandığımızda ise $i_{1}$ merkezi değişken cinsinden yazılmış olur. \begin{equation*} {\rm A2a}: \ i_{1}(\tau) = i_{3}(\tau) + \alpha \beta i_{3}^{\prime}(\tau). \end{equation*} Şimdi (A2a) ve (V3b) denklemlerini (V2) denkleminde kullanırsak, o zaman \begin{equation*} {\rm V2a}: \ q_{2}(\tau) = \frac{L}{R} (i_{3}^{\prime}(\tau) + \alpha \beta i_{3}^{\prime \prime}(\tau)) + \beta RC i_{3}(\tau) \end{equation*} elde ediliyor. Denklemin her iki tarafının $t$'ye göre türevini aldığımızda ise, $i_{2}$ değişkenini de $i_{3}$ cinsinden ifade etmiş olacağız. \begin{equation*} {\rm V2b}: \ i_{2}(\tau) = \alpha \beta \kappa i_{3}^{\prime \prime \prime} (\tau) + \kappa i_{3}^{\prime \prime}(\tau) + \beta i_{3}^{\prime}(\tau) . \end{equation*} Burada $\kappa := L/R^{2}C$ problemi tanımlayan bütün parametrelerin öbeklendiği birimsiz bir parametredir ve devrenin hem nitel hem de nicel davranışı $\kappa$ vasıtasıyla belirlenebilir. (A2a) ve (V2b) denklemlerinden faydalanarak (C1) etiketli denkleme geri döndüğümüzde aşağıdaki sonucu elde ediyoruz. \begin{equation*} {\rm A1a}:\ i(\tau) = \alpha \beta \kappa i_{3}^{\prime \prime \prime}(\tau) + \kappa i_{3}^{\prime \prime}(\tau) + (1+\alpha)\beta i_{3}^{\prime}(\tau) + i_{3}(\tau) \end{equation*}

Bu manipulasyonlar neticesinde tüm durum değişkenlerini $i_{3}$ ve türevleri cinsinden ifade etmiş olduk. Dolayısıyla $i_{3}$ değişkenini çözmek aslında problemi çözmekle eş anlamlıdır. Eğer yalnızca $i_{3}$ değişkenini içeren bir denklem elde etmek istiyorsak, o zaman şimdiye kadar hiç kullanmadığımız bir Kirchhoff denklemi kullanmalıyız ki o denklem de (V1) etiketli denklemdir. (A1a) ve (V2a) söz konusu denklemde kullanıldığında aşağıdaki lineer, üçüncü dereceden, homojen olmayan, sabit katsayılı adi diferansiyel denklem elde edilir. \begin{equation*} \boxed{ \alpha \beta \kappa i_{3}^{\prime \prime \prime}(\tau) + \kappa (1+ \alpha \beta) i_{3}^{\prime \prime} (\tau) +(\alpha \beta + \beta + \kappa) i_{3}^{\prime}(\tau) + (1+\beta) i_{3}(\tau) = \frac{1}{R} V_{I}(\tau) } \end{equation*}

Ancak başlangıç şartlarını vermek suretiyle bu diferansiyel denklemin tam çözümü talep edebiliriz. $i_{3}(\tau)$ için verilen adi diferansiyel denklem üçüncü derece olduğundan, üç tane başlangıç şartına ihtiyacımız vardır. Bunlar $i_{3}(0)$, (V3a) denklemi kanalıyla $i_{3}^{\prime}(0) = i_{4}(0)/\alpha \beta$ ve (V2a) denklemiyle $i_{3}^{\prime \prime}(0) = \frac{1}{\alpha \beta} \left( \frac{R}{L}q_{2}(0)-\frac{1}{\alpha \beta}i_{4}(0)-\frac{\beta}{\kappa}i_{3}(0)\right)$ olarak temin edilebilirler. Başlangıçta boş sığaçlarla ve yavaşça sıfır değerinden artan ya da azalan girdi gerilimiyle başlamak her zaman için iyi bir uygulamadır. Dolayısıyla genelliği biraz kaybederek $i_{3}(0)=i_{3}^{\prime}(0) = i_{3}^{\prime \prime}(0)=0$ olduğunu kabul edebiliriz.

Homojen çözüm ve bu çözümün kararlılığı

Homojen olmayan bir adi diferansiyel denklemin çözüm sürecinin ilk basamağı sanki denklem homojenmiş gibi davranıp tekabül eden homojen denklemin çözülmesidir, diğer bir deyişle kutudaki denklemin sağ tarafını sıfıra eşitleyip çözümüne bakacağız. Tekabül eden homojen denklem sabit katsayılı ve lineer olduğundan, $\lambda$ problemin özdeğerlerini temsil etmek üzere, ilgili denklemin $i_{3}(\tau) = \exp(\lambda \tau)$ formunda çözümlerini arayacağız. Böylesi bir deneme çözümü denklemin sıfatını diferansiyelden cebirsele, yani $\tilde{p}(\lambda)=0$ denkleminin kök bulma problemine dönüştürür. Burada $\tilde{p}$ aşağıdaki gibi tanımlanan kübik bir polinomdur. \begin{equation*} \tilde{p}(x) := \alpha \beta \kappa x^{3} + \kappa (1+ \alpha \beta) x^{2} + (\alpha \beta + \beta + \kappa) x + 1 + \beta \end{equation*} $\lambda_{i}$, $\tilde{p}$ polinomunun köklerini temsil etmek üzere, tekabül eden homojen çözüm aşağıdaki gibi olur. \begin{equation*} {\rm H1}: \ i_{3}^{\rm h}(\tau) = C_{1} \exp(\lambda_{1}\tau) + C_{2} \exp(\lambda_{2}\tau) + C_{3} \exp(\lambda_{3}\tau) \end{equation*} Burada $C_{i}$ yukarıda bahsi geçen başlangıç şartlarından tayin edilmesi gereken katsayılardır. Şu an için yozlaşma (degeneracy) yani köklerin çokkatlılığı ihtimalini dikkate almayacağız. Bir polinomdaki bütün katsayıları baş katsayıya böldüğümüzde elde edilen (monik) polinomla asıl polinomun kökleri aynıdır. Bu yüzden bundan sonra $\tilde{p} =: \alpha \beta \kappa p$ denklemiyle verilen $p$ polinomunu inceleyeceğiz. \begin{equation*} p(x) := x^{3} + \left( 1+\frac{1}{\alpha \beta} \right) x^{2} + \left( \frac{1}{\kappa}+\frac{1}{\alpha \kappa}+\frac{1}{\alpha \beta}\right) x +\frac{1}{\alpha \beta \kappa} + \frac{1}{\alpha \kappa} \end{equation*}

$\lambda_{i}$ değerlerinden en az birisinin pozitif gerçel kısmı olması halinde, ilgili katsayı eğer sıfır değilse -ki bu aşırı derecede nadir bir durumdur-, o zaman homojen çözüm asimptotik olarak kararsız olur. Diğer bir deyişle $\lim_{\tau \to \infty} i_{3}^{\rm h}(\tau) \to \pm \infty$ olacaktır. Böylesi devrelerin çok yüksek akım sonucu eriyecekleri aşikardır. Öte yandan eğer tüm $\lambda_{i}$ değerlerinin gerçel kısımlar negatif ise, o zaman $\lim_{\tau \to \infty} i_{3}^{\rm h} \to 0$ olur. Başka bir deyişle devre asimptotik olarak kararlıdır ve homojen çözümün katkısı ancak geçici bir süreliğine varlığını hissettirir. Yani $\tau > |1/\min \Re[\lambda_{i}]|$ ise, o zaman $i^{\rm h}_{3}(\tau)$ ihmal edilebilir. (Burada $\Re[z]$ ile karmaşık bir sayının gerçel kısmını temsil ediyoruz.)

Bu noktaya değin sunduğumuz sezgisel argümanlar bizi Routh-Hurwitz kararlılık problemine getiriyor: Katsayıları gerçel bir polinomun bütün köklerinin gerçel kısımlarının negatif olup olmadığını kat'i olarak tayin etmek için sadece sonlu sayıda aritmetik işlemden faydalanan bir algoritma icat ediniz. 19. yy'da ifade edilen ve de çözümü verilen bu problemin muhtelif çözümleri mevcut. Daha önce yerölçüsünde de şerhli tercümesini verdiğimiz Strelitz algoritması ile bu problemi çözeceğiz. Teoremin burada yer vereceğimiz lafzı Prasolov'un kitabından alınmadır.

Teorem: (Sh. Strelitz) $p(z):=z^{n}+a_{1}z^{n-1}+\cdots+a_{n}$ katsayıları gerçel, monik bir polinom olsun. $q(z):=z^{m}+b_{m-1}z^{m-1}+\cdots+b_{m}$ ise $m:=\tfrac{1}{2}n(n-1)$ olmak üzere, kökleri $p$ polinomunun köklerinin farklı çiftlerinin toplamı olan polinom olsun. $p$ polinomunun Routh-Hurwitz kararlı olması için gerek ve yeter şart, hem $p$ hem de $q$ polinomlarının katsayılarının pozitif olmasıdır.

Bu teoremin ispatını yukarıda verdiğimiz linkten takip edebilirsiniz. Her ne kadar burada kullanmayacak olsak da, Strelitz'in makalesinde $q$ polinomunun katsayılarını Newton toplamlarından faydalanarak, $p$ polinomunun katsayıları cinsinden ifade eden bir algoritma da mevcut. Bizim problemimizde hem $p$ hem de $q$ polinomlarının derecesi 3 gibi küçük olduğundan genel bir algoritmaya girmeyeceğiz.

Strelitz algoritmasındaki kilit nokta $b_{i}$ katsayılarının kurulumudur. $p(z):=z^{3}+a_{1}z^{2}+a_{2}z+a_{3}$ olsun ve kökleri de $\lambda_{1}$, $\lambda_{2}$ ve $\lambda_{3}$ ile temsil edilsin. O zaman tanım gereği $q$ polinomunun kökleri $\lambda_{1}+\lambda_{2}$, $\lambda_{1}+\lambda_{3}$ ve $\lambda_{2}+\lambda_{3}$ olmalıdır. Dolayısıyla \begin{equation*} b_{1}=-(\lambda_{1}+\lambda_{2}+\lambda_{1}+\lambda_{3}+\lambda_{2}+\lambda_{3})=-2(\lambda_{1}+\lambda_{2}+\lambda_{3})=2a_{1} \end{equation*} sonucunu hemen elde ederiz. $b_{2}$ katsayısı için aşağıdaki küçük cebirsel manipulasyonlara başvuracağız. \begin{eqnarray} \nonumber b_{2} &=& (\lambda_{1}+\lambda_{2})(\lambda_{1}+\lambda_{3}) + (\lambda_{1}+\lambda_{2})(\lambda_{2}+\lambda_{3}) + (\lambda_{1}+\lambda_{3})(\lambda_{2}+\lambda_{3}) \\ \nonumber &=& \lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2}+3(\lambda_{1}\lambda_{2}+\lambda_{1}\lambda_{3}+\lambda_{2}\lambda_{3}) \\ \nonumber &=& (\lambda_{1}+\lambda_{2}+\lambda_{3})^{2} + \lambda_{1}\lambda_{2}+\lambda_{1}\lambda_{3}+\lambda_{2}\lambda_{3} \\ \nonumber &=& a_{1}^{2}+a_{2} \end{eqnarray} $a_{1}$ ve $a_{2}$ pozitif ise, o zaman $b_{1}$ ve $b_{2}$ değerlerinin de pozitif olacağı barizdir. Gerçekten de ele aldığımız problemde $a_{1} := 1+(\alpha\beta)^{-1}$, $a_{2}:=\kappa^{-1}+(\alpha\kappa)^{-1}+(\alpha\beta)^{-1}$ ve $a_{3}:=(\alpha\beta\kappa)^{-1}+(\alpha\kappa)^{-1}$ hep pozitif niceliklerdir. Stabilite analizini $b_{3}$ için bir ifade vererek tamamlayacağız. \begin{eqnarray}\nonumber b_{3} &=& - (\lambda_{1}+\lambda_{2})(\lambda_{1}+\lambda_{3})(\lambda_{2}+\lambda_{3}) \\ \nonumber &=&-(\lambda_{1}+\lambda_{2}+\lambda_{3}-\lambda_{3})(\lambda_{1}+\lambda_{2}+\lambda_{3}-\lambda_{2})(\lambda_{1}+\lambda_{2}+\lambda_{3}-\lambda_{1}) \\ \nonumber &=&-(-a_{1}-\lambda_{3})(-a_{1}-\lambda_{2})(-a_{1}-\lambda_{1}) \\ \nonumber &=&a_{1}^{3}+a_{1}^{2}(\lambda_{1}+\lambda_{2}+\lambda_{3}) + a_{1}(\lambda_{1}\lambda_{2}+\lambda_{1}\lambda_{3}+\lambda_{2}\lambda_{3}) + \lambda_{1}\lambda_{2}\lambda_{3} \\ \nonumber &=&a_{1}a_{2}-a_{3} \end{eqnarray} Problemimizdeki yegane aşikar olmayan işaret sorgusu da budur. Değerleri yerine koyduğumuzda \begin{equation*} b_{3} = \kappa^{-1} + (\alpha \beta)^{-1} + (\alpha^{2}\beta\kappa)^{-1} + (\alpha^{2}\beta^{2})^{-1} > 0 \end{equation*} olduğunu da gösterebiliriz. Dolayısıyla Strelitz teoremi uyarınca devremizin kayıtsız şartsız Routh-Hurwitz kararlı olduğunu ve homojen çözümün varlığını ancak geçici bir süreliğine gösterdiğini ispatlamış olduk.

Kısmi ve toplam çözümler

Homojen olmayan adi diferansiyel denklemlerin çözümlerinde, homojen çözüm başlangıç şartlarını çözüme yansıtırken, kısmi çözüm ise tahrik gücünden gelen katkıyı gösterir. Girdi voltajının alternatif/sinusoidal olduğunu kabul edeceğiz. Söz konusu voltaj birimsiz zaman ve frekans değişkenlerinde aşağıdaki gibidir. \begin{equation*} V_{I}(\tau) = V_{\circ} \sin (\Omega \tau) \end{equation*} $V_{I}$ fonksiyonunun trigonometrik formundan ötürü, aşağıdaki şablona uyan kısmi çözümler arayacağız. \begin{equation*} {\rm P1}: \ i_{3}^{\rm p}(\tau) := A\sin(\Omega \tau) + B\cos(\Omega \tau) \end{equation*} Buradaki $A$ ve $B$ katsayıları $i_{3}^{\rm p}(\tau)$ fonksiyonunun kutu içindeki diferansiyel denklemi sağlaması gerektiği şartından temin edilecektir. (Dikkat: başlangıç şartlarından değil!) Müteakip işlemler usandırıcı olsa da, işlem hatası yapmamak kaydıyla, kolaydır. (P1) etiketli deneme çözümünü kutudaki denkleme koyup türevleri aldıktan sonra cosinus ve sinus terimlerini aşağıdaki gibi ayrıştırıyoruz. \begin{eqnarray}\nonumber &&\left\{-\alpha\beta\kappa \Omega^{3}A-\kappa(1+\alpha\beta)\Omega^{2}B+(\alpha\beta+\beta+\kappa)\Omega A+(1+\beta)B \right\} \cos(\Omega \tau) \\ \nonumber &+& \left\{\alpha\beta\kappa \Omega^{3}B-\kappa(1+\alpha\beta)\Omega^{2}A-(\alpha\beta+\beta+\kappa)\Omega B+(1+\beta)A - \frac{V_{\circ}}{R} \right\} \sin(\Omega \tau)= 0 \end{eqnarray} Bu aşamada yapılması gereken kilit gözlem şudur: sinus ve cosinus birbirlerinden lineer olarak bağımsız iki fonksiyon olduklarından, bu denklemdeki katsayıları -her ikisi de- sıfır olmak zorundadır. Böylece iki adet lineer denklem elde ediyoruz. \begin{eqnarray} \nonumber && (\alpha \beta \kappa \Omega^{3} - (\alpha \beta + \beta + \kappa) \Omega)A + (\kappa (1 + \alpha \beta)\Omega^{2} - (1+\beta))B = 0 \\ \nonumber &&-(\kappa (1 + \alpha \beta)\Omega^{2} - (1+\beta))A +(\alpha \beta \kappa \Omega^{3} - (\alpha \beta + \beta + \kappa) \Omega)B = \frac{V_{\rm o}}{R} \end{eqnarray} Bu denklemlerde cebirsel bataklığa dalıp boğulmak çok kolay. Bu yüzden önce aşağıdaki parametreleri tanımlayarak \begin{equation*} x := \alpha \beta \kappa \Omega^{3} - (\alpha \beta + \beta + \kappa) \Omega \ \ \ {\rm ve} \ \ \ y := \kappa (1 + \alpha \beta)\Omega^{2} - (1+\beta) \end{equation*} lineer denklem sistemini daha sade bir hale getiriyor \begin{equation*} xA+yB = 0 \ \ \ {\rm ve} \ \ \ -yA+xB = \frac{V_{\circ}}{R} \end{equation*} ardından da çözümünü veriyoruz. \begin{equation*} A = -\frac{y}{x^{2}+y^{2}} \frac{V_{\circ}}{R} \ \ \ {\rm ve} \ \ \ B = \frac{x}{x^{2}+y^{2}} \frac{V_{\circ}}{R} \end{equation*} $\varphi \in [0,2\pi)$ aralığında yer alan faz açısını $\cos \varphi := -y/\sqrt{x^{2}+y^{2}}$ ve $\sin \varphi = x/\sqrt{x^{2}+y^{2}}$ ile tanımlarsak, o zaman kısmi çözüm şöyle olacaktır: \begin{equation*} i^{\rm p}_{3}(\tau) = \frac{V_{\circ}/R}{\sqrt{x^{2}+y^{2}}} \sin(\Omega \tau + \varphi). \end{equation*} Bu denklemi yazarken $\sin(a+b)=\sin a \cos b + \cos a \sin b$ trigonometrik özdeşlikten faydalandık.

Hem kısmi hem de homojen çözümleri elde ettiğimize göre başlangıç değer probleminin toplam çözümünü verebiliriz. \begin{equation*} {\rm HP}: \ i_{3}(\tau) = i_{3}^{\rm h}(\tau) + i^{\rm p}_{3}(\tau) = C_{1} \exp(\lambda_{1}\tau) + C_{2} \exp(\lambda_{2}\tau) + C_{3} \exp(\lambda_{3}\tau) + \frac{V_{\circ}/R}{\sqrt{x^{2}+y^{2}}} \sin(\Omega \tau + \varphi) \end{equation*} Burada göstermeyeceğiz ama $C_{i}$ katsayılarını elde etmek için daha önce $i_{3}(0)$, $i^{\prime}_{3}(0)$ ve $i^{\prime \prime}_{3}(0)$ için verdiğimiz başlangıç şartlarından faydalanabiliriz. (Kuşkusuz tam çözümün bulunması kübik polinomun köklerinin hesaplanmasını da gerektirecektir.) Devremiz Routh-Hurwitz kararlı olduğundan, belli bir zaman sonra geçici (homojen) terimler etkisini yitirecek ve toplam çözüm asimptotik/pratik olarak $i^{\rm p}_{3}(\tau)$ ile aynileşecektir.

Geçirme bandı düz, düşük frekans bileşenlerini geçiren filtre

Her ölçüm cihazında olduğu gibi, spektroskopik ekipmanda da ya bir devre unsurundaki voltaj düşmesi ya da oradan geçen akım günlük dildeki kullanımıyla sinyal olarak addedilir ve bu sinyalin tefsiriyle incelenen malzemenin fizikokimyasal boyutundaki bir özelliğe ulaşılır. Diyelim ki spektrometre tasarımcımız $\beta R$ direncinden geçen akımı okumak istiyor. Bu akımın genliği devre parametrelerine ve girdi voltajının frekansına güçlü bir biçimde bağlılık gösterecektir. \begin{equation*} \max i_{3}(\tau) \sim \frac{V_{\circ}/R}{\sqrt{x^{2}+y^{2}}} \end{equation*} Daha açık yazıldığında $x^{2}+y^{2}$ ifadesi aşağıdaki gibidir. \begin{equation*} (\alpha \beta \kappa)^{2} \Omega^{6} + (\kappa^{2}(1+\alpha\beta)^{2}-2\alpha\beta\kappa(\alpha\beta+\beta+\kappa))\Omega^{4} + ((\alpha\beta+\beta+\kappa)^{2}-2\kappa(1+\alpha\beta)(1+\beta))\Omega^{2} + (1+\beta)^{2} \end{equation*} Geçici davranışın katkıları sönümlendikten sonra $i_{3}$ aşağıdaki gibi olur. \begin{eqnarray} \nonumber \max i_{3}(\tau) &\sim& \frac{V_{\circ}/R}{\alpha\beta\kappa \Omega^{3}}, \ \ \ \ \ \Omega \to \infty. \\ \nonumber \max i_{3}(\tau) &\sim& \frac{V_{\circ}/R}{1+\beta}, \ \ \ \ \ \ \Omega \to 0. \end{eqnarray} Açıkça görülüyor ki girdi voltajının frekansı eğer yüksek ise, o zaman tekabül eden sinyal $\beta R$ etrafında bastırılacağı gibi, frekansın düşük olması halinde akım genliği $V_{\circ}/(1+\beta)R$ olacak şekilde okunacaktır. İşte spektroskopide istediğimiz filtre türü tam da böylesi bir filtredir. Bu tip filtrelere literatürde düşük geçiren (lowpass) filtre denir ve $\Omega \ll 1$ şartını sağlayan frekansların bulunduğu bölge geçirme bandı (passband) olarak adlandırılır.

Devre parametrelerinin özensiz bir biçimde seçilmesi halinde, filtreden geçen akım genliği (passband bölgesinde) frekansa güçlü bir bağımlılık sergiler. Spektroskopicilerin deyişiyle okunan spektrumda yapay parlaklık/yoğunluk/şiddet (intensity) gözlenir. Akım genliğini tayin eden $(x^{2}+y^{2})^{-1/2}$ niceliğine transfer diyeceğiz. Yanda yer alan panoda iki ayrı parametre kümesi için transfer fonksiyonunun frekans davranışı sergilenmektedir. Bir parametre kümesi kullanıldığında arzulanan düz passband elde edilirken, ötekinde ise genlikte yapay dalgalanmalar bulunmaktadır. $x^{2}+y^{2}$ ifadesinin, $\Omega$ cinsinden yazıldığında 6. dereceden bir cebirsel fonksiyon olduğunu gözleyiniz. Passband bölgesini düzleştirmenin akıllı çözümlerinden birisi (ve naçizane kanaatime göre en güzeli) $x^{2}+y^{2}$ fonksiyonunun parametrelerinin, $n$ pozitif bir tam sayı ve $a$ pozitif bir gerçel sayı olmak üzere, $a+\Omega^{2n}$ formunu tutturacak şekilde ayarlanmasıdır. Bu durum ancak $\Omega^{2}$ ve $\Omega^{4}$ terimlerinin katsayılarının sıfır olmasıyla yakalanabilir. Dolayısıyla devre parametreleri olan $\alpha,\beta$ ve $\kappa$ aşağıdaki denklemleri sağlar ise \begin{eqnarray}\nonumber &&\kappa^{2}(1+\alpha\beta)^{2}-2\alpha\beta\kappa(\alpha\beta+\beta+\kappa)=0 \ \ \ {\rm ve} \\ \nonumber &&(\alpha\beta+\beta+\kappa)^{2}-2\kappa(1+\alpha\beta)(1+\beta)=0, \end{eqnarray} o zaman düz bir passband bölgesi garanti edilmiş olur. $\alpha=1$, $\beta=1$ ve $\kappa=2$ seçilerek bu denklemlerin gerçekten de sağlandığını not ediniz.

Daha fazla bilgi için

1. Sh. Strelitz, On the Routh-Hurwitz Problem, Amer. Math. Monthly 84, 542-544 (1977).
2. V. V. Prasolov, Polynomials, Springer-Verlag Berlin Heidelberg (2004).
3. S. Butterworth, On the Theory of Filter Amplifiers, Experimental Wireless and the Wireless Engineer 7, 536-541 (1930). [Okur, bu postada yer alan filtrelemeyle ilgili fikirlerin genelleştirildiği Butterworth'un çığır açıcı makalesine bakabilir.]

Strelitz algoritmasıyla Routh-Hurwitz kararlılık probleminin çözümü

Motivasyon

Fen bilimlerinde ve mühendislik uygulamalarında pek çok problem aşağıdaki lineer başlangıç değer problemine kadar indirgenebilir. \begin{equation*} c_{n} \frac{{\rm d}y^{n}(t)}{{\rm d}t^{n}} + \cdots + c_{1} \frac{{\rm d}y(t)}{{\rm d}t} + c_{0}y(t) = f(t) \end{equation*} Burada $y(t)$ çalıştığımız sistemi tarif eden bir durum değişkeni ve $c_{n},\ldots,c_{0}$ ise sabit değerli (ve genellikle gerçel) katsayılardır. $f(t)$ ise sistemin davranışını yönlendirmeye çalışan ve dışarıdan uygulanan bir çeşit itici güç (driving force) olarak adlandırılır ve genellikle amaç sistemin bu itici güce verdiği tepkiyi (response) $y$ niceliği kanalıyla hesaplamaktır. Örneğin, ne kadar karmaşık olursa olsun, pasif, lineer devre elemanlarıyla kurulan ve sözgelimi sinusoidal bir güç kaynağı ile tahrik edilen bütün analog devrelerin dinamiği bu şablona uyar. $f$, her $t$ için sıfır değerini aldığında bu denklemin homojen olduğu söylenir ve çözümü $y_{\rm h}(t)=\exp(\lambda t)$ şablonunda fonksiyonlarla aranır. İlgili homojen denkleme bu deneme fonksiyonunu koyup gerekli sadeleştirmeleri yaptığımızda karşımıza cebirsel bir denklem çıkar. \begin{equation*} c_{n}\lambda^{n} + \cdots + c_{1} \lambda + c_{0} = 0 \end{equation*} Demek ki aradığımız deneme fonksiyonlarındaki $\lambda$ değerleri, derecesi $n$ olan cebirsel bir denklemin kökleriymiş. Buraya kadar her şey son derece tek düze giderken, özellikle kararlılık analizi aşamasında $t \to \infty$ asimptotiğinde $|y_{\rm h}| < y_{\rm max} < \infty$ olmasını hatta $|y_{\rm h}| \to 0$ olmasını isteriz. Sistem tasarımcılarının en büyük kabusu $t \to \infty$ asimptotiğinde $|y| \to \infty$ davranışına rastlamaktır. (Devre elemanından sonsuz akım geçmesi o devrenin -eğer patlamazsa- aşırı ısınıp buharlaşması anlamına gelir.) $t \to \infty$ asimptotiğinde kararlı bir davranış elde etmemizin, deneme fonksiyonlarının şablonuna bakıldığında $\Re [\lambda_{i}] < 0$ şartında düğümlendiği barizdir. Burada $\Re[z]$ ile karmaşık bir sayının gerçel kısmını simgeliyoruz. Bu ve benzeri endişeler Routh-Hurwitz kararlılık problemini tanımlar.

Routh-Hurwitz problemi: Katsayıları gerçel bir polinomun tüm köklerinin gerçel kısımlarının negatif olması için gerekli ve yeterli şartları bulunuz.

Strelitz çözümü

Routh-Hurwitz problemi, cebirsel kök tasnifi (root classification) konusunun abecesini oluşturur ve bu alanda öğrencinin ilk rastladığı konulardan birisidir. Bu problem 19. yy'da ortaya atılmış ve Sturm teoremi kullanılarak çözülmüştür. 20. yy'da ise Strelitz tarafından Sturm teoremini kullanmadan, daha düşük cebirsel işlem hammaliyesiyle başka bir çözüm önerilmiştir. Biz bu postada Strelitz'in çözümünü orijinal makalesini takip ve yer yer şerh ederek sunmaya çalışacağız.

Polinomun kökleri, polinomdaki bütün katsayıları baş katsayıya böldüğümüzde değişmeyeceğinden, genelliği kaybetmeden, aşağıdaki monik polinomu ele alıyoruz. \begin{equation} P(z) := z^{n} + a_{1}z^{n-1} + \cdots + a_{n} \ \ \ \ \ (1) \end{equation} Katsayı indisleri ile terimlerin dereceleri arasındaki zıt yönlerde artışa dikkat ediniz. (Literatürdeki çoğu formül maalesef bizim burada uyguladığımız konvansiyonu takip etmiyor ama yazarla paralel gidebilmek için biz bu konvansiyonu tercih etmek zorunda kaldık.) Strelitz bu noktada bariz olduğunu düşünerek ispatsız şöyle bir argüman sunuyor:

(1) nolu denklemdeki katsayıların gerçel olmaları şartının zaruri olmadığı barizdir, zira eğer karmaşık katsayılı bir $P_{0}(z)$ polinomunun bütün köklerinin gerçel kısımları negatif ise, o zaman $P_{0}(z)\overline{P_{0}(\overline{z})}$ polinomunun tüm köklerinin gerçel kısımları da negatiftir.

(Burada $\overline{z}$ ile $z$ sayısının karmaşık eşleniği kastediliyor.) Bu önermeyi ispatlayalım. Monik $P_{0}$ polinomunun derecesi $k$, kökleri ise $z_{1},\ldots,z_{k}$ olsun. O zaman cebirin temel teoremi uyarınca $P_{0}(z) = (z-z_{1}) \cdots (z-z_{k})$ yazabiliriz. Yine benzer şekilde $P_{0}(z)\overline{P_{0}(\overline{z})}$ polinomunun köklerinin de $z_{1},\overline{z_{1}}, \ldots , z_{k}, \overline{z_{k}}$ olması gerektiği barizdir. Dolayısıyla cebirin temel teoremini kullanarak \begin{eqnarray}\nonumber P_{0}(z)\overline{P_{0}(\overline{z})} &=& (z-z_{1}) (z-\overline{z_{1}}) \cdots (z-z_{k}) (z-\overline{z_{k}}) \\ \nonumber &=& (z^{2} - 2 \Re [z_{1}]z + |z_{1}|^{2}) \cdots (z^{2} - 2 \Re [z_{k}]z + |z_{k}|^{2}) \end{eqnarray} sonucunu elde ederiz. Bu denklemin sağ tarafının sadece gerçel katsayılardan oluştuğu barizdir. Uzun lafın kısası Strelitz şunu demeye çalışıyor: Eğer $P_{0}$ polinomunun katsayıları karmaşık ise, o zaman bu polinomun Routh-Hurwitz sorgusu yerine katsayıları tamamen gerçel olan $P_{0}(z)\overline{P_{0}(\overline{z})}$ polinomunun Routh-Hurwitz sorgusuna bakabiliriz. Her iki polinomun da köklerinin gerçel kısımları aynı olduğundan, Routh-Hurwitz sorguları da aynı sonucu verecektir.

(1) nolu denklemde verilen polinoma ek olarak ikinci bir polinom tanımlıyoruz. \begin{equation} Q(z) := z^{m} + b_{1}z^{m-1} + \cdots + b_{m}, \ \ \ m := \frac{n(n-1)}{2} \ \ \ \ \ (2) \end{equation} (Strelitz'in makalesinde bu polinomun baş teriminin derecesi yazılırken baskı hatası yapılmış, yanlışlıkla $n$ konulmuş.) $P$ polinomun kökleri $z_{k}$ (ve $k=1,\ldots,n$) ise, o zaman $Q$, kökleri \begin{equation*} z_{i}+z_{j}; \ \ \ i < j, \ \ \ i,j=1,\ldots,n \end{equation*} ile verilen polinom olacak şekilde tanımlanıyor. Makalenin ilerleyen bölümlerinde, $P$'nin köklerini hesaplamadan $Q$ polinomunu türetmenin yollarını gösteriyor yazar. Ayrıca $n$ adet kökten seçilebilecek ikili çiftlerin sayısı da $m = n(n-1)/2$ eşitliğini izah etmektedir. Bu ön hazırlıktan sonra Routh-Hurwitz probleminin kesin çözümünü verebiliriz.

Teorem: Katsayıları gerçel $P(z)$ polinomunun tüm köklerinin karmaşık düzlemin sol yarısında kalması için gerek ve yeter şart hem $P(z)$ hem de $Q(z)$ polinomlarının katsayılarının pozitif olmasıdır. Diğer bir deyişle \begin{eqnarray} \nonumber &&a_{j}>0, \ \ \ j=1,\ldots, n; \\ \nonumber &&b_{k}>0, \ \ \ k=1,\ldots, m; \ \ \ m = n(n-1)/2. \end{eqnarray} İspat: (Strelitz) Gereklilik. $P$'nin bütün köklerinin gerçel kısmı karmaşık düzlemin sol tarafında yer alsın. Bu köklere $z_{1},\ldots,z_{n}$ diyelim. Söz konusu polinomun katsayıları gerçel olduğundan $P(z_{i})=0 \iff P(\overline{z_{i}})=0$ olduğu barizdir. $z_{p}=\alpha+\beta i$ olsun. (Hipotez gereği $\alpha < 0$ olduğunu unutmayınız.) O zaman \begin{equation*} (z-z_{p})(z-\overline{z_{p}}) = z^{2} - 2\alpha z + \alpha^{2} +\beta^{2} \ \ \ \ \ (3) \end{equation*} ifadesindeki bütün katsayılar pozitif olur. $z_{q}$ gerçel bir kök olsun. Yine hipotez gereği $z_{q}<0$ olması gerektiğinden, $(z-z_{q})$ ifadesindeki bütün katsayılar yine pozitif olacaktır. $k=1,\ldots,s$ için $\alpha_{k}+i\beta_{k}$ ve $\alpha_{k}-i\beta_{k}$, $P$ polinomunun tüm karmaşık köklerini temsil ederken, $j=1,\ldots,t=n-2s$ için $\gamma_{j}$ aynı polinomun gerçel köklerini temsil etsin. O zaman \begin{equation*} P(z) = \prod_{k=1}^{n} (z-z_{k}) = \prod_{k=1}^{s} (z^{2} - 2\alpha _{k} z + \alpha _{k}^{2} +\beta _{k}^{2}) \prod_{l=1}^{t} (z-\gamma_{l}) \ \ \ \ \ (4) \end{equation*} yazabiliriz. (4) nolu denklemdeki ifadeleri çarptıktan sonra $-\alpha_{k}$ ve $-\gamma_{k}$ değerlerinin negatifliğine ilişkin yukarıdaki mülahazalar muvacehesinde, (1) nolu denklemdeki bütün katsayıların pozitif olduğu sonucuna varabiliriz: $a_{i} > 0$.

Ayrıca, $k=1,2,\ldots, n$ için $\Re [z_{k}] < 0$ olması halinde, o zaman $Q(z)$ polinomunun bütün kökleri karmaşık düzlemin sol tarafında kalır ve böylece $k=1,2,\ldots,m$ için $b_{k} > 0$ olduğunu da göstermiş oluruz.

Yeterlilik. $\alpha \pm i \beta$, $P(z)$ polinomunun $\alpha$ ve $\beta$ sıfırdan farklı olacak şekilde iki kökü olsun. O zaman $2\alpha$ da $Q(z)$ polinomunun bir köküdür: $Q(2\alpha) = 0$. Ancak, $Q(z)$ polinomunun bütün katsayıları, $b_{k}$, pozitif olduğundan, $Q(z)=0$ denkleminin tüm gerçel kökleri negatiftir ve sonuç itibariyle $P(z)=0$ denkleminin tüm karmaşık kökleri karmaşık düzlemin sol yarısına düşer. Ayrıca, $P(z)$ polinomunun da tüm gerçel kökleri gerçel eksenin sol tarafındadır, zira $j=1,2,\ldots,n$ için $a_{j} > 0$.

Teoremin ispatı burada bitmiştir.

Strelitz algoritması

Yukarıdaki ispat güzel ama işe yarayabilmesi için $Q(z)$ polinomunun katsayılarını, $P(z)$ polinomunun köklerini hesaplamadan, hızlı bir biçimde hesaplayacak bir yönteme ihtiyacımız var. Makalenin geri kalan kısmında Strelitz ilgili katsayıları hesaplamak için Newton toplamlarına dayanan bir algoritma kurmakla uğraşıyor. Bu amaçla $k=1,\ldots,n$ olmak üzere $P(z)$ polinomunun köklerine $z_{k}$ diyelim ve \begin{equation*} \sigma_{j} := \sum_{k=1}^{n} z_{k}^{j}; \ \ \ s_{j} := \sum_{p=1}^{n-1} \sum_{q=p+1}^{n} (z_{p}+z_{q})^{j}, \ \ \ j=0,1,2,3,\ldots \end{equation*} niceliklerini tanımlayalım. Bu tanımların her ikisi de bildiğiniz Newton toplamı: ilki $P$, ikincisi $Q$ polinomlarının kökleriyle hesaplanmış. $j \geq n$ için Newton toplamını ilk $n-1$ Newton toplamı verilmişse hesaplamak çok kolay. Zira $0=P(z_{k})=z_{k}^{n}+a_{1}z_{k}^{n-1}+\cdots+a_{n}$ eşitliğinde her iki tarafı da önce $z_{k}^{j-n}$ ile çarpıp, sonra bütün kökler üzerinden toplama yaparsak \begin{equation*} \sigma_{j} + a_{1}\sigma_{j-1} + \cdots + a_{n}\sigma_{j-n} = 0, \ \ \ j \geq n, \ \ \ (5a) \end{equation*} denklemine ulaşıyoruz. (Strelitz'in makalesinde yine bir baskı hatası var. Bu denklemdeki $a_{n}$ terimi unutulmuş. bak. 544. sayfadaki ilk denklem.) Bu denklem sayesinde kendisinden önceki $n$ Newton toplamının verilmesi halinde bir sonraki Newton toplamının da sadece çarpma ve toplama işlemlerini kullanarak hesaplanabileceğini gözleyiniz. $j < n$ için Newton toplamlarının nasıl hesaplanılacağı konusunda biraz düşünmek gerekiyor. Gerçi temel cebirsel özdeşlikler kanalıyla $\sigma_{0}=n$ ve $\sigma_{1}=-a_{1}$ (basitçe kökler toplamı) olduğu kolayca ortaya çıksa da diğer Newton toplamları o kadar kolay değil. Biz ispatını daha sonra yapmak üzere sadece bu durum için Newton toplamlarının formülünü vermekle yetineceğiz. \begin{equation*} \sigma_{j} + a_{1}\sigma_{j-1} + \cdots + a_{j-1}\sigma_{1} + j a_{j} = 0, \ \ \ j < n, \ \ \ (5b) \end{equation*}

Algoritmanın ikinci aşamasında $s_{j}$ niceliklerini $\sigma_{j}$ cinsinden ifade etmeye çalışacağız. Bu amaçla \begin{equation*} \left( \sum_{k=1}^{n} e^{z_{k}t} \right)^{2} = e^{2z_{1}t}+e^{2z_{2}t}+ \cdots +e^{2z_{n}t} + 2 \sum_{p=1}^{n-1}\sum_{q=p+1}^{n} e^{(z_{p}+z_{q})t} \end{equation*} özdeşliğinin her iki tarafını da Taylor serilerini kullanarak açacağız. (Strelitz'in makalesinde yine bir baskı hatası var. 544. sayfada yukarıdan dördüncü satırda yer alan bu denklemde $e^{2z_{n}t}$ terimi yerine $e^{2z_{k}t}$ yazılmış.) Sol tarafı açtığımızda \begin{eqnarray}\nonumber \left( \sum_{k=1}^{n} e^{z_{k}t} \right)^{2} &=& \left( \sum_{j=0}^{\infty} \sum_{k=1}^{n} z_{k}^{j} \frac{t^{j}}{j!} \right)^{2} = \left( \sum_{j=0}^{\infty} \sigma_{j} \frac{t^{j}}{j!} \right)^{2} = \sum_{j=0}^{\infty}\sum_{i=0}^{\infty} \sigma_{j}\sigma_{i} \frac{t^{i+j}}{i!j!} \\ \nonumber &=& \frac{\sigma_{0}\sigma_{0}}{0!0!} + \left( \frac{\sigma_{1}\sigma_{0}}{1!0!} + \frac{\sigma_{0}\sigma_{1}}{0!1!} \right) t + \left( \frac{\sigma_{2}\sigma_{0}}{2!0!} + \frac{\sigma_{1}\sigma_{1}}{1!1!} + \frac{\sigma_{0}\sigma_{2}}{0!2!}\right) t^{2} + \cdots \\ \nonumber &=& \sum_{j=0}^{\infty} \left( \sum_{i=0}^{j} \frac{\sigma_{i}\sigma_{j-i}}{i!(j-i)!} \right) t^{j} \end{eqnarray} eşitliği elde ediliyor. Aynı özdeşliğin sağ tarafını açtığımızda ise \begin{eqnarray}\nonumber \sum_{k=1}^{n} e^{2z_{k}t} + 2 \sum_{p=1}^{n-1}\sum_{q=p+1}^{n} e^{(z_{p}+z_{q})t} &=& \sum_{j=0}^{\infty}\sum_{k=1}^{n} z_{k}^{j} \frac{2^{j}t^{j}}{j!} + 2 \sum_{j=0}^{\infty} \sum_{p=1}^{n-1}\sum_{q=p+1}^{n} (z_{p}+z_{q})^{j} \frac{t^{j}}{j!} \\ \nonumber &=& \sum_{j=0}^{\infty} \sigma_{j} \frac{2^{j}t^{j}}{j!} + 2 \sum_{j=0}^{\infty} s_{j} \frac{t^{j}}{j!} \end{eqnarray} sonucunu elde ediyoruz. Matematikte iş bu aşamaya geldiğinde genellikle denklemin sol ve sağ tarafında yer alan aynı kuvvetteki terimlerin ($t^{j}$) katsayıları eşitlenir. (Neden?) Yeniden düzenleme yaptığımızda $Q$ polinomunun Newton toplamlarını $P$ polinomunun Newton toplamlarına bağlayan özdeşliği elde etmiş oluyoruz. \begin{equation*} s_{j} = -2^{j-1}\sigma_{j} + \frac{1}{2} \sum_{i=0}^{j} \frac{j!}{i!(j-i)!} \sigma_{i}\sigma_{j-i} \end{equation*} (Sterlitz'in makalesinde bu denklemle ilgili yine bir işlem hatası var. Onun verdiği formülü uyguladığınızda $s_{0}=n(2n-1)/4$ çıkıyor. Bu sonuç yanlış. Bizim burada verdiğimiz formül kullanıldığında $s_{0}=n(n-1)/2$ elde ediliyor, olması gerektiği gibi. Zira derecesi sıfır olan Newton toplamı her zaman polinomun derecesine eşittir.)

Algoritmanın üçüncü ve son aşamasında $Q$ polinomunun katsayılarını elde etmek var. Bu ise çok kolay. Öncelikle $b_{0}=1$ olduğunu not edelim. Sonra (5) nolu denklemler kanalıyla, bu denklemlerde $a$ yerine $b$, $\sigma$ yerine $s$ ve $n$ yerine $m$ koyarak, önce $s_{1} + b_{1} = 0$ ilişkisinden $b_{1}=-s_{1}$ hesaplanabileceğini gözleyelim. Ardından $s_{2} + b_{1}s_{1} + 2b_{2}=0$ denkleminin bize $b_{2}$ değerini vereceğini. Böyle devam ettiğimizde bütün $b_{j}$ katsayılarının hesaplanacağı da barizdir.

Çift Atwood makinesi

Problem: Şekilde $m_{1}$ ve $m_{2}$ kütleleri, kütlesiz, sürtünmesiz ve esnemez, kısaca ideal, $A$ sicimiyle hareketli (yine kütlesiz ve sürtünmesiz) $B$ makarası üzerinden birbirlerine bağlıdırlar. $B$ makarası ise ideal $C$ sicimiyle kütlesiz ve sürtünmesiz, tavana sabitlenmiş, $D$ makarası ile $m_{3}$ kütlesine bağlıdır. Sistem başlangıçta hareketsiz olarak $t=0$ anında resimde yukarıdan aşağıya olan $g$ yerçekimi ivmesi ile harekete başlamaktadır. Tüm kütlelerin ve hareketli $B$ makarasının ivmesiyle, $A$ ve $C$ sicimlerindeki gerilmeyi soruda verilenler cinsinden ifade ediniz. (Resimin üzerine tıklayarak daha yüksek çözünürlükteki versiyonunu görebilirsiniz.)

Bu soruyu bir fizik ders kitabının arkasındaki sorular arasında işte meydan sınıfında bir soru olarak gördüm ve epeyce ilgimi çekti. Küçücük (hem de ideal) bir basit makinenin insanı nasıl otuzüçleteceğini görmek isterseniz, bu soruda şansınızı denemelisiniz.

Soruyu çözmeye başlamadan önce bizden istenilenleri, yani bilinmeyenleri, yazarak en az kaç adet denkleme ihtiyacımız olduğunu tespit edelim. Kütlelerin ivmeleri $a_{1}$, $a_{2}$ ve $a_{3}$ ve hareketli $B$ makarasının ivmesi $a_{B}$ olmak üzere 4 adet ivme. Buna ek olarak $A$ ve $C$ sicimlerindeki gerilimler $T_{A}$ ve $T_{B}$ var. Dolayısıyla toplamda 6 adet bilinmeyenimiz mevcut. Bunun için altı tane denklem ya da şart yazmamız gerekiyor.

Şimdi eğer Newton'ın ikinci hareket kanununu üç kütleye ve kütlesiz makaraya uygularsak, o zaman elimizde dört tane denklem olmaktadır. Bu ise yeterli değil. O zaman fiziğin bize sağladığı denklemlerin/şartların dışında geometrik ya da matematiksel başka şartlar bulunmalı ki, problemin çözümü tespit edilebilsin. Bizim çalıştığımız problemde sicimlerin uzunluğunun sabit olması bir geometrik şart olarak ortaya çıkmaktadır. Ben cisimlerin $D$ makarasının merkezine göre olan uzaklıklarını sırasıyla $x_{1},$ $x_{2},$ $x_{3}$ ve $x_{B}$ nicelikleri ile ölçtüm. Kuşkusuz başka bir eylemsiz referans noktasına göre de, örneğin tavana göre, bu ölçümleri yapabilirdik. Ama $B$ makarasının merkezine göre bu ölçümleri yapsaydık çıkan sonuç, bu makara ivmelendiği için yanlış olacaktı.

$B$ ve $D$ makaralarının sabit olan çevre uzunluklarına sırasıyla $s_{B}$ ve $s_{D}$ diyelim. O zaman basitçe şekile bakarak $A$ siciminin uzunluğunun $l_{A} = (x_{1}-x_{B}) + s_{B}/2 + (x_{2}-x_{B})$ ya da \begin{equation*} l_{A} = x_{1} - 2x_{B} + x_{2} + \frac{1}{2} s_{B} \end{equation*} olduğunu tespit ederiz. Resime göre aşağı yönü pozitif aldığımızdan, kinematiğin temel tanımlarını kullanarak bu denklemin zamana göre ikinci türevini alırsak, o zaman \begin{equation*} {\rm K1:} \ \ \ \boxed{a_{1} - 2a_{B} + a_{2} = 0} \end{equation*} denklemini elde ediyoruz. Tamamen benzer şekilde $C$ siciminin uzunluğu için de \begin{equation*} l_{C} = x_{B} + x_{3} + \frac{1}{2}s_{D} \end{equation*} denklemi geçerlidir. Yine her iki tarafın zamana göre ikinci türevi alındığında \begin{equation*} {\rm K2:} \ \ \ \boxed{a_{B}+a_{3} = 0} \end{equation*} denklemini çıkarmış oluyoruz.

Problemin zor noktasını aştık. Bundan sonra artık sadece Newton'ın ikinci hareket kanununu sistemde hareket eden dört cisme de uygulayacağız ve mekanik içerikli dört denklem elde edeceğiz. (Aşağı yönü pozitif aldığımızı tekrar vurgulayalım.) \begin{eqnarray} \nonumber &&{\rm M1:} \ \ \ \boxed{m_{1}g - T_{A} = m_{1}a_{1}} \\ \nonumber &&{\rm M2:} \ \ \ \boxed{m_{2}g - T_{A} = m_{2}a_{2}} \\ \nonumber &&{\rm M3:} \ \ \ \boxed{m_{3}g - T_{C} = m_{3}a_{3}} \\ \nonumber &&{\rm M4:} \ \ \ \boxed{2T_{A} - T_{C} = 0} \end{eqnarray} Bu mekanik denklemlerden ilk üçü aşırı derecede kolay anlaşılabilir. Lakin dördüncü mekanik denklemin sağ tarafının sıfır olması denge şartından kaynaklanmamakta, tam tersine hareketli olan $B$ makarasının kütlesinin sıfır olmasından kaynaklanmaktadır. $B$ makarasının kütlesi olsaydı, o zaman M4 etiketli denklemin sağ tarafı $m_{B}a_{B}$ ile verilecekti.

Soruyla ilgili yazdığımız denklemlerde bütün bilinmeyenler birinci dereceden yer alıyor. Bu tip denklemlere literatürde lineer denklemler denir. Lineer cebirde Gauss eliminasyonu denilen bir yöntem ile bütün lineer denklemler -eğer çözümleri varsa- sadece aritmetik işlemler kullanılarak, sonlu sayıda işlem sonrasında çözülebilirler. Bu kuşkusuz bilgisayarlar için, hele hele sadece altı tane bilinmeyeni olan bir sistemde, çocuk oyuncağıdır ama iş Gauss eliminasyonunu kağıt kalemle uygulamaya geldiğinde bu biraz kıllı... Denklem sistemimizdeki altı denklemin hemen hemen hepsinde sadece iki değişken yer alıyor. Normalde altı değişken/bilinmeyen de yer alabilir. Bu tip az sayıda değişkenin kendini denklemlerde gösterdiği sistemlere literatürde sparse (seyrek) sistemler denir. Seyrek sistemlerin çözümleri nisbeten daha hızlıdır. Çalıştığımız sistem seyrek olduğu için onu kağıt kalemle çözmeye cüret edebiliriz. Normalde doğrudan bir bilgisayar cebir yazılımını (mesela Maxima) kullanmamız gerekirdi bu aşamada.

Şimdi (K1) etiketli denklemden $a_{B} = (a_{1}+a_{2})/2$ ve (K2) etiketli denklemden de $a_{3} = -a_{B} = -(a_{1}+a_{2})/2$ ara sonuçlarını elde edebiliriz. Devamla (M1) ve (M2) denklemlerini sırasıyla $m_{1}$ ve $m_{2}$ ile böler ve çıkan sonuçları taraf tarafa toplarsak, o zaman aşağıdaki ara sonuca ulaşırız. \begin{equation*} {\rm AS1:} \ \ \ 2g - \frac{m_{1}+m_{2}}{m_{1}m_{2}}T_{A} = a_{1}+a_{2} = -2a_{3} \end{equation*} Şimdi (AS1) ile (M3) etiketli denklemler taraf tarafa bölünür ve (M4) etiketli denklemden $T_{C}=2T_{A}$ yazılırsa, o zaman \begin{equation*} {\rm AS2:} \ \ \ \frac{m_{3}g-2T_{A}}{2g - \frac{m_{1}+m_{2}}{m_{1}m_{2}}T_{A}} = -\frac{m_{3}}{2} \end{equation*} denklemine ulaşıyoruz. Sadece bir tane bilinmeyen içeren (AS2) etiketli denklem o kadar basit ki ben bile çözebilirim! \begin{equation*} T_{A} = \frac{4m_{1}m_{2}m_{3}g}{4m_{1}m_{2}+m_{1}m_{3}+m_{2}m_{3}} \end{equation*} $T_{C}$ ise bu gerilimin tam iki katı kadar olacaktır. (Gerilim değerlerinin, olması gerektiği gibi, çalıştığımız referans sisteminde her türlü kütle seçimi için pozitif çıktığına dikkat ediniz.) (M1) etiketli denklem yoluyla \begin{equation*} a_{1} = g - \frac{1}{m_{1}}T_{A} = \frac{4m_{1}m_{2}+m_{1}m_{3}-3m_{2}m_{3}}{4m_{1}m_{2}+m_{1}m_{3}+m_{2}m_{3}}g \end{equation*} hesaplanılır. Benzer şekilde (M2) etiketli denklem bize \begin{equation*} a_{2} = g - \frac{1}{m_{2}}T_{A} = \frac{4m_{1}m_{2}-3m_{1}m_{3}+m_{2}m_{3}}{4m_{1}m_{2}+m_{1}m_{3}+m_{2}m_{3}}g \end{equation*} formülünü verecektir. Son olarak \begin{equation*} a_{3} = -\frac{1}{2}(a_{1}+a_{2}) = \frac{4m_{1}m_{2}-m_{1}m_{3}-m_{2}m_{3}}{4m_{1}m_{2}+m_{1}m_{3}+m_{2}m_{3}}g \end{equation*} ifadesi hesaplanılır. $a_{B} = -a_{3}$ olduğundan problemin çözümü burada bitmiştir.

İki boyutta noktasal tanecikler arası esnek çarpışmanın trigonometrisi

Kütlesi $m$, hızı $+x$ yönünde $v_{\circ}$ olan bir cisim başlangıçta ilgili referans/koordinat sisteminde duran bir tanecikle esnek olarak çapışıyor. Saçılmadan sonra hareketin bir düzlem üzerinde kalacağı barizdir. (Neden?) Saçılmanın geometrisi öyle ki $m$ cismi $+x$ ekseni ile $\alpha$ açısı yapan bir $v$ hızı ile saçılırken, kütlesi $M$ olan öteki cisim ise yine aynı eksenle $\beta$ açısı yapan bir $w$ hızıyla saçılıyor. Amacımız $\alpha$ ve $\beta$ arasında trigonometrik bir bağıntı bulmak.

Saçılma ister elastik (esnek) isterse inelastik (esnek olmayan) bir biçimde gerçekleşsin, $m$ ve $M$ kütlelerinden oluşan sisteme dışarıdan bir kuvvet etki etmeyeceği için, toplam momentum vektörü korunacaktır. Hareket iki boyutta gerçekleştiği için bu bize iki tane denklem sunar. Öncelikle toplam momentumun $x$ bileşeninden başlayalım. \begin{equation*} mv_{\circ} = mv \cos \alpha + Mw \cos \beta \end{equation*} Bu denklemde görünüşte iki parametre var. Ama her iki tarafı $m$ ile böler ve $\kappa := M/m$ oranını tanımlarsak, o zaman denklemimizdeki parametre sayısını ikiden bire düşürmüş oluruz. \begin{equation*} \boxed{ v_{\circ} = v \cos \alpha + \kappa w \cos \beta } \end{equation*} İkinci olarak $y$ yönündeki toplam momentum vektörünü yazalım. \begin{equation*} 0 = mv \sin \alpha - M w \sin \beta \end{equation*} Her iki tarafı $m$ ile bölüp $\kappa$ tanımını uygularsak aşağıdaki denkleme ulaşırız. \begin{equation*} \boxed{0 = v \sin \alpha - \kappa w \sin \beta} \end{equation*} Kutu içindeki her iki formülde de her iki tarafında karelerini alıp taraf tarafa toplama yapar ve $\sin^{2} \varphi + \cos^{2} \varphi =1$ ve $\cos (\varphi + \vartheta) = \cos \varphi \cos \vartheta - \sin \varphi \sin \vartheta$ trigonometrik özdeşliklerinden faydalanırsak o zaman aşağıdaki nisbeten daha sade ara sonuca (as1) ulaşırız. \begin{equation*} {\rm as1:} \ \ \ v_{\circ}^{2} = v^{2} + \kappa^{2} w^{2} + 2 \kappa v w \cos (\alpha + \beta) \end{equation*}

Çarpışma elastik olduğu için toplam kinetik enerji de korunacaktır. O zaman \begin{equation*} \frac{m}{2}v_{\circ}^{2} = \frac{m}{2}v^{2} + \frac{M}{2}w^{2} \end{equation*} denklemi geçerli olur. Yine her iki tarafı $m/2$ ile bölüp $\kappa$ tanımını uygularsak, o zaman \begin{equation*} \boxed{v_{\circ}^{2} = v^{2} + \kappa w^{2}} \end{equation*} denklemini elde ediyoruz.

Bu noktaya değin yapageldiğimiz analizde bir şey dikkatimizi çekmeli. Problemimizin bilinmeyenleri $\alpha,\beta,v,w$ olmak üzere dört adet. Öte yandan biz kutu içindeki denklemleri, ki fiziksel prensiplerden türettik bunları, saydığımızda üç adet denklem yazabiliyoruz. Varacağımız yargı kaçınılmaz olarak iki boyutta, iki noktasal cismin elastik çarpışmasına ait problemin çözülemeyeceğidir. Problemin çözümü derken, bilinmeyenlerin tamamen bilinen nicelikler cinsinden ifade edilmesini kastediyoruz.

Şimdi (as1) etiketli denklem ile kutu içindeki üçüncü denklemi taraf tarafa çıkartır ve sonucu yeniden düzenlersek \begin{equation*} \cos (\alpha + \beta) = \frac{1-\kappa}{2 \kappa} \frac{w}{v} \end{equation*} bağıntısına ulaşırız. Burada eğer kütleler çok özel olarak $m=M$ şeklinde seçildiğinde, örneğin nötron-nötron çarpışması gibi, $\kappa=1$ ve $\cos (\alpha+\beta)=0$ sonucu çıkacaktır. Bu ise $\alpha + \beta = \tfrac{\pi}{2}$ olmasını gerektirir. Demek ki özdeş noktasal taneciklerin elastik çarpışmalarında saçılma açısı her zaman $90^{\circ}$ oluyormuş. Bizim pu postada amacımız bu durumu genelleştirmek.

Şimdi ikinci kutudaki denklem yeniden düzenlenirse, o zaman $w/v = \sin (\alpha ) / \kappa \sin (\beta)$ ara sonucuna ulaşılır. Bu ara sonuç $\cos (\alpha + \beta)$ için verilen denklemde yerine konursa aşağıdaki ifadeye ulaşırız. \begin{equation*} \cos(\alpha + \beta) = \frac{1-\kappa}{2\kappa} \frac{\sin \alpha}{\sin \beta} \end{equation*} Artık bu aşamadan sonra amacımız $\alpha$ ve $\beta$ arasında daha net bir ilişki türetilmek için trigonometri bilgimizi kullanmak. \begin{eqnarray} \nonumber \cos \alpha \cos \beta - \sin \alpha \sin \beta &=& \frac{1-\kappa}{2\kappa} \frac{\sin \alpha}{\sin \beta} \\ \nonumber \frac{1}{2} \cot \alpha \sin (2\beta) - \sin^{2} \beta &=& \frac{1-\kappa}{2\kappa} \\ \nonumber \cot \alpha \sin (2\beta) + \cos (2\beta) &=& \frac{2-\kappa}{2\kappa} \\ \nonumber \cos \alpha \sin (2\beta) + \sin \alpha \cos (2\beta) &=& \frac{2-\kappa}{2\kappa} \sin \alpha \\ \nonumber \sin(2\beta + \alpha) &=& \frac{2-\kappa}{2\kappa} \sin \alpha \end{eqnarray} Son olarak her iki tarafa da arcsin fonksiyonunu uygularsak, saçılma açılarından birisini ötekisi cinsinden ifade etmiş olur ve çözümü tamamlarız. \begin{equation*} \beta = -\frac{\alpha}{2} + \frac{1}{2} \arcsin \left( \frac{2-\kappa}{2\kappa} \sin \alpha \right) \end{equation*}

A rigorous explanation of Rüchhardt's approximation

The ingenious method of Rüchhardt (1929) determines the ratio of heat capacities of an ideal gas through elementary mechanics as an alternative to the well established methods of thermal physics. One can find the description of the method elsewhere on the net as well as in Zemansky's influential Heat and Thermodynamics. I was not satisfied by the exposition of the method in Zemansky and decided to present it in my own way.

Adiabats

The method relies neither the first nor the second laws of thermodynamics. We start with The Doctrine of Latent and Specific Heats which states that the instant heating of a thermodynamic system can be described by the following equation. \begin{equation} h(t) := \Lambda_{V}(T(t),V(t)) \dot{V}(t) + C_{V}(T(t),V(t)) \dot{T}(t) \label{doctrine} \end{equation} Here $\dot{x}(t)$ signifies the time derivative of $x(t)$, $\Lambda_{V}(T,V)$ is the latent heat function of the material, and the positive definite function $C_{V}(T,V)$ is the specific heat at constant volume. $t$, as always signifies time. All physical theories rely on the existence of an entity called state and in thermodynamics we assume that the state of a closed sytem can be described by two independent variables, say temperature and volume. The collection of all states is given by the set $\Sigma$ and we assume that $\Sigma$ is simply connected. A thermodynamic process can be represented by a (piecewise) differentiable curve in $\Sigma$. Thus the time parametric function $\mathbf{\alpha} : [a,b] \to \Sigma$ gives all the points $(T(t),V(t))$ of the path traced by a thermodynamic process as time runs from $a$ to $b$.

Among all processes, adiabatic ones stand out for two reasons: They are indispensible in the description of Carnot's heat engine as well as the experiments of Joule which ultimately gave us the second and first laws of thermodynamics, respectively. A process is adiabatic if, and only if, the instant heating function vanishes at all times. Set $h(t)=0$ and rearrange (\ref{doctrine}) to read \begin{equation} \frac{{\rm d}\psi(V)}{{\rm d}V} = -\frac{\Lambda_{V}(\psi(V),V)}{C_{V}(\psi(V),V)}. \label{adiabat} \end{equation} The solutions of the initial value problem (\ref{adiabat}) are called adiabats, they are signified by $T = \psi(V)$, and satisfy the initial condition $T_{\circ} = \psi(V_{\circ})$. The uniqueness of the solutions of ordinary differential equations guarantee that two adiabats either agree on all points in $\Sigma$ or they never cross.

A second hidden postulate of the founding fathers of thermodynamics was their relience on the existence of a pressure function $p=\varpi(T,V)$ that gives the equilibrium pressure of a material in state $(T,V)$. Pressure function is further assumed to satisfy the mechanical stability condition expressed by $\tfrac{\partial \varpi}{\partial V} < 0$. One can invert the functional relation $p=\varpi(T,V)$ to solve for volume $V=\varphi(T,p)$ and feed this to Eq. (\ref{doctrine}) to arrive at \begin{eqnarray}\nonumber h&=&\Lambda_{V}(T,V)\left\{ \frac{\partial V}{\partial T}\dot{T}+\frac{\partial V}{\partial p}\dot{p}\right\} + C_{V}(T,V)\dot{T}\\ \nonumber &=& \Lambda_{V}(T,V)\frac{\partial V}{\partial p}\dot{p} + \left\{ \Lambda_{V}(T,V)\frac{\partial V}{\partial T} + C_{V}(T,V)\right\} \dot{T} \\ &=:& \Lambda_{p}(T,V) \dot{p} + C_{p}(T,V) \dot{T} \end{eqnarray} where the specific heat at constant pressure is $C_{p}(T,V):=C_{V}(T,V)+\Lambda_{V}(T,V)\tfrac{\partial V}{\partial T}$ and $\Lambda_{p}$ is analogously defined. Obviously, if one plugs in $V=\varphi(T,p)$ to these relations then everything can be transformed to $Tp$ space.

A quantity that is of much interest to the studies of thermodynamics and fluid mechanics is the ratio of two heat capacities, $\gamma := C_{p}/C_{V}$. Clearly $\gamma = 1 + \tfrac{\Lambda_{V}}{C_{V}} \tfrac{\partial V}{\partial T}$. This can be exploited in Eq. (\ref{adiabat}) to read \begin{equation} \frac{{\rm d}\psi(V)}{{\rm d}V} = -(\gamma (\psi(V),V)-1) \frac{\partial T}{\partial V} (\psi(V),V) . \label{adiabat2} \end{equation} This equation is exact in the sense that it only postulates the validity of The Doctrine of Latent and Specific Heats and the existence of a pressure function of a mechanically stable material. Now we are in a position to state a bold and simplifying assumption which requires $\gamma$ to be a constant for an ideal gas. This adscititious hypothesis must be supported by experimental evidence and measurements suggest that for monatomic gases $\gamma \approx 5/3$. Since we are also working with an ideal gas fluid $\varpi(T,V) := K T/V$ where $K$ is a positive constant. Then $\tfrac{\partial T}{\partial V} = p/K = T/V$. Substitute this to Eq. (\ref{adiabat2}) to read \begin{equation} \frac{{\rm d}\psi(V)}{{\rm d}V} = -(\gamma -1) \frac{\psi(V)}{V} . \label{adiabatideal} \end{equation} Solution of (\ref{adiabatideal}) is an elementary exercise of calculus and left to the reader. It is \begin{equation} T = \psi(V) = \frac{T_{\circ}V^{\gamma-1}_{\circ}}{V^{\gamma-1}} \end{equation} By substituting $T = pV/K$ for an ideal gas one can represent the adiabat in $pV$ space more succinctly as \begin{equation} pV^{\gamma} = p_{\circ}V_{\circ}^{\gamma}. \label{adiabatcurve} \end{equation}

Rüchhardt's method

Rüchhardt's method can be described by the figure on the right which is taken from wikipedia. There is a spherical ball of mass $m$ and radius $r$ that just fits the glass pipe which is immersed into the vessel that encloses an ideal gas inside. Rüchhardt's assumptions are as follows: (1) The ball is frictionless and can move up and down in the pipe with no dissipation of its mechanical energy. (2) During the motion of the ball no gas escapes. (3) There is not enough time for the system to thermally equilibrate with its surroundings, ergo the entire process may be regarded as adiabatic if the duration of the experiment is not too long.

As a preliminary lemma we start by calculating the total (net) force acting on a half sphere which suffers a constant gas pressure $p$. It is agreed that the north pole of the sphere is aligned in $\mathbf{k}$ direction. We have ${\rm d}{\mathbf F} = -p {\mathbf {\hat{n}}}{\rm d}S$ where ${\rm d}S$ is the area element and $\mathbf{\hat{n}}$ is the unit normal to that area element. For a spherical surface of radius $r$ they are ${\rm d}S := r^{2}\sin \theta {\rm d}\theta {\rm d}\phi$ and $\mathbf{\hat{n}} := \sin \theta \cos \phi \mathbf{i} + \sin \theta \sin \phi \mathbf{j} + \cos \theta \mathbf{k}$. Then the net force acting on the half sphere can be calculated by an elementary exercise of calculus to read \begin{equation*} \mathbf{F} = -pr^{2}\int\limits_{0}^{2\pi}{\rm d}\phi \int\limits_{0}^{\pi/2} {\rm d}\theta \sin \theta \mathbf{\hat{n}} = -p\pi r^{2} \mathbf{k} . \end{equation*} Note that due to the symmetry of the surface, only the vertical component of the net force is non vanishing.

At equilibrium let the internal pressure of the gas be $p_{\rm eq}$ and the atmospheric pressure which is constant throughout the experiment be $p_{\circ}$. Then the condition of mechanical equilibrium reads \begin{equation*} p_{\circ} \pi r^{2} + mg = p_{\rm eq} \pi r^{2} \end{equation*} which says that the pressure inside the vessel is slightly larger than the atmospheric pressure. Let the volume of the gas at equilibrium be $V_{\rm eq}$. If one moves the ball from equilibrium by a distance $z$, then the volume of the gas becomes $V = V_{\circ} + z \pi r^{2}$. Consequently, by the assumption of adiabaticity the gas pressure must change into \begin{equation*} p = p_{\rm eq}\frac{V^{\gamma}_{\rm eq}}{(V_{\rm eq} + z \pi r^{2})^{\gamma}} \approx p_{\rm eq} \left( 1 - \frac{z\pi r^{2}}{V_{\rm eq}}\right)^{\gamma} \approx p_{\rm eq} \left( 1-\gamma\frac{z\pi r^{2}}{V_{\rm eq}}\right) \end{equation*} where the assumption of small deviations $z\pi r^{2}/V_{\rm eq} \ll 1$ and application of well known Neuman series and Taylor approximations are used to linearize the pressure of the gas.

Upon the displacement of the ball the net force acting on it can be linearized accordingly. \begin{equation*} \mathbf{F}_{\rm net} = (p\pi r^{2} - mg - p_{\circ}\pi r^{2}) \mathbf{k} \approx p_{\rm eq} \left( 1-\gamma\frac{z\pi r^{2}}{V_{\rm eq}}\right)\pi r^{2}\mathbf{k} - p_{\rm eq} \pi r^{2}\mathbf{k} = -\frac{\gamma p_{\rm eq}\pi^{2} r^{4}}{V_{\rm eq}}z \mathbf{k} \end{equation*} It is time to invoke the second law of motion in the form $\mathbf{a} = \tfrac{1}{m} \mathbf{F}_{\rm net}$ where $\mathbf{a} := \ddot{z} \mathbf{k}$. Then the equation of motion for $z$ is \begin{equation} \ddot{z}(t) \approx -\frac{\gamma p_{\rm eq}\pi^{2} r^{4}}{mV_{\rm eq}}z(t) . \label{sho} \end{equation} Our efforts to linearize the gas pressure pays off and we recognize the equation of motion for a simple harmonic oscillator of natural frequency $\omega$ in Eq. (\ref{sho}) where \begin{equation} \omega :=\sqrt{\frac{\gamma p_{\rm eq} A^{2}}{mV_{\rm eq}}} \label{freq} \end{equation} and the cross section of the glass pipe is substituted as $A := \pi r^{2}$. Let the period of the oscillations be $\tau$. Then the relation $\omega = 2 \pi / \tau$ can be solved for $\gamma$ to read as follows. \begin{equation} \boxed{ \gamma = \frac{4 \pi^{2}mV_{\rm eq}}{p_{\rm eq}A^{2}\tau^{2}} } \label{ruchhardt} \end{equation} This is the Rüchhardt formula to measure $\gamma$ via the measurement of the period of oscillations of the spherical ball $m$.

In the desing of the experiments the laboratory worker needs a $\tau$ that is not too long for otherwise the neglected friction and leakage can become pertinent. Nor a very short $\tau$ makes the detection of a periodic motion feasable with human eye. For the optimal design we should chose the density $d := m/(\tfrac{4}{3}\pi r^{3})$ and the radius of the ball such that $\tau \sim 1$ s or so. Rearrange the Rüchhardt formula to express the period of the motion as \begin{equation} \tau = \sqrt{\frac{16 \pi V_{\rm eq}}{3 p_{\rm eq} \gamma} \frac{d}{r}} . \end{equation} Keep in mind that $\gamma \sim 1$ for almost every material and the apparatus dependence of $\tau$ is given as $(d/r)^{1/2}$. It is up to the designer of the experiment to play with $r,d$, and $V_{\rm eq}$ to bring $\tau$ into a sensible interval. (We do not have much freedom on $p_{\rm eq}$ as it must be very close to the atmospheric pressure.)

Further reading

Present post is obviously influenced by the works of the philomath Clifford A. Truesdell. Owen's textbook is a good source for a beginner of rational thermodynamics. Finally Chapter 5 of Zemansky is where the Rüchhardt method is described in a less quantitative way.

1. Clifford A. Truesdell and Subramanyam Bharatha, The Concepts and Logic of Classical Thermodynamics as a Theory of Heat Engines Rigorously Constructed upon the Foundation Laid by S. Carnot and F. Reech, Springer-Verlag, New York, 1977.
2. Dawid R. Owen, A First Course in the Mathematical Foundations of Thermodynamics, Springer-Verlag, New York, 1984.
3. Mark W. Zemansky and Richard H. Dittman, Heat and Thermodynamics, McGraw-Hill, New York, 1996. 7th edition.