SDP - Atış Yöntemi

Yazar

Taygun Bulmus

Yayınlanma Tarihi

26 Ağustos 2024

Atış Yöntemi

Aşağıdaki denklemi çözelim.

\[ y'' = f(x,y(x),y'(x)), \qquad y(a)=\alpha, \qquad y(b)=\beta \text{.} \]

Bu denklemi başlangıç değer problemine çevirmeye çalışalım.

\[ y'' = f(x,y(x),y'(x)), \qquad y(a)=\alpha, \qquad y'(a)=u \text{.} \]

Bulmamız gereken değer \(u\) değeridir. Bu değeri tahmin edeceğiz. Buradan diferansiyel denklemi tıpkı başlangıç değer problemi gibi \(a\)’dan \(b\)’ye kadar çözeceğiz. Bulduğumuz çözümün son noktasındaki değerde yani \(y_{cozum}(b)\), \(\beta\)’yu bulana kadar bu işleme devam edeceğiz.

Sistematik bir tahmin yürütme mekanizması kurmamız lazım. Bunun için \(y(b)\), \(u\)’nun bir fonksiyonu olsun tanımlaması yapalım.

\[ y(b)= \theta(u) \text{.} \]

Buradan \(u\) değerini bulmak için yukarıda verilen denklemin kökünü bulmamız lazım.

\[ \theta(u)- \beta \equiv r(u) = 0 \text{.} \]

Buradaki \(r(u)\) fonksiyonuna artık fonksiyon (residual function) denir. Bu fonksiyonun kökü \(u\) değeridir. \(u\) değerini biliyorsak sınır değer probleminden çevirdiğimiz başlangıç değer problemini çözebiliriz.

Artık fonksiyonunun kökü \(u\) değerini verecek. Kök bulmak için Ridder’in yöntemi kullanabiliriz. (Bisection veya Newton-Raphson yöntemleri değil.)

Çözüme Giden Adımları

  1. İkinci dereceden diferansiyel denklemimiz var. \(y''=f(x,y(x),y')(x)\).
  2. Çözüm kümesinin sınır değerlerini biliyoruz. \(y(a)=\alpha\), \(y(b)=\beta\).
  3. Bizim çözebildiğimiz denklem tipinde başlangıç koşulları \(y(a)=\alpha\) ve \(y'(a)=u\) şeklinde olmalı.
  4. \(u\) değerini tahmin edelim.
  5. Tahmin ettiğimiz \(u\) değerini kullanarak \(y''=f(x,y(x),y'(x))\) denklemini \(b\) noktasına kadar çözelim.
  6. Elde ettiğimiz çözümün \(b\)’deki değeri \(y(b)=\beta\) değerini veriyor mu? Eğer vermiyorsa \(u\) değerini değiştirerek tekrar deneyelim.
  7. \(u\)’nun bu değişimi sanki bir fonksiyon gibi olsun. Yukarıda bahsettiğimiz \(\theta(u)\) fonksiyonu budur.
  8. \(\theta(u)\) fonksiyonu \(\beta\)’ya eşit olursa \(u\) değerini bulduk demektir. Yani \(\theta(u)- \beta = 0\) denkleminin kökünü bulmaya çalışıyoruz.
  9. \(y(b)=\beta\) değerine yeteri kadar yaklaştı ise \(u\) değerini yani başlangıç koşulunu yani \(y(x)\) fonksiyonunu bulduk demektir.
  10. Denklem çözüldü.

Alıştırma 1

Aşağıda verilen difarensiyel denklemi atış yöntemi ile çözelim. \[ y''(t) = -g, \qquad y(0)=0, \qquad y(5)=50 \text{.} \]

Çözüm

Çözüm adımlarından dördüncü olan ile başlayalım. \(u\) değerini tahmin edelim. \(u=25\) olsun. Bu değerle denklemi çözelim. Çözeceğimiz denklem aşağıdaki gibi olacak.

\[ y''(t) = -g, \qquad y(0)=0, \qquad y'(5)=25 \text{.} \]

Bu denklemin çözümü

\[ \begin{align*} y'(t)&= -gt + C_{1} \text{,} \qquad y'(5)=25 \text{ ,} \\ &= -gt + 25 \text{ ,} \\ y(x)&= -\frac{gt^2}{2} + 25t + C_{2} \text{,} \qquad y(0)=0 \text{,} \\ &= -\frac{gt^2}{2} + 25t \text{.} \end{align*} \]

Şimdi tahminimiz doğru mu deneyelim. Yani \(y(t=5)\) değeri \(50\) mi olacak?

\[ y(t=5)= -\frac{10\times 25}{2} + 25\times 5 = 0 \ne 50 \text{.} \]

Yanlış! Farklı bir \(u\) değeri almamız gerekecek. Aynı işlemi \(u=40\) ile tekrarlayalım. Sonuç:

\[ y(t)=-\frac{gt^2}{2} + 40t \]

Olacak denersek eğer \(y(5)=75\) olacak ve yanlış olacak.

Teker teker denemek yerine \(u\) değeri sanki bir fonksiyonmuş gibi davranalım.

\[ \theta(u) - 50 = 0\text{.} \]

O halde \(\theta(u)\) fonksiyonu başlangıç değer probleminin çözümüne bağlı. Bunu bir fonksiyon gibi tanımlayıp kökünü bulalım. Bu fonksiyon kodda şöyle tanımlanacak. u değeri girdi olacak. Bu u değerine göre başlangıç değer problemini çözecek. Çözümün son değeri \(y(5)\) değeri olacak. Bu değer ile \(50\) değeri arasındaki farkı döndürecek. Bu farkın kökünü bulmak için scipy.optimize.ridder(fonksiyon, kok_araligi_baslangic, kok_araligi_bitis) fonksiyonunu kullanacağız.

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as scInteg
import scipy.optimize as scOpt
# Diff. Denk:
# y'' = -10 , y(0)=0 , y(5)=50
# Fonksiyon:
def fonk_SDP(t, y):
    # y[0]= y
    # y[1]= v
    return np.array([y[1], -10])
# Başlangıç Koşulları
t0=0.0
tSon=5.0
y0= 0.0  # y(0) = 0
ySon= 50. # y(5) = 50
# v(0) = u (Bilmiyoruz)
# Theta fonksiyonu:
def thetaFonk(u):
    #solve_ivp(fonksiyon, t_span, y0)
    # t_span=[t0, tStop] (y(0) ve y(1) değerini biliyoruz.)
    cozum = scInteg.solve_ivp(fonk_SDP, [t0, tSon], [y0, u])
    
    # İkinci derece diferansiyel denklemin çözümü
    # cozum.y[0] = y(t)
    # cozum.y[1] = v(t)
    yCozum = cozum.y[0]
    # yCozum bir array. Son elemanı yani y(t=5) değeri bizim için önemli.
    # Bu da yCozum[-1] ile alınır.
    # Eğer bulduğumuz yCozum[-1] değeri ySon değerine eşitse u değerini bulduk demektir.
    return yCozum[-1]- ySon
# ---------------------------
# Theta fonksiyonunun içerini şöyle düşünelim.
# Polinom tipi bir g(x)=x^2+4x+4 fonksiyonu olsun. Bunun kökünü şöyle bulurduk.
'''
def g(x):
    return x**2 + 4*x + 4

# Kökünü 0 ile 10 arasında bul.
x0, = scOpt.ridder(g, 0, 10)
'''
# ---------------------------
# theta fonksiyonunun kökünü bulalım.
thetaU_kok_v0= scOpt.ridder(thetaFonk, 0, 50)
# Olması gereken başlangıç değeri:
print(f"v0= {thetaU_kok_v0} m/s")
# Deneyelim. Gerçekten v0 bu değerdeyken y(5)= 50 oluyor mu?
cozum= scInteg.solve_ivp(fonk_SDP, [t0, tSon], [y0, thetaU_kok_v0])
# Son değerleri karşılaştıralım.
print(f"y(t=5) Bulduğumuz Çözüm: {cozum.y[0][-1]:.1f} m")
print(f"y(t=5) Gerçek Değer    : {ySon} m")
v0= 35.000000000001016 m/s
y(t=5) Bulduğumuz Çözüm: 50.0 m
y(t=5) Gerçek Değer    : 50.0 m

Python ile ne yaptık?

  1. Sanki başlangıç değer problemi çözer gibi fonksiyon tanımladık.
  2. Başlangıç koşullarını tanımladık. İkinci başlangıç değeri yani \(y(b)\) değerini tanımlamadık.
  3. Ayrı bir fonksiyon ile başlangıç değer problemini çözdük. Bu fonksiyonun girdi değeri u, dönüş dğeri ise verilen \(y(b)\) değeri ile çözümün son değerinin farkıdır.
  4. Başlangıç değer çözücü denklemi scipy.optimize.ridder(fonksiyon, kok_araligi_baslangic, kok_araligi_bitis) ile çözdük.

Problemler

Problem 1

Aşağıdaki denklemi atış yöntemi ile çözün.

\[ y''(t)=t, \quad y(0)=1, \quad y(3)=9 \]

  1. Elde ettiğiniz çözümü yani \(y(t)\) fonksiyonunu çizdirin.
  2. Bu denklemi analitik olarak çözün ve çözümü de çizdirin.

Problem 2

Aşağıdaki diferansiyel denklemi atış yöntemini kullanarak çözün. Adım aralığını \(h=0.05\) alın. Analitik sonuçla karşılaştırın.

\[ y''(t) = t^{3}+t+5, \quad y(0)=1, \quad y(6)=545.8 \]

Problem 3

  1. Aşağıdaki diferansiyel denklemi atış (shooting) yöntemini kullanarak çözün.
  2. Çözümün grafiğini çizdirin. Sayısal çözüm noktalarını kırmızı nokta olarak belirleyin.

\[ \frac{d^{2}y(x)}{dx^{2}}+ \frac{dy(x)}{dx} = -4 \text{ ,} \]

Sınır koşulları: \(y(0)=0\), \(y(3)=9\).

Bu soruyu çözerken kökü \(-30\) ile \(30\) arasında arayın.

  1. Aşağıda verilen analitik çözümü de aynı grafikte çizdirin.

\[ y(x) = \frac{\left(-4x + 21e^{3 - x} + e^{3}(4x - 21)\right)}{1 - e^3} \text{ .} \]

Kaynaklar

  1. Numerical Methods in Engineering with Python 3, Jaan Kiusalaas, Cambridge University Press, 2013, Syf:293
  2. Python Programming and Numerical Methods: A Guide for Engineers and Scientists, Qingkai Kong & Timmy Siauw & Alexandre Bayen, 2020, chapter23.02-The-Shooting-Method