10 Kök Bulma ve Optimizasyon
Bir gezegenin yörüngedeki konumunu veren Kepler denklemi \(E - e \sin E = M\), ya da \(\cos x = x\) ve \(\tan x = x\) gibi denklemlerin çözümü kapalı bir formülle yazılamaz. Beşinci ve daha yüksek dereceden polinomlar için de genel bir kök formülü yoktur. Bu tür denklemleri sayısal olarak çözeriz: bir başlangıç tahmininden yola çıkar, onu adım adım iyileştirir ve istediğimiz kadar basamağa ulaşınca dururuz.
Optimizasyon, yani bir fonksiyonun en küçük ya da en büyük değerini aramak, kök bulmanın yakın akrabasıdır. Türevlenebilir bir fonksiyonun bir iç noktadaki yerel ekstremumunda türev sıfırdır; dolayısıyla minimum aramak çoğu zaman türevin kökünü aramak demektir. Bir kutunun hacmini en büyük yapan ölçüler, bir noktanın bir eğriye en yakın noktası, bir modelin verilere en iyi uyan parametreleri hep bu türden sorulardır.
Bu bölümde önce ikiye bölme ve Newton yöntemlerini kendimiz yazıyoruz. Yöntemlerin teorisi Nümerik Analiz notlarındadır; burada onların Python’daki karşılığını kuruyoruz. Ardından aynı işleri scipy.optimize modülüyle yapıyoruz: tek değişkenli denklemler, denklem sistemleri, tek ve çok değişkenli minimizasyon ve kısıtlı problemler.
10.1 Kökleri Ayırmak: İşaret Değişimi
Bir kökü sayısal olarak aramaya başlamadan önce onun yaklaşık nerede olduğunu bilmemiz gerekir.
Tanım 10.1 (Kök Parantezi) \(f\), \([a, b]\) aralığında sürekli bir fonksiyon olsun. \(f(a)\) ile \(f(b)\) zıt işaretliyse, yani \(f(a)\,f(b) < 0\) ise \([a, b]\) aralığına \(f\) için bir kök parantezi (bracket) denir.
Yani parantez, içinde en az bir kök bulunduğunu kesin olarak bildiğimiz bir aralıktır. Bunu Bolzano Teoremi garanti eder (bkz. Analiz 1): sürekli bir fonksiyon negatif bir değerden pozitif bir değere geçerken sıfırdan geçmek zorundadır. Parantezin içinde birden çok kök de olabilir; teorem yalnız en az bir kök bulunduğunu söyler.
Parantez bulmanın en kolay yolu, \(f\)’yi sık bir ızgarada hesaplayıp ardışık iki değerin işaretine bakmaktır. NumPy ile bu iş tek satırdır: y[:-1] * y[1:] dizisi komşu değerlerin çarpımlarını verir ve bir çarpım negatifse o hücrede işaret değişir (Vektörizasyon ve Broadcasting bölümündeki dilimleme ve vektörel işlemler). Örnek olarak \(f(x) = x^3 - 3x + 1\) polinomunu \([-2{,}5;\ 2{,}5]\) aralığında, adımı \(0{,}5\) olan bir ızgarada tarayalım.
import numpy as np
def f(x):
return x**3 - 3*x + 1
x = np.linspace(-2.5, 2.5, 11) # adım 0.5 olan ızgara
y = f(x)
# komşu iki noktada f'nin işareti değişiyorsa arada bir kök vardır
i = np.nonzero(y[:-1] * y[1:] < 0)[0]
brackets = list(zip(x[i], x[i + 1]))
for a, b in brackets:
print(f"[{a}, {b}] aralığında işaret değişiyor")Çıktı:
[-2.0, -1.5] aralığında işaret değişiyor
[0.0, 0.5] aralığında işaret değişiyor
[1.5, 2.0] aralığında işaret değişiyornp.nonzero(...)[0] koşulun doğru olduğu indisleri verir; i indisi, işaretin x[i] ile x[i + 1] arasında değiştiğini söyler. Taramayı bir grafikle de görmek iyi bir alışkanlıktır. Aşağıdaki kod önceki kodun devamıdır ve Matplotlib ile Grafik Çizimi bölümündeki nesne yönelimli arayüzü kullanır; axvspan işaretin değiştiği aralıkları gölgeler.
import matplotlib.pyplot as plt
xx = np.linspace(-2.5, 2.5, 401)
fig, ax = plt.subplots(figsize=(6, 4))
ax.plot(xx, f(xx), label="f(x) = x³ − 3x + 1")
ax.plot(x, y, "o", label="ızgara noktaları")
ax.axhline(0, color="gray", linewidth=0.8)
for a, b in brackets:
ax.axvspan(a, b, alpha=0.2) # işaret değişen aralıklar
ax.set_xlabel("x")
ax.legend()
fig.savefig("isaret_tarama.png", dpi=150)
print("grafik isaret_tarama.png dosyasına kaydedildi")Çıktı:
grafik isaret_tarama.png dosyasına kaydedildiÜç parantez bulduk. Kübik bir polinomun en çok üç kökü olduğundan bütün kökler bu aralıklardadır. Bu polinomun köklerini kapalı biçimde de biliyoruz: \(x = 2\cos\theta\) yazılırsa \(x^3 - 3x = 2(4\cos^3\theta - 3\cos\theta)\) olur ve bu ifade \(2\cos 3\theta\)’ya eşittir. Denklem \(\cos 3\theta = -\tfrac{1}{2}\) hâline gelir ve kökler
\[ 2\cos\frac{2\pi}{9} \approx 1{,}532, \qquad 2\cos\frac{4\pi}{9} \approx 0{,}347, \qquad 2\cos\frac{8\pi}{9} \approx -1{,}879 \]
bulunur. Bu kesin değerler, bölüm boyunca yöntemlerin gerçek hatasını ölçmemizi sağlayacak.
İşaret taraması iki durumda kök kaçırır. Birincisi, iki kök aynı hücreye düşerse: \((x - 0{,}1)(x - 0{,}3)\) fonksiyonunun iki kökü de \([0;\ 0{,}5]\) hücresindedir, ama \(f(0) = 0{,}03\) ve \(f(0{,}5) = 0{,}08\) aynı işaretli olduğundan tarama hiçbir şey bulmaz. İkincisi, çift katlı köklerde: \((x - 1)^2\) fonksiyonu \(x = 1\)’de sıfır olur ama işaret değiştirmez. Izgarayı sıklaştırmak ilk sorunu hafifletir, ikincisini çözmez. Bir de kök tam bir ızgara noktasına denk gelirse çarpım \(0\) olur ve < 0 testi onu atlar. Grafiğe bakmak her durumda iyi bir kontroldür.
10.2 İkiye Bölme Yöntemini Yazmak
Elimizde bir parantez varsa onu sürekli yarıya indirerek köke istediğimiz kadar yaklaşabiliriz. Yöntemin ayrıntılı anlatımı Nümerik Analiz notlarındadır; burada onu bir Python fonksiyonuna çeviriyoruz.
- \(f(a)\,f(b) < 0\) olduğunu denetle; değilse dur, çünkü parantez yoktur.
- Orta nokta \(p = (a + b)/2\)’yi ve \(f(p)\)’yi hesapla.
- \(f(a)\) ile \(f(p)\) zıt işaretliyse kök \([a, p]\) aralığındadır, \(b = p\) al; değilse kök \([p, b]\) aralığındadır, \(a = p\) al.
- Yarı uzunluk \((b - a)/2\) toleransın altına inene kadar 2. ve 3. adımları tekrarla; son orta noktayı döndür.
Bu tarifi \(f(x) = x^3 - 3x + 1\) polinomunun \([1{,}5;\ 2]\) parantezindeki köküne, \(10^{-10}\) toleransla uygulayalım.
import math
def bisection(f, a, b, tol=1e-10, max_iter=100):
"""f'nin [a, b] parantezindeki bir kökünü ikiye bölmeyle bulur.
(p, n) döndürür: p yaklaşık kök, n hesaplanan orta nokta sayısı.
"""
fa = f(a)
if fa * f(b) > 0:
raise ValueError("f(a) ile f(b) aynı işaretli: parantez yok")
for n in range(1, max_iter + 1):
p = (a + b) / 2
fp = f(p)
if fp == 0 or (b - a) / 2 < tol:
return p, n
if fa * fp < 0: # kök [a, p] içinde
b = p
else: # kök [p, b] içinde
a, fa = p, fp
raise RuntimeError("istenen hassasiyete ulaşılamadı")
def f(x):
return x**3 - 3*x + 1
p, n = bisection(f, 1.5, 2.0)
exact = 2 * math.cos(2 * math.pi / 9)
print(f"kök ≈ {p:.12f}, {n} orta nokta")
print(f"gerçek hata: {abs(p - exact):.1e}")
print("öngörülen adım sayısı:", math.ceil(math.log2(0.5 / 1e-10)))Çıktı:
kök ≈ 1.532088886190, 33 orta nokta
gerçek hata: 4.8e-11
öngörülen adım sayısı: 33Fonksiyon parantez koşulunu en başta denetler ve sağlanmıyorsa bir ValueError fırlatır (Sınıflar, Hata Yakalama ve Dosyalar). fa değerini saklayarak her adımda \(f\)’yi yalnız bir kez çağırırız; \(f\)’nin hesabı pahalıysa asıl maliyet çağrı sayısıdır. max_iter sınırı, bir hata yüzünden sonsuz döngüye girmeyi önler.
Kodun son satırı, kaç orta nokta gerektiğini önceden hesaplıyor. Bunu şu sınır sağlar.
Önerme 10.1 (İkiye Bölme Yönteminin Hata Sınırı) \(f \in C[a, b]\) ve \(f(a)\,f(b) < 0\) olsun. İkiye bölme yönteminin \(n\)-inci orta noktası \(p_n\) ise \(f\)’nin \([a, b]\) aralığındaki bir \(p\) kökü için
\[ |p_n - p| \le \frac{b - a}{2^n} \]
olur.
İspatı için bkz. Nümerik Analiz. Hatanın bir \(\varepsilon\) toleransından küçük olması için \(\frac{b - a}{2^n} < \varepsilon\), yani \(n > \log_2 \frac{b - a}{\varepsilon}\) yeterlidir. Örneğimizde \(b - a = 0{,}5\) ve \(\varepsilon = 10^{-10}\) olduğundan \(n > \log_2(5 \cdot 10^{9}) \approx 32{,}2\), yani \(n = 33\) çıkar. Fonksiyonun saydığı orta nokta sayısı da tam budur ve gerçek hata \(4{,}8 \cdot 10^{-11}\) ile sınırın altında kalır.
10.3 Newton Yöntemini Yazmak
İkiye bölme güvenilirdir ama yavaştır: hata sınırı her adımda yalnız yarıya iner. Fonksiyonun türevini de kullanırsak çok daha hızlı ilerleyebiliriz. Kontrol Yapıları ve Fonksiyonlar bölümünde karekökü hesaplarken kullandığımız iterasyon, aslında Newton yönteminin \(x^2 - a = 0\) denklemine uygulanmasıydı. Genel biçimde yöntem
\[ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} \]
kuralıyla çalışır: \(f\) yerine onun \(x_n\) noktasındaki teğetini koyar ve teğetin kökünü bir sonraki yaklaşım olarak alır (bkz. Nümerik Analiz). Yöntemlerin hızını karşılaştırabilmek için önce hız kavramını kesinleştirelim.
Tanım 10.2 (Yakınsama Mertebesi) \((x_n)\) dizisi \(p\)’ye yakınsasın ve \(e_n = |x_n - p|\) olsun. Bir \(\alpha \ge 1\) ve bir \(\lambda > 0\) sabiti için
\[ \lim_{n \to \infty} \frac{e_{n+1}}{e_n^{\alpha}} = \lambda \]
ise dizi \(p\)’ye \(\alpha\) mertebesinden yakınsar denir. \(\alpha = 1\) ise (bu durumda \(\lambda < 1\) istenir) yakınsama doğrusal, \(\alpha = 2\) ise kareseldir.
Yani doğrusal yakınsamada hata her adımda kabaca sabit bir oranla küçülür ve doğru basamak sayısı her adımda sabit bir miktar artar. Karesel yakınsamada ise yeni hata eski hatanın karesiyle orantılıdır; hata \(10^{-4}\) iken bir sonraki adımda \(10^{-8}\) civarına iner, yani doğru basamak sayısı her adımda yaklaşık iki katına çıkar.
Teorem 10.1 (Newton Yönteminin Karesel Yakınsaması) \(f\), \(p\)’yi içeren bir aralıkta iki kez sürekli türevlenebilir, \(f(p) = 0\) ve \(f'(p) \ne 0\) olsun. \(x_0\), \(p\)’ye yeterince yakın seçilirse Newton dizisi \(p\)’ye yakınsar ve hiçbir adımda \(x_n = p\) olmuyorsa
\[ \lim_{n \to \infty} \frac{|x_{n+1} - p|}{|x_n - p|^2} = \left|\frac{f''(p)}{2f'(p)}\right| \]
olur. Dolayısıyla \(f''(p) \ne 0\) ise yakınsama karesel, \(f''(p) = 0\) ise karesel yakınsamadan da hızlıdır.
İspat
Yakınsama Nümerik Analiz notlarında ispatlanmıştır; burada yalnız hızı gösteriyoruz. Taylor Teoremi’ne göre \(p\) ile \(x_n\) arasında bir \(\xi_n\) için
\[ 0 = f(p) = f(x_n) + f'(x_n)(p - x_n) + \frac{f''(\xi_n)}{2}(p - x_n)^2 \]
olur. Bu eşitliği \(f'(x_n)\)’e bölüp \(\frac{f(x_n)}{f'(x_n)} = x_n - x_{n+1}\) eşitliğini kullanırsak
\[ x_{n+1} - p = \frac{f''(\xi_n)}{2f'(x_n)}\,(x_n - p)^2 \]
bulunur. \(x_n \to p\) olduğundan \(\xi_n \to p\)’dir. \(f'\) ile \(f''\) sürekli olduğundan \(f'(x_n) \to f'(p)\) ve \(f''(\xi_n) \to f''(p)\) olur. Her iki tarafın mutlak değerini \(|x_n - p|^2\)’ye bölüp limite geçmek iddiayı verir. \(\blacksquare\)
Teoremi kodla gözleyelim. Aşağıdaki newton_raphson fonksiyonu bütün iterasyonları bir listede saklar; böylece her adımın hatasını gerçek kök \(p = 2\cos\frac{2\pi}{9}\) ile karşılaştırabiliriz.
import math
def newton_raphson(f, df, x0, tol=1e-12, max_iter=50):
"""Newton yöntemi: (yaklaşık kök, bütün iterasyonlar) döndürür."""
xs = [x0]
for _ in range(max_iter):
x = xs[-1]
d = df(x)
if d == 0:
raise ZeroDivisionError(f"x = {x} noktasında türev sıfır")
xs.append(x - f(x) / d)
if abs(xs[-1] - x) < tol:
return xs[-1], xs
raise RuntimeError("Newton yöntemi yakınsamadı")
def f(x):
return x**3 - 3*x + 1
def df(x):
return 3*x**2 - 3
p = 2 * math.cos(2 * math.pi / 9) # gerçek kök
root, xs = newton_raphson(f, df, 2.0)
e = [abs(x - p) for x in xs]
for n, (x, err) in enumerate(zip(xs, e)):
print(f"n={n} x_n={x:.15f} hata={err:.1e}")
for n in range(1, 4):
print(f"e_{n+1} / e_{n}^2 = {e[n+1] / e[n]**2:.4f}")
print("öngörülen oran:", round(p / (p**2 - 1), 4))Çıktı:
n=0 x_n=2.000000000000000 hata=4.7e-01
n=1 x_n=1.666666666666667 hata=1.3e-01
n=2 x_n=1.548611111111111 hata=1.7e-02
n=3 x_n=1.532390161865380 hata=3.0e-04
n=4 x_n=1.532088989397224 hata=1.0e-07
n=5 x_n=1.532088886237968 hata=1.2e-14
n=6 x_n=1.532088886237956 hata=0.0e+00
e_2 / e_1^2 = 0.9123
e_3 / e_2^2 = 1.1036
e_4 / e_3^2 = 1.1365
öngörülen oran: 1.1372Hata sütunu \(10^{-1}\), \(10^{-2}\), \(10^{-4}\), \(10^{-7}\), \(10^{-14}\) mertebelerinden geçiyor: üs her adımda yaklaşık iki katına çıkıyor. Altıncı adımda hata \(0\)’dır: \(x_6\), gerçek kökün float olarak hesaplanan değerine eşittir ve makine duyarlığına ulaşılmıştır. \(e_{n+1}/e_n^2\) oranları da teoremin söylediği sabite yaklaşıyor: \(f''(x) = 6x\) ve \(f'(x) = 3x^2 - 3\) olduğundan
\[ \left|\frac{f''(p)}{2f'(p)}\right| = \frac{6p}{2(3p^2 - 3)} = \frac{p}{p^2 - 1} \approx 1{,}1372 \]
bulunur.
Teorem 10.1, \(x_0\)’ın köke yeterince yakın olmasını ister. Uzak bir başlangıçta yöntem döngüye girebilir ya da ıraksayabilir. Aşağıdaki kodun ilk yarısı \(f(x) = x^3 - 2x + 2\) için \(x_0 = 0\)’dan, ikinci yarısı \(f(x) = \arctan x\) için \(x_0 = 1{,}5\)’ten başlar.
import math
def f(x):
return x**3 - 2*x + 2
def df(x):
return 3*x**2 - 2
x = 0.0
for n in range(1, 5):
x = x - f(x) / df(x)
print(f"x_{n} = {x}")
x = 1.5 # arctan için Newton
for n in range(1, 6):
x = x - math.atan(x) * (1 + x**2)
print(f"x_{n} = {x:.6g}")Çıktı:
x_1 = 1.0
x_2 = 0.0
x_3 = 1.0
x_4 = 0.0
x_1 = -1.69408
x_2 = 2.32113
x_3 = -5.11409
x_4 = 32.2957
x_5 = -1575.32İlk denklemde \(x_1 = 0 - \frac{2}{-2} = 1\) ve \(x_2 = 1 - \frac{1}{1} = 0\) olduğundan iterasyon \(0\) ile \(1\) arasında sonsuza kadar döner. İkincisinde \(\arctan\)’ın türevi \(\frac{1}{1 + x^2}\) uzaklarda çok küçük olduğundan teğetler neredeyse yataydır ve her adım bir öncekinden daha uzağa savrulur. (\(|x_0|\), yaklaşık \(1{,}3917452\) olan bir eşiğin üstündeyse iterasyon ıraksar, altındaysa \(0\)’a yakınsar.)
Çare, Newton’u bir parantezle desteklemektir; aşağıda tanıtacağımız brentq tam bunu yapar.
10.4 SciPy ile Tek Değişkenli Denklemler
Kendi fonksiyonlarımız yöntemleri anlamak için değerlidir; gerçek işte ise scipy.optimize modülünün test edilmiş ve kenar durumları düşünülmüş fonksiyonlarını kullanırız. Tek değişkenli denklemler için üç araç tanıyacak, sonra yöntemlerin hızlarını karşılaştıracağız.
brentq: parantezli ve hızlı
brentq(f, a, b) çağrısı Brent yöntemini uygular. Yöntem kökü hiçbir zaman parantezin dışına bırakmaz, ama mümkün olduğunda yarıya bölmek yerine sekant ve ters karesel interpolasyon adımları atar. Böylece ikiye bölmenin güvenilirliğini sekant yönteminin hızıyla birleştirir. Durma toleransları xtol (varsayılan \(2 \cdot 10^{-12}\)) ve rtol (varsayılan değeri makine epsilonunun dört katı, yaklaşık \(8{,}9 \cdot 10^{-16}\)) parametreleridir. full_output=True verilirse kökün yanında, iterasyon ve fonksiyon çağrısı sayılarını taşıyan bir RootResults nesnesi de döner.
import numpy as np
from scipy.optimize import brentq
def f(x):
return x**3 - 3*x + 1
brackets = [(-2.0, -1.5), (0.0, 0.5), (1.5, 2.0)]
exact = 2 * np.cos(2 * np.array([4, 2, 1]) * np.pi / 9)
for (a, b), p in zip(brackets, exact):
r, info = brentq(f, a, b, full_output=True)
print(f"[{a}, {b}]: kök = {r:.15f}, hata = {abs(r - p):.1e}, "
f"{info.function_calls} f çağrısı")Çıktı:
[-2.0, -1.5]: kök = -1.879385241571423, hata = 3.9e-13, 8 f çağrısı
[0.0, 0.5]: kök = 0.347296355333861, hata = 1.1e-16, 8 f çağrısı
[1.5, 2.0]: kök = 1.532088886237956, hata = 0.0e+00, 8 f çağrısıÜç kökün her biri yalnız 8 fonksiyon çağrısıyla bulundu. Kendi ikiye bölme fonksiyonumuz \(10^{-10}\) hassasiyet için iki uç nokta ve 33 orta nokta, toplam 35 çağrı harcamıştı. İlk kökteki \(3{,}9 \cdot 10^{-13}\) hata xtol toleransının altındadır; daha fazla basamak gerekirse xtol küçültülebilir.
root_scalar: bütün yöntemlere tek kapı
root_scalar fonksiyonu farklı yöntemlere tek bir arayüzden ulaştırır. Parantezli yöntemler (bisect, brentq) bracket ister, açık yöntemler ise bir başlangıç noktası x0 ister. Sekant yöntemi isteğe bağlı ikinci bir nokta x1 de alır; verilmezse SciPy onu x0’ın hemen yanında kendisi seçer. Newton yönteminde türev fprime ile verilir; verilmezse SciPy türevi sonlu farklarla yaklaşık hesaplar. Halley yöntemi ise hem fprime hem de ikinci türev fprime2 ister. Aşağıda beş yöntemi aynı köke uyguluyoruz. root_scalar(f, **kw) yazımı sözlükteki anahtarları isimli argümanlar olarak geçirir; örneğin ilk satır root_scalar(f, method="bisect", bracket=[1.5, 2.0]) çağrısıyla aynıdır.
from scipy.optimize import root_scalar
def f(x):
return x**3 - 3*x + 1
def df(x):
return 3*x**2 - 3
def d2f(x):
return 6*x
# her yöntemin istediği bilgiler bir sözlükte
calls = [
{"method": "bisect", "bracket": [1.5, 2.0]},
{"method": "brentq", "bracket": [1.5, 2.0]},
{"method": "secant", "x0": 2.0, "x1": 1.9},
{"method": "newton", "x0": 2.0, "fprime": df},
{"method": "halley", "x0": 2.0, "fprime": df, "fprime2": d2f},
]
for kw in calls:
sol = root_scalar(f, **kw)
print(f"{kw['method']:7s} kök={sol.root:.15f} "
f"iterasyon={sol.iterations:2d} çağrı={sol.function_calls}")Çıktı:
bisect kök=1.532088886238853 iterasyon=38 çağrı=40
brentq kök=1.532088886237956 iterasyon= 7 çağrı=8
secant kök=1.532088886237956 iterasyon= 7 çağrı=8
newton kök=1.532088886237956 iterasyon= 6 çağrı=12
halley kök=1.532088886237956 iterasyon= 4 çağrı=12SciPy’ın bisect yöntemi varsayılan xtol \(= 2 \cdot 10^{-12}\) toleransına inmek için 38 adım atıyor. Newton 6 iterasyonda bitiyor, ama her iterasyonda hem \(f\)’yi hem \(f'\)’yü çağırdığından 12 çağrı yapıyor. Halley yöntemi ikinci türevi de kullanır ve üçüncü mertebeden yakınsar; 4 iterasyonda, her birinde üç çağrıyla biter. Hangi yöntemin “ucuz” olduğu, türevleri hesaplamanın \(f\)’yi hesaplamaya göre ne kadar pahalı olduğuna bağlıdır.
newton: birçok başlangıç noktası birden
scipy.optimize.newton(f, x0, fprime=None, fprime2=None) fonksiyonu türev verilmezse sekant, fprime verilirse Newton, ikisi de verilirse Halley yöntemini kullanır. Asıl gücü şudur: x0 bir NumPy dizisi olursa her bileşen için ayrı bir Newton iterasyonu, hepsi birlikte ve vektörel olarak yürütülür.
import numpy as np
from scipy.optimize import newton
def f(x):
return x**3 - 3*x + 1
def df(x):
return 3*x**2 - 3
print(newton(f, 2.0, fprime=df)) # tek başlangıç noktası
x0 = np.linspace(-3, 3, 12) # 12 başlangıç noktası birden
roots = newton(f, x0, fprime=df)
for a, r in zip(x0, roots):
print(f"x0 = {a:6.3f} -> {r:9.6f}")Çıktı:
1.532088886237956
x0 = -3.000 -> -1.879385
x0 = -2.455 -> -1.879385
x0 = -1.909 -> -1.879385
x0 = -1.364 -> -1.879385
x0 = -0.818 -> 1.532089
x0 = -0.273 -> 0.347296
x0 = 0.273 -> 0.347296
x0 = 0.818 -> 0.347296
x0 = 1.364 -> 1.532089
x0 = 1.909 -> 1.532089
x0 = 2.455 -> 1.532089
x0 = 3.000 -> 1.532089Başlangıç noktalarının çoğu kendisine en yakın köke gidiyor, ama \(x_0 = -0{,}818\) en yakın kök olan \(-1{,}879\)’a değil, en uzaktaki \(1{,}532\)’ye gidiyor. Sebep, bu noktanın \(f'(x) = 0\) olan \(x = -1\)’e yakın olmasıdır: orada teğet neredeyse yataydır ve \(x\) eksenini çok uzakta keser. Aynı kodu \([-3;\ 3]\) aralığındaki 1201 noktalık bir ızgarada (türevin sıfır olduğu ve Newton’un hiç başlayamadığı \(x = \pm 1\) noktaları çıkarılarak) çalıştırıp her başlangıç noktasını vardığı kökün rengine boyarsak aşağıdaki şekli elde ederiz.
newton fonksiyonunun oradan vardığı kökün rengine boyanmıştır (1201 noktalık ızgara); şeridin altındaki 12 nokta, metindeki kodun başlangıç noktalarıdır. Türevin sıfır olduğu x = ±1 yakınında renkler birbirine karışır.Newton yönteminde başlangıç noktası yalnız yakınsayıp yakınsamamayı değil, hangi köke yakınsandığını da belirler. Belirli bir kökü istiyorsak parantezli bir yöntem kullanmak en güvenlisidir.
Yakınsama hızlarının karşılaştırılması
Tanım 10.2 ile tanımladığımız hızı üç yöntem için yan yana görelim. Aşağıdaki kod ikiye bölmenin, sekantın (\(x_0 = 2\), \(x_1 = 1{,}9\)) ve Newton’un (\(x_0 = 2\)) ilk altı adımının hatalarını, yine \(p = 2\cos\frac{2\pi}{9}\) köküne göre yazdırır. Sekant yöntemi Newton’un türevi yerine son iki noktadan geçen kirişin eğimini kullanır (bkz. Nümerik Analiz). İki başlangıç noktasıyla başladığından sekantın \(n\)-inci adımı \(x_{n+1}\)’i üretir; bu yüzden kod sec[n + 1] değerini yazdırır ve tablonun her satırında üç yöntemin de \(n\) adım sonraki hatası yer alır.
import math
def f(x):
return x**3 - 3*x + 1
def df(x):
return 3*x**2 - 3
p = 2 * math.cos(2 * math.pi / 9) # gerçek kök
bis, a, b = [], 1.5, 2.0 # ikiye bölme: orta noktalar
for _ in range(6):
m = (a + b) / 2
bis.append(m)
if f(a) * f(m) < 0:
b = m
else:
a = m
sec = [2.0, 1.9] # sekant: x0 = 2, x1 = 1.9
for _ in range(6):
u, v = sec[-2], sec[-1]
sec.append(v - f(v) * (v - u) / (f(v) - f(u)))
nwt = [2.0] # Newton: x0 = 2
for _ in range(6):
x = nwt[-1]
nwt.append(x - f(x) / df(x))
print(" n ikiye bölme sekant Newton")
for n in range(1, 7):
print(f"{n:2d} {abs(bis[n - 1] - p):9.1e} {abs(sec[n + 1] - p):9.1e}"
f" {abs(nwt[n] - p):9.1e}")Çıktı:
n ikiye bölme sekant Newton
1 2.2e-01 1.1e-01 1.3e-01
2 9.3e-02 3.2e-02 1.7e-02
3 3.0e-02 3.6e-03 3.0e-04
4 8.4e-04 1.3e-04 1.0e-07
5 1.5e-02 5.2e-07 1.2e-14
6 7.0e-03 7.6e-11 0.0e+00Logaritmik eksende doğrusal yakınsama bir doğru, daha yüksek mertebeden yakınsama ise giderek dikleşen bir eğri verir. İkiye bölmenin hatası tek tek adımlarda büyüyebilir (4. adımda orta nokta tesadüfen köke çok yakın düşmüştür), ama hiçbir zaman Önerme 10.1 sınırını aşmaz ve sınır her adımda yarıya iner. Sekant yönteminin mertebesi altın oran \(\frac{1 + \sqrt{5}}{2} \approx 1{,}618\)’dir; Newton’un mertebesi ise 2’dir.
Örnek 10.1 (Kepler Denklemi) Bir gezegenin yörüngedeki konumu, yörüngenin dış merkezliği \(e\) ve ortalama anomali \(M\) verildiğinde Kepler denkleminin
\[ E - e \sin E = M \]
çözümü olan \(E\) açısından (dış merkez anomaliden) hesaplanır. \(e = 0{,}5\) için \(M = 0, \frac{\pi}{4}, \frac{\pi}{2}, \frac{3\pi}{4}, \pi\) değerlerine karşılık gelen \(E\) değerlerini tek bir newton çağrısıyla bulunuz.
Çözüm
Denklemi kök problemine çevirmek. \(g(E) = E - e\sin E - M\) diyelim. \(g'(E) = 1 - e\cos E\) olduğundan \(g'(E) \ge 1 - e = 0{,}5 > 0\)’dır ve \(g\) kesin artandır; her \(M\) için tek bir kök vardır ve Newton yönteminde türev hiçbir zaman sıfır olmaz.
Başlangıç noktası. \(e \sin E\) terimi en çok \(0{,}5\) olduğundan \(E\), \(M\)’den fazla uzaklaşmaz; doğal başlangıç \(E_0 = M\)’dir.
Kod. kepler fonksiyonu M dizisini kullandığından dizi hâlindeki bir E için bütün hesaplar bileşen bileşen yapılır; newton her bileşene kendi iterasyonunu uygular.
import numpy as np
from scipy.optimize import newton
e = 0.5 # dış merkezlik
M = np.linspace(0, np.pi, 5) # ortalama anomaliler
def kepler(E):
return E - e * np.sin(E) - M
def dkepler(E):
return 1 - e * np.cos(E)
E = newton(kepler, M, fprime=dkepler) # başlangıç: E0 = M
for m, val in zip(M, E):
print(f"M = {m:.6f} -> E = {val:.12f}")
print(f"en büyük artık: {np.max(np.abs(kepler(E))):.1e}")Çıktı:
M = 0.000000 -> E = 0.000000000000
M = 0.785398 -> E = 1.261703055253
M = 1.570796 -> E = 2.020979938090
M = 2.356194 -> E = 2.609753954400
M = 3.141593 -> E = 3.141592653590
en büyük artık: 2.2e-16Kontrol. \(M = 0\) ve \(M = \pi\) için \(\sin E = 0\) olduğundan \(E = M\) kesin çözümdür ve kod da bunları verdi. \(M = \frac{\pi}{2}\) için
\[ 2{,}020980 - 0{,}5 \sin(2{,}020980) \approx 1{,}570796 \approx \frac{\pi}{2} \]
bulunur. Bütün denklemlerin artığı \(2{,}2 \cdot 10^{-16}\) civarındadır. \(\blacksquare\)
10.5 Denklem Sistemleri
Birden çok bilinmeyenli denklem sistemlerinde de aynı fikir işler; türevin yerini Jacobi matrisi alır. \(F: \mathbb{R}^n \to \mathbb{R}^n\) için \(F(\mathbf{x}) = \mathbf{0}\) sistemini çözmek istiyoruz. \(F\)’nin \(\mathbf{x}_k\) noktasındaki doğrusal yaklaşımı \(F(\mathbf{x}_k) + J(\mathbf{x}_k)(\mathbf{x} - \mathbf{x}_k)\)’dir; burada \(J\), \(F\)’nin Jacobi matrisidir (bkz. Analiz 4). Bu yaklaşımı sıfır yapan nokta bir sonraki iterasyondur:
\[ J(\mathbf{x}_k)\,\mathbf{s}_k = -F(\mathbf{x}_k), \qquad \mathbf{x}_{k+1} = \mathbf{x}_k + \mathbf{s}_k. \]
Her adımda ters matris hesaplamak yerine bir lineer sistem çözeriz; NumPy ile Lineer Cebir bölümünde solve’un inv’den neden daha iyi olduğunu görmüştük. Örnek olarak
\[ x^2 + y^2 = 4, \qquad e^x + y = 1 \]
sistemini ele alalım. \(F(x, y) = (x^2 + y^2 - 4,\ e^x + y - 1)\) için Jacobi matrisi
\[ J(x, y) = \begin{pmatrix} 2x & 2y \\ e^x & 1 \end{pmatrix} \]
olur. \((1, -1)\) noktasından başlayalım.
import numpy as np
def F(v):
x, y = v
return np.array([x**2 + y**2 - 4, np.exp(x) + y - 1])
def J(v):
x, y = v
return np.array([[2*x, 2*y],
[np.exp(x), 1.0]])
v = np.array([1.0, -1.0]) # başlangıç noktası
for k in range(1, 7):
s = np.linalg.solve(J(v), -F(v)) # J s = -F sistemini çöz
v = v + s
print(f"{k}: x = {v[0]:.12f}, y = {v[1]:.12f}, "
f"|s| = {np.linalg.norm(s):.1e}")Çıktı:
1: x = 1.075765685480, y = -1.924234314520, |s| = 9.3e-01
2: x = 1.009470792611, y = -1.737844860783, |s| = 2.0e-01
3: x = 1.004188579849, y = -1.729653231335, |s| = 9.7e-03
4: x = 1.004168738694, y = -1.729637287086, |s| = 2.5e-05
5: x = 1.004168738475, y = -1.729637287026, |s| = 2.3e-10
6: x = 1.004168738475, y = -1.729637287026, |s| = 0.0e+00Adım uzunlukları \(|\mathbf{s}_k|\) tek değişkenli durumdaki gibi karesel hızla küçülüyor: \(10^{-2}\), \(10^{-5}\), \(10^{-10}\).
SciPy’da aynı işi fsolve ve root yapar. fsolve(F, x0, fprime=J), Newton yöntemini uzak başlangıçlarda da sağlam kalacak biçimde değiştiren “hybrid” (Powell) yöntemini kullanır; Jacobi matrisi verilmezse onu sonlu farklarla kendisi hesaplar. root(F, x0, jac=J) aynı yöntemi daha düzenli bir arayüzle sunar ve sonucu, başarının ve fonksiyon çağrısı sayısının da içinde olduğu bir OptimizeResult nesnesi olarak döndürür.
import numpy as np
from scipy.optimize import fsolve, root
def F(v):
x, y = v
return np.array([x**2 + y**2 - 4, np.exp(x) + y - 1])
def J(v):
x, y = v
return np.array([[2*x, 2*y],
[np.exp(x), 1.0]])
for start in ([1.0, -1.0], [-1.0, 1.0]):
v = fsolve(F, start, fprime=J)
print(f"fsolve, başlangıç {start}: x = {v[0]:.12f}, "
f"y = {v[1]:.12f}")
sol = root(F, [-1.0, 1.0], jac=J) # varsayılan yöntem: "hybr"
print(sol.success, sol.message)
print("x =", sol.x, " nfev =", sol.nfev)
print("F(x) =", F(sol.x))Çıktı:
fsolve, başlangıç [1.0, -1.0]: x = 1.004168738475, y = -1.729637287026
fsolve, başlangıç [-1.0, 1.0]: x = -1.816264068824, y = 0.837367799891
True The solution converged.
x = [-1.81626407 0.8373678 ] nfev = 10
F(x) = [-2.91411340e-12 -2.17825757e-13]İki farklı başlangıç noktası, sistemin iki farklı çözümünü verdi; şekil, çemberle eğrinin tam iki noktada kesiştiğini gösteriyor. fsolve yakınsamadığında yalnız bir uyarı yazar ve yine bir dizi döndürür. Bu yüzden sonucu \(F\)’ye koyup artığı denetlemek ya da root kullanıp sol.success değerine bakmak iyi bir alışkanlıktır.
10.6 Tek Değişkenli Optimizasyon
Şimdi kök bulmadan minimum bulmaya geçiyoruz. Yerel ve mutlak minimum kavramları Analiz 2 notlarındadır. Fermat Teoremi’ne göre türevlenebilir bir fonksiyonun bir iç noktadaki yerel ekstremumunda türev sıfırdır (bkz. Analiz 2).
minimize_scalar(f) tek değişkenli bir fonksiyonun yerel minimumunu türev kullanmadan arar. Varsayılan Brent yöntemi, altın oran aramasını parabolik interpolasyonla birleştirir. Aramanın nereden başlayacağı bracket ile seçilebilir; method="bounded" ve bounds=(a, b) ise aramayı bir aralığa hapseder. Sonuç bir OptimizeResult nesnesidir: x minimum noktası, fun oradaki değer, nfev fonksiyon çağrısı sayısıdır. SciPy’da yalnız minimize eden fonksiyonlar vardır; \(f\)’nin maksimumu, \(-f\)’nin minimumu olarak bulunur.
Örnek olarak iki yerel minimumu olan \(f(x) = x^4 - 4x^2 + x\) fonksiyonunu üç farklı çağrıyla minimize edelim.
from scipy.optimize import minimize_scalar
def f(x):
return x**4 - 4*x**2 + x
r1 = minimize_scalar(f) # varsayılan: Brent
r2 = minimize_scalar(f, bracket=(-3, -2))
r3 = minimize_scalar(f, bounds=(0, 3), method="bounded")
for name, r in [("varsayılan", r1), ("bracket=(-3, -2)", r2),
("bounds=(0, 3)", r3)]:
print(f"{name:16s} x = {r.x:.6f} f(x) = {r.fun:.6f} "
f"nfev = {r.nfev}")Çıktı:
varsayılan x = 1.346997 f(x) = -2.618556 nfev = 15
bracket=(-3, -2) x = -1.472998 f(x) = -5.444192 nfev = 14
bounds=(0, 3) x = 1.346998 f(x) = -2.618556 nfev = 12bracket=(-3, -2) ile yapılan arama soldaki mutlak minimumu, varsayılan arama ve bounds=(0, 3) ile yapılan arama sağdaki yerel minimumu bulur.Varsayılan arama \((0, 1)\) çiftinden başlayıp yokuş aşağı ilerler ve sağdaki vadiye düşer: \(f(1{,}347) \approx -2{,}619\) yalnız bir yerel minimumdur. bracket=(-3, -2) aramayı soldan başlatır ve asıl, mutlak minimumu \(f(-1{,}473) \approx -5{,}444\)’ü bulur. bounds=(0, 3) ise aramayı \([0, 3]\) aralığıyla sınırladığı için yine sağdaki minimumu verir. Yöntem, başladığı vadinin dibini bulur; vadiyi seçmek bize kalır.
Kök bulma ile optimizasyon arasındaki bağı kodda da görelim: \(f'(x) = 4x^3 - 8x + 1\) türevinin köklerini işaret taraması ve brentq ile bulup ikinci türev \(f''(x) = 12x^2 - 8\) ile sınıflandıralım.
import numpy as np
from scipy.optimize import brentq
def f(x):
return x**4 - 4*x**2 + x
def df(x):
return 4*x**3 - 8*x + 1
def d2f(x):
return 12*x**2 - 8
x = np.linspace(-3, 3, 61)
d = df(x)
for i in np.nonzero(d[:-1] * d[1:] < 0)[0]:
c = brentq(df, x[i], x[i + 1])
kind = "minimum" if d2f(c) > 0 else "maksimum"
print(f"x = {c:9.6f} f(x) = {f(c):9.6f} yerel {kind}")Çıktı:
x = -1.472998 f(x) = -5.444192 yerel minimum
x = 0.126000 f(x) = 0.062748 yerel maksimum
x = 1.346997 f(x) = -2.618556 yerel minimumTürevin üç kökü vardır: iki yerel minimum ve aralarında bir yerel maksimum. minimize_scalar sonuçları bu tablodaki minimumlarla beş ondalık basamağa kadar örtüşüyor.
10.7 Çok Değişkenli Optimizasyon
Çok değişkenli bir \(f: \mathbb{R}^n \to \mathbb{R}\) fonksiyonunu minimize etmek için minimize fonksiyonu kullanılır. Bir yerel minimumda gradyan sıfırdır (bkz. Analiz 4); Hessian matrisi pozitif tanımlı olan bir kritik nokta ise kesin yerel minimumdur (bkz. Analiz 4).
Yöntemleri sınamak için klasik bir örnek Rosenbrock fonksiyonudur:
\[ f(x, y) = (1 - x)^2 + 100\,(y - x^2)^2. \]
İki karenin toplamı olduğundan \(f \ge 0\)’dır ve minimum değer \(f(1, 1) = 0\)’dır. Sayısal olarak zor olmasının sebebi, minimumun \(y = x^2\) parabolü boyunca uzanan dar ve kıvrık bir vadinin dibinde olmasıdır: yöntemler vadiye çabuk iner, ama vadinin içinde minimuma doğru yavaş ilerler. Alışılmış başlangıç noktası \((-1{,}2;\ 1)\)’dir.
İlk yöntem Nelder–Mead’dir. Türev kullanmaz: düzlemde bir üçgeni (genel olarak \(n + 1\) köşeli bir simpleksi) yansıtma, genişletme ve büzme hamleleriyle fonksiyonun azaldığı yöne kaydırır.
import numpy as np
from scipy.optimize import minimize
def f(v):
x, y = v
return (1 - x)**2 + 100 * (y - x**2)**2
x0 = np.array([-1.2, 1.0]) # alışılmış başlangıç noktası
res = minimize(f, x0, method="Nelder-Mead")
print(res.message)
print("x =", res.x, " f(x) =", res.fun)
print("iterasyon:", res.nit, " f çağrısı:", res.nfev)Çıktı:
Optimization terminated successfully.
x = [1.00002202 1.00004222] f(x) = 8.177661197416674e-10
iterasyon: 85 f çağrısı: 159Nelder–Mead minimumu buldu, ama 159 fonksiyon çağrısıyla ve ancak \(10^{-4}\) mertebesinde bir doğrulukla; varsayılan durma toleransları xatol ve fatol \(10^{-4}\)’tür. İkinci yöntem BFGS’tir: gradyanı kullanır ve Hessian matrisinin bir yaklaşımını adım adım güncelleyen bir yarı-Newton yöntemidir. Gradyan
\[ \nabla f(x, y) = \bigl(-2(1 - x) - 400x(y - x^2),\ \ 200(y - x^2)\bigr) \]
olarak jac parametresiyle verilebilir; verilmezse SciPy onu sonlu farklarla yaklaşık hesaplar. Aşağıdaki kod önceki kodun devamıdır.
def grad(v):
x, y = v
return np.array([-2 * (1 - x) - 400 * x * (y - x**2),
200 * (y - x**2)])
r1 = minimize(f, x0, method="BFGS") # gradyan verilmedi
r2 = minimize(f, x0, method="BFGS", jac=grad) # gradyan verildi
for name, r in [("gradyansız", r1), ("gradyanlı", r2)]:
print(f"{name:10s} x = {r.x} nit = {r.nit} "
f"nfev = {r.nfev} njev = {r.njev}")Çıktı:
gradyansız x = [0.99999328 0.99998655] nit = 31 nfev = 114 njev = 38
gradyanlı x = [0.99999997 0.99999995] nit = 32 nfev = 39 njev = 39Sonuçtaki nit iterasyon sayısı, njev ise gradyanın kaç kez hesaplandığıdır. Gradyan verilmeyince SciPy her gradyanı sonlu farklarla yaklaşık hesaplar ve bunun için iki ek fonksiyon çağrısı yapar. Burada \(f\) 38 noktada hesaplanmış, her noktada gradyan için 2 çağrı daha yapılmıştır: \(38 + 2 \cdot 38 = 114\). Gradyan verilince fonksiyon çağrısı 39’a düşer ve sonuç daha doğrudur. Gradyanı biliyorsak vermek neredeyse her zaman kazançtır; elle türev almak zahmetliyse SymPy ile Sembolik Hesap bölümündeki diff ve lambdify bu işi yapar.
Yöntemlerin izlediği yolu görmek için callback parametresi kullanılır: SciPy her iterasyondan sonra bu fonksiyonu o anki noktayla çağırır. Aşağıdaki kod yine öncekilerin devamıdır.
paths = {}
for method, jac in [("Nelder-Mead", None), ("BFGS", grad)]:
path = [x0]
minimize(f, x0, method=method, jac=jac,
callback=lambda xk: path.append(xk.copy()))
paths[method] = np.array(path)
print(f"{method:11s}: {len(path)} nokta, son nokta {path[-1]}")Çıktı:
Nelder-Mead: 85 nokta, son nokta [1.00002202 1.00004222]
BFGS : 33 nokta, son nokta [0.99999997 0.99999995]Nokta sayılarında küçük bir ayrıntı var. BFGS yolunda 32 iterasyonun noktalarına x0 eklenince 33 nokta olur. Nelder–Mead ise başlangıç simpleksini de bir iterasyon sayar; bu yüzden 85 iterasyonda callback 84 kez çağrılır ve x0 ile birlikte 85 nokta olur.
xk.copy() önemlidir: optimizasyon fonksiyonu aynı diziyi yerinde güncelleyebilir, kopyalamadan eklenen bütün liste elemanları sonunda aynı diziyi gösterebilir.
Son olarak bulunan noktanın gerçekten minimum olduğunu denetleyelim: gradyan sıfıra yakın olmalı, Hessian pozitif tanımlı olmalıdır. Bu kod da öncekilerin devamıdır.
def hessian(v):
x, y = v
return np.array([[2 - 400 * (y - x**2) + 800 * x**2, -400 * x],
[-400 * x, 200.0]])
print("gradyan:", grad(r2.x))
H = hessian(np.array([1.0, 1.0]))
print(H)
print("özdeğerler:", np.linalg.eigvalsh(H))
print("koşul sayısı:", np.linalg.cond(H))Çıktı:
gradyan: [-1.70451579e-06 8.23259305e-07]
[[ 802. -400.]
[-400. 200.]]
özdeğerler: [3.99360767e-01 1.00160064e+03]
koşul sayısı: 2508.009601277225Gradyan \(10^{-6}\) mertebesindedir ve Hessian’ın iki özdeğeri de pozitiftir; \((1, 1)\) kesin yerel minimumdur. Koşul sayısının yaklaşık \(2508\) olması, eş yükselti eğrilerinin bir yönde diğerine göre çok uzamış olduğunu söyler. Dar vadinin sayısal karşılığı budur (NumPy ile Lineer Cebir bölümündeki koşul sayısı).
Bu bölümdeki bütün yöntemler başlangıç noktasının çevresindeki bir yerel minimumu arar. Fonksiyonun birden çok minimumu varsa hangisinin bulunacağı x0’a bağlıdır (alıştırmalardaki Himmelblau fonksiyonu buna bir örnektir). Mutlak minimum için birkaç farklı başlangıçtan çalıştırıp en iyi sonucu seçmek basit ve etkili bir yoldur; SciPy’da ayrıca differential_evolution, basinhopping ve shgo gibi küresel arama yöntemleri de vardır. Konveks fonksiyonlarda ise her yerel minimum mutlak minimumdur.
10.8 Kısıtlı Optimizasyon
Gerçek problemlerde değişkenler çoğu zaman serbest değildir: uzunluklar pozitif olmalı, bütçe aşılmamalı, hacim sabit kalmalıdır. Bu tür problemleri önce genel bir biçimde yazalım.
Tanım 10.3 (Kısıtlı Optimizasyon Problemi) \(f, g_1, \dots, g_m, h_1, \dots, h_k : \mathbb{R}^n \to \mathbb{R}\) fonksiyonları verilsin.
\[ \begin{aligned} \text{en küçük yap:}\quad & f(\mathbf{x}) \\[1mm] \text{kısıtlar:}\quad & g_i(\mathbf{x}) \ge 0, \quad i = 1, \dots, m, \\[1mm] & h_j(\mathbf{x}) = 0, \quad j = 1, \dots, k \end{aligned} \]
biçimindeki probleme kısıtlı optimizasyon problemi denir. \(f\) amaç fonksiyonu, \(g_i(\mathbf{x}) \ge 0\) koşulları eşitsizlik kısıtları, \(h_j(\mathbf{x}) = 0\) koşulları eşitlik kısıtlarıdır. Bütün kısıtları sağlayan noktaların kümesine uygun bölge denir.
Yani amaç fonksiyonunu yalnız uygun bölgedeki noktalar arasında en küçük yapmaya çalışırız. SciPy’ın minimize fonksiyonu kısıtları tam bu biçimde bekler: {"type": "ineq", "fun": g} sözlüğü \(g(\mathbf{x}) \ge 0\), {"type": "eq", "fun": h} sözlüğü \(h(\mathbf{x}) = 0\) anlamına gelir. Değişkenlerin alt ve üst sınırları ayrıca bounds ile, her değişken için bir (alt, üst) çifti olarak verilir; None o yönde sınır olmadığını söyler. Kısıtlı problemler için en çok kullanılan yöntem, her adımda problemi karesel bir alt probleme yaklaştıran SLSQP’dir; trust-constr yöntemi de aynı işi LinearConstraint ve NonlinearConstraint nesneleriyle yapar.
Örnek olarak \(f(x, y) = (x - 2)^2 + (y - 1)^2\) fonksiyonunu \(x + y \le 2\), \(x \ge 0\), \(y \ge 0\) üçgeninde minimize edelim. Geometrik olarak bu, üçgenin \((2, 1)\) noktasına en yakın noktasını aramaktır. \(x + y \le 2\) kısıtını önce SciPy’ın istediği biçime, \(2 - x - y \ge 0\) biçimine çeviriyoruz.
import numpy as np
from scipy.optimize import minimize
def f(v):
x, y = v
return (x - 2)**2 + (y - 1)**2
def g(v): # x + y <= 2 yerine 2 - x - y >= 0
return 2 - v[0] - v[1]
cons = [{"type": "ineq", "fun": g}]
bounds = [(0, None), (0, None)] # x >= 0, y >= 0
res = minimize(f, x0=[0.0, 0.0], method="SLSQP",
bounds=bounds, constraints=cons)
print(res.message)
print(f"x = {res.x[0]:.6f}, y = {res.x[1]:.6f}, f = {res.fun:.6f}")
print(f"g(x, y) = {g(res.x):.1e}")
print("f'nin gradyanı:", 2 * (res.x - [2, 1]))Çıktı:
Optimization terminated successfully
x = 1.500000, y = 0.500000, f = 0.500000
g(x, y) = -3.3e-16
f'nin gradyanı: [-1. -1.]\((2, 1)\) noktası uygun bölgede değildir (\(2 + 1 > 2\)). Bu yüzden minimum üçgenin \(x + y = 2\) kenarı üzerindedir: \((2, 1)\)’in bu doğru üzerindeki dik izdüşümü \((1{,}5;\ 0{,}5)\) noktasıdır ve \(f = 0{,}5\) olur. Kısıtın değeri \(0\)’dır (yazdırılan \(-3{,}3 \cdot 10^{-16}\) yuvarlama hatasıdır), yani kısıt sınırdadır. Optimum noktada \(\nabla f = (-1, -1)\), kısıt doğrusunun normali \((1, 1)\) ile paraleldir: \(\nabla f + \lambda \nabla(x + y - 2) = \mathbf{0}\) eşitliği \(\lambda = 1\) ile sağlanır. Bu, Lagrange çarpanları koşuludur (bkz. Analiz 4).
SciPy’da "ineq" her zaman \(g(\mathbf{x}) \ge 0\) demektir. \(x + y \le 2\) kısıtını çevirmeden "fun": lambda v: v[0] + v[1] - 2 diye yazarsak \(x + y \ge 2\) problemini çözmüş oluruz. Bu yeni problemde \((2, 1)\) uygundur, çözücü kısıtsız minimumu \((2, 1)\)’i bulur ve “başarılı” der. Hata mesajı çıkmaz; yanlış problemin doğru çözümünü alırız.
Örnek 10.2 (En Az Malzemeli Silindir Kutu) Hacmi \(1000\ \text{cm}^3\) olan kapaklı bir silindir kutunun yüzey alanını en küçük yapan \(r\) yarıçapını ve \(h\) yüksekliğini minimize ile bulunuz.
Çözüm
Problemi kurmak. Yüzey alanı iki daire ile yan yüzeyin toplamıdır:
\[ A(r, h) = 2\pi r^2 + 2\pi r h. \]
\(V = 1000\) olmak üzere hacim koşulu \(\pi r^2 h - V = 0\) eşitlik kısıtıdır. Uzunluklar pozitif olmalıdır; bounds ile \(r, h \ge 0{,}1\) istiyoruz.
Gradyan. \(\frac{\partial A}{\partial r} = 4\pi r + 2\pi h\) ve \(\frac{\partial A}{\partial h} = 2\pi r\)’dir. Gradyanı jac ile veriyoruz.
import numpy as np
from scipy.optimize import minimize
V = 1000.0 # hacim, cm³
def area(v):
r, h = v
return 2 * np.pi * r**2 + 2 * np.pi * r * h
def area_grad(v):
r, h = v
return np.array([4 * np.pi * r + 2 * np.pi * h, 2 * np.pi * r])
def volume_gap(v): # π r² h - V = 0 olmalı
r, h = v
return np.pi * r**2 * h - V
res = minimize(area, x0=[5.0, 10.0], jac=area_grad, method="SLSQP",
bounds=[(0.1, None), (0.1, None)],
constraints=[{"type": "eq", "fun": volume_gap}])
r, h = res.x
print(f"r = {r:.4f} cm, h = {h:.4f} cm, alan = {res.fun:.2f} cm²")
print(f"h / r = {h / r:.4f}")
r_exact = (V / (2 * np.pi)) ** (1 / 3)
print(f"kapalı biçim: r = {r_exact:.4f}, h = {2 * r_exact:.4f}")Çıktı:
r = 5.4193 cm, h = 10.8385 cm, alan = 553.58 cm²
h / r = 2.0000
kapalı biçim: r = 5.4193, h = 10.8385Kapalı biçimle karşılaştırma. Kısıttan \(h = \frac{V}{\pi r^2}\) çekilirse \(A(r) = 2\pi r^2 + \frac{2V}{r}\) olur. \(A'(r) = 4\pi r - \frac{2V}{r^2} = 0\) denkleminden \(r^3 = \frac{V}{2\pi}\), yani \(r = \sqrt[3]{V/(2\pi)} \approx 5{,}4193\) bulunur. Bu durumda \(\pi r^3 = \frac{V}{2}\) olduğundan \(h = \frac{V}{\pi r^2} = 2r\)’dir: en ekonomik kutunun yüksekliği çapına eşittir. Sayısal sonuç bununla yazdırılan bütün basamaklarda örtüşüyor.
Bir uyarı. Aynı çağrı jac=area_grad olmadan, yani gradyan sonlu farklarla hesaplanarak yapılırsa SLSQP bu başlangıç noktasından \(r \approx 5{,}5242\), \(h \approx 10{,}4306\) noktasında durur ve yine “başarılı” der. Sayısal bir optimizasyonun sonucunu her zaman bağımsız bir yolla, burada kapalı biçimle, denetlemek gerekir. \(\blacksquare\)
10.9 Alıştırmalar
Aşağıdaki alıştırmalarda bölümün araçlarını yeni denklemlere ve optimizasyon problemlerine uyguluyoruz; her çözümde kod ve çalıştırılmış çıktısı verilmiştir.
Alıştırma 10.1 (İkiye Bölmede Adım Sayısı) \(x^3 - x - 2 = 0\) denkleminin \([1, 2]\) aralığındaki kökünü \(10^{-8}\) hassasiyetle bulmak için ikiye bölme yönteminin kaç orta nokta hesaplaması gerektiğini Önerme 10.1 ile önceden belirleyiniz ve bölümdeki bisection fonksiyonuyla doğrulayınız.
Çözüm
Parantez. \(f(x) = x^3 - x - 2\) için \(f(1) = -2 < 0\) ve \(f(2) = 4 > 0\) olduğundan \([1, 2]\) bir kök parantezidir.
Öngörü. \(\frac{2 - 1}{2^n} < 10^{-8}\), yani \(2^n > 10^8\) olmalıdır. \(\log_2 10^8 \approx 26{,}58\) olduğundan \(n = 27\) orta nokta yeterlidir.
Doğrulama.
import math
def bisection(f, a, b, tol=1e-10, max_iter=100):
"""f'nin [a, b] parantezindeki bir kökünü ikiye bölmeyle bulur."""
fa = f(a)
if fa * f(b) > 0:
raise ValueError("f(a) ile f(b) aynı işaretli: parantez yok")
for n in range(1, max_iter + 1):
p = (a + b) / 2
fp = f(p)
if fp == 0 or (b - a) / 2 < tol:
return p, n
if fa * fp < 0:
b = p
else:
a, fa = p, fp
raise RuntimeError("istenen hassasiyete ulaşılamadı")
def f(x):
return x**3 - x - 2
print("öngörülen n:", math.ceil(math.log2((2 - 1) / 1e-8)))
p, n = bisection(f, 1.0, 2.0, tol=1e-8)
print(f"p = {p:.10f}, n = {n}, f(p) = {f(p):.1e}")Çıktı:
öngörülen n: 27
p = 1.5213797018, n = 27, f(p) = -3.0e-08Fonksiyon da 27 orta noktada durdu. Denklemin Cardano formülünden gelen gerçek kökü
\[ \sqrt[3]{1 + \sqrt{26/27}} + \sqrt[3]{1 - \sqrt{26/27}} \approx 1{,}5213797068 \]
sayısıdır; bulunan değerin hatası yaklaşık \(5 \cdot 10^{-9}\) olup \(10^{-8}\)’in altındadır. \(\blacksquare\)
Alıştırma 10.2 (Tanjant Denkleminin İlk Pozitif Kökü) \(\tan x = x\) denkleminin en küçük pozitif kökünü brentq ile bulunuz.
Çözüm
Kökün yeri. \(g(x) = \tan x - x\) diyelim. \(\left(0, \frac{\pi}{2}\right)\) aralığında \(\tan x > x\) olduğundan kök yoktur. \(\left(\frac{\pi}{2}, \pi\right)\) aralığında \(\tan x < 0 < x\) olduğundan yine kök yoktur. \(\left(\pi, \frac{3\pi}{2}\right)\) aralığında \(g(\pi) = -\pi < 0\)’dır ve \(x \to \frac{3\pi}{2}^-\) iken \(g(x) \to +\infty\) olur; \(g\) bu aralıkta sürekli ve artan olduğundan tam bir kök vardır.
Tuzak. \([1, 2]\) aralığında da işaret değişir: \(g(1) \approx 0{,}557\) ve \(g(2) \approx -4{,}185\). Ama bu aralıkta \(\tan x\)’in \(\frac{\pi}{2}\) kutbu vardır; \(g\) sürekli olmadığından Tanım 10.1 koşulu sağlanmaz.
Kod. Önce yanlış parantezin ne verdiğine, sonra doğru paranteze bakalım. Sağ uç, kutbun \(10^{-6}\) solunda seçildi.
import numpy as np
from scipy.optimize import brentq
def g(x):
return np.tan(x) - x
# yanlış parantez: [1, 2] içinde tan x'in kutbu π/2 var
bad = brentq(g, 1.0, 2.0)
print(f"[1, 2]: x = {bad:.10f}, g(x) = {g(bad):.2e}")
# doğru parantez: (π, 3π/2), sağ uç kutbun hemen solunda
a, b = np.pi, 1.5 * np.pi - 1e-6
print(f"g(a) = {g(a):.4f}, g(b) = {g(b):.4e}")
r = brentq(g, a, b)
print(f"x = {r:.12f}, g(x) = {g(r):.1e}")Çıktı:
[1, 2]: x = 1.5707963268, g(x) = 1.44e+12
g(a) = -3.1416, g(b) = 1.0000e+06
x = 4.493409457909, g(x) = 8.9e-16brentq yanlış parantezde \(\frac{\pi}{2} \approx 1{,}5707963268\) noktasını döndürdü; orada \(g\)’nin değeri \(10^{12}\) mertebesindedir, yani bu bir kök değil kutuptur. İşaret değişimi ancak sürekli bir fonksiyon için kök garantisi verir; sonucu \(g\)’ye koyup denetlemek bu tür hataları yakalar. Doğru parantezde en küçük pozitif kök \(x \approx 4{,}493409457909\) bulunur. \(\blacksquare\)
Alıştırma 10.3 (Kosinüs Denklemi için Üç Yöntem) \(\cos x = x\) denkleminin kökünü root_scalar ile brentq (\([0, 1]\) parantezi), sekant (\(x_0 = 0\), \(x_1 = 1\)) ve Newton (\(x_0 = 1\)) yöntemleriyle bulunuz ve iterasyon ile fonksiyon çağrısı sayılarını karşılaştırınız.
Çözüm
Hazırlık. \(f(x) = \cos x - x\) için \(f(0) = 1 > 0\) ve \(f(1) = \cos 1 - 1 \approx -0{,}460 < 0\) olduğundan \([0, 1]\) bir parantezdir. Türev \(f'(x) = -\sin x - 1\)’dir ve \([0, 1]\) aralığında \(-1{,}85\) ile \(-1\) arasında kalır, hiç sıfır olmaz.
Kod.
import numpy as np
from scipy.optimize import root_scalar
def f(x):
return np.cos(x) - x
def df(x):
return -np.sin(x) - 1
for kw in [{"method": "brentq", "bracket": [0.0, 1.0]},
{"method": "secant", "x0": 0.0, "x1": 1.0},
{"method": "newton", "x0": 1.0, "fprime": df}]:
sol = root_scalar(f, **kw)
print(f"{kw['method']:7s} x = {sol.root:.15f} "
f"iterasyon = {sol.iterations} çağrı = {sol.function_calls}")Çıktı:
brentq x = 0.739085133215156 iterasyon = 7 çağrı = 8
secant x = 0.739085133215161 iterasyon = 6 çağrı = 7
newton x = 0.739085133215161 iterasyon = 4 çağrı = 8Karşılaştırma. Üç yöntem de kök \(x \approx 0{,}739085133215\) için anlaşıyor; brentq sonucunun son basamağındaki fark xtol toleransının içindedir. Newton en az iterasyonu (4) yapıyor, ama her iterasyonda \(f\) ve \(f'\) çağrıldığından toplam 8 çağrı harcıyor. Sekant 6 iterasyonda 7 çağrıyla bitiyor; türevi hesaplamak pahalıysa sekant en ekonomik seçenektir. \(\blacksquare\)
Alıştırma 10.4 (Sinüs ile Doğrunun Bütün Kesişimleri) \(\sin x = \dfrac{x}{10}\) denkleminin bütün reel köklerini bulunuz.
Çözüm
Kökler nerede olabilir? \(|\sin x| \le 1\) olduğundan bir kökte \(\left|\frac{x}{10}\right| \le 1\), yani \(x \in [-10, 10]\) olmalıdır. Bu sonlu aralığı taramak yeterlidir.
Izgara seçimi. \(x = 0\) bir köktür. 201 noktalı bir ızgara (\(-10, -9{,}9, \dots, 10\)) \(0\)’ı tam içerir; o noktada \(f = 0\) olur, çarpım testi < 0 sağlanmaz ve kök kaçar. Bu yüzden 0’ı içermeyen 200 noktalık bir ızgara alıyoruz. Adım yaklaşık \(0{,}1\) olduğundan iki kökün aynı hücreye düşmesi de beklenmez.
import numpy as np
from scipy.optimize import brentq
def f(x):
return np.sin(x) - x / 10
x = np.linspace(-10, 10, 200) # adım 20/199 ≈ 0.1
y = f(x)
i = np.nonzero(y[:-1] * y[1:] < 0)[0]
roots = [brentq(f, x[k], x[k + 1]) for k in i]
print(len(roots), "kök:")
for r in roots:
print(f" {r: .10f} f = {f(r): .1e}")Çıktı:
7 kök:
-8.4232039324 f = -3.2e-15
-7.0681743581 f = -2.2e-16
-2.8523418945 f = 9.2e-14
0.0000000000 f = 1.2e-17
2.8523418945 f = -9.2e-14
7.0681743581 f = 2.2e-16
8.4232039324 f = 3.2e-15Sonuç. Denklemin 7 kökü vardır: \(0\), \(\pm 2{,}8523418945\), \(\pm 7{,}0681743581\) ve \(\pm 8{,}4232039324\). \(f\) tek fonksiyon olduğundan kökler \(0\)’a göre simetriktir. \(\blacksquare\)
Alıştırma 10.5 (Döngüye Giren Newton İterasyonunu Kurtarmak) \(x^3 - 2x + 2 = 0\) denkleminin reel kökünü, SciPy’ın newton fonksiyonunun \(x_0 = 0\) başlangıcında neden başarısız olduğunu da göstererek güvenilir bir yöntemle bulunuz.
Çözüm
Kaç reel kök var? \(f(x) = x^3 - 2x + 2\) için \(f'(x) = 3x^2 - 2\) türevi \(x = \pm\sqrt{2/3} \approx \pm 0{,}816\) noktalarında sıfırdır. Yerel minimum değeri \(f\bigl(\sqrt{2/3}\bigr) \approx 0{,}911 > 0\) olduğundan \(f\), \(x > -0{,}816\) için hiç sıfır olmaz; tek reel kök sol taraftadır.
Newton neden başarısız? “Newton her zaman yakınsamaz” uyarısında gördüğümüz gibi, \(x_0 = 0\)’dan başlayan iterasyon \(0, 1, 0, 1, \dots\) döngüsüne girer. SciPy 50 iterasyon sonra bir RuntimeError fırlatır; bunu try/except ile yakalıyoruz.
Parantezle çözüm. \(f(-2) = -2 < 0\) ve \(f(-1) = 3 > 0\) olduğundan \([-2, -1]\) bir parantezdir.
from scipy.optimize import brentq, newton
def f(x):
return x**3 - 2*x + 2
def df(x):
return 3*x**2 - 2
try:
newton(f, 0.0, fprime=df)
except RuntimeError as err:
print("Newton, x0 = 0:", err)
print("f(-2) =", f(-2.0), " f(-1) =", f(-1.0))
r = brentq(f, -2.0, -1.0)
print(f"brentq: {r:.15f}")
print(f"Newton, x0 = -2: {newton(f, -2.0, fprime=df):.15f}")Çıktı:
Newton, x0 = 0: Failed to converge after 50 iterations, value is 0.0.
f(-2) = -2.0 f(-1) = 3.0
brentq: -1.769292354238632
Newton, x0 = -2: -1.769292354238631Kök \(x \approx -1{,}769292354239\)’dur. Parantezin bir ucundan, \(x_0 = -2\)’den başlatılan Newton da aynı köke yakınsıyor: başlangıç doğru seçilince Newton yine hızlı ve doğrudur. \(\blacksquare\)
Alıştırma 10.6 (Elips ile Parabolün Kesişimi) \(\dfrac{x^2}{4} + y^2 = 1\) elipsi ile \(y = x^2 - 1\) parabolünün bütün kesişim noktalarını fsolve ile bulunuz ve sonucu kesin değerlerle karşılaştırınız.
Çözüm
Kesin çözüm. \(y = x^2 - 1\) elips denklemine konursa \(\frac{x^2}{4} + (x^2 - 1)^2 = 1\) olur. \(u = x^2\) dersek \(u^2 - \frac{7}{4}u = 0\), yani \(u = 0\) ya da \(u = \frac{7}{4}\) bulunur. Kesişim noktaları \((0, -1)\) ve \(\left(\pm\frac{\sqrt{7}}{2}, \frac{3}{4}\right)\)’tür.
Sayısal çözüm. Başlangıç noktalarını kaba bir çizimden okuyoruz: kesişimlerden biri elipsin alt tepesi yakınında, ikisi sağda ve solda üst yarı düzlemdedir. Bu yüzden \((0{,}5;\ -1{,}5)\), \((2, 1)\) ve \((-2, 1)\) noktalarından başlıyoruz.
import numpy as np
from scipy.optimize import fsolve
def F(v):
x, y = v
return np.array([x**2 / 4 + y**2 - 1, y - x**2 + 1])
found = []
for start in ([0.5, -1.5], [2.0, 1.0], [-2.0, 1.0]):
v = fsolve(F, start)
found.append(v)
print(f"başlangıç {start}: x = {v[0]: .12f}, y = {v[1]: .12f}")
print("√7/2 =", np.sqrt(7) / 2)Çıktı:
başlangıç [0.5, -1.5]: x = 0.000000014406, y = -1.000000000000
başlangıç [2.0, 1.0]: x = 1.322875655532, y = 0.750000000000
başlangıç [-2.0, 1.0]: x = -1.322875655532, y = 0.750000000000
√7/2 = 1.3228756555322954Karşılaştırma. Yan noktalar \(\frac{\sqrt{7}}{2} \approx 1{,}322875655532\) ve \(y = 0{,}75\) ile bütün basamaklarda örtüşüyor. \((0, -1)\) noktasında ise \(x\) yalnız \(1{,}4 \cdot 10^{-8}\) doğrulukla bulundu. Sebep, iki eğrinin bu noktada birbirine teğet olmasıdır: elipsin alt tepesi ile parabolün tepesi aynı yatay teğete sahiptir. Jacobi matrisi
\[ J(x, y) = \begin{pmatrix} x/2 & 2y \\ -2x & 1 \end{pmatrix} \]
\((0, -1)\) noktasında tekildir ve Newton türü yöntemler orada yalnız doğrusal hızla yakınsar. Bu, tek değişkenli durumdaki çift katlı kökün (\(u = x^2 = 0\)) sistemlerdeki karşılığıdır. \(\blacksquare\)
Alıştırma 10.7 (En Büyük Hacimli Açık Kutu) \(30 \times 20\) cm’lik bir kartonun dört köşesinden kenarı \(x\) olan kareler kesilip kenarlar yukarı katlanarak üstü açık bir kutu yapılıyor. Kutunun hacmini en büyük yapan \(x\) değerini minimize_scalar ile bulunuz.
Çözüm
Model. Kutunun tabanı \((30 - 2x) \times (20 - 2x)\), yüksekliği \(x\)’tir:
\[ V(x) = x(30 - 2x)(20 - 2x), \qquad 0 < x < 10. \]
Maksimum için minimum. \(V\)’nin maksimumu \(-V\)’nin minimumudur. Aralık bilindiğinden method="bounded" ve bounds=(0, 10) kullanıyoruz.
import numpy as np
from scipy.optimize import minimize_scalar
def volume(x):
return x * (30 - 2*x) * (20 - 2*x)
res = minimize_scalar(lambda x: -volume(x), bounds=(0, 10),
method="bounded")
print(f"x = {res.x:.6f} cm, hacim = {volume(res.x):.4f} cm³")
x_exact = (25 - 5 * np.sqrt(7)) / 3
print(f"kapalı biçim: x = {x_exact:.6f}, hacim = {volume(x_exact):.4f}")Çıktı:
x = 3.923748 cm, hacim = 1056.3059 cm³
kapalı biçim: x = 3.923748, hacim = 1056.3059Kapalı biçim. \(V(x) = 4x^3 - 100x^2 + 600x\) ve \(V'(x) = 12x^2 - 200x + 600\)’dür. \(V'(x) = 0\) denklemi \(3x^2 - 50x + 150 = 0\) olur; kökleri \(x = \frac{25 \pm 5\sqrt{7}}{3}\)’tür. \(\frac{25 + 5\sqrt{7}}{3} \approx 12{,}74\) aralığın dışındadır. \(V(0) = V(10) = 0\) ve aralığın içinde \(V > 0\) olduğundan kalan kritik nokta \(x = \frac{25 - 5\sqrt{7}}{3} \approx 3{,}923748\) cm maksimumdur; hacim yaklaşık \(1056{,}31\ \text{cm}^3\)’tür. Sayısal sonuç bununla örtüşüyor. \(\blacksquare\)
Alıştırma 10.8 (Parabole En Yakın Nokta) \(y = x^2\) parabolünün \((3, 0)\) noktasına en yakın noktasını ve bu en kısa uzaklığı bulunuz.
Çözüm
Amaç fonksiyonu. Parabolün \((x, x^2)\) noktasının \((3, 0)\)’a uzaklığının karesi \(d(x) = (x - 3)^2 + x^4\)’tür. Karekök artan bir fonksiyon olduğundan uzaklığı ve karesini aynı \(x\) en küçük yapar; karekökten kurtulmak hesabı kolaylaştırır.
Kod. minimize_scalar ile minimumu buluyor, sonra aynı noktayı \(d'(x) = 4x^3 + 2(x - 3)\) türevinin kökü olarak brentq ile doğruluyoruz.
import numpy as np
from scipy.optimize import brentq, minimize_scalar
def dist2(x): # (x, x²) ile (3, 0) arası uzaklığın karesi
return (x - 3)**2 + x**4
res = minimize_scalar(dist2)
print(f"x = {res.x:.8f}, nokta = ({res.x:.6f}, {res.x**2:.6f})")
print(f"uzaklık = {np.sqrt(res.fun):.8f}, √5 = {np.sqrt(5):.8f}")
c = brentq(lambda x: 4*x**3 + 2*(x - 3), 0, 2) # dist2'nin türevi
print("türevin kökü:", c)Çıktı:
x = 1.00000000, nokta = (1.000000, 1.000000)
uzaklık = 2.23606798, √5 = 2.23606798
türevin kökü: 1.0Kesin sonuç. \(d'(x) = 0\) denklemi \(2x^3 + x - 3 = 0\), yani \((x - 1)(2x^2 + 2x + 3) = 0\) olur. İkinci çarpanın diskriminantı \(4 - 24 < 0\) olduğundan tek kritik nokta \(x = 1\)’dir. En yakın nokta \((1, 1)\), en kısa uzaklık \(\sqrt{(1 - 3)^2 + 1^2} = \sqrt{5} \approx 2{,}2361\)’dir. \(\blacksquare\)
Alıştırma 10.9 (Karesel Bir Fonksiyonun Minimumu) \(f(x, y) = x^2 + xy + y^2 - 3x\) fonksiyonunun minimumunu gradyanı vererek minimize ile bulunuz ve sonucu \(\nabla f = \mathbf{0}\) lineer sistemini çözerek doğrulayınız.
Çözüm
Gradyan ve Hessian. \(\nabla f(x, y) = (2x + y - 3,\ x + 2y)\)’dir. \(\nabla f = \mathbf{0}\) koşulu
\[ \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} 3 \\ 0 \end{pmatrix} \]
lineer sistemidir. Katsayı matrisi aynı zamanda \(f\)’nin Hessian matrisidir.
Kod.
import numpy as np
from scipy.optimize import minimize
def f(v):
x, y = v
return x**2 + x*y + y**2 - 3*x
def grad(v):
x, y = v
return np.array([2*x + y - 3, x + 2*y])
res = minimize(f, x0=[0.0, 0.0], jac=grad, method="BFGS")
print("minimize:", res.x, " f =", res.fun, " nit =", res.nit)
A = np.array([[2.0, 1.0], [1.0, 2.0]]) # ∇f = 0 <=> A v = b
b = np.array([3.0, 0.0])
print("solve :", np.linalg.solve(A, b))
print("Hessian özdeğerleri:", np.linalg.eigvalsh(A))Çıktı:
minimize: [ 2.00000101 -1.00000124] f = -2.9999999999986944 nit = 5
solve : [ 2. -1.]
Hessian özdeğerleri: [1. 3.]Sonuç. Lineer sistemin çözümü \((2, -1)\)’dir ve \(f(2, -1) = 4 - 2 + 1 - 6 = -3\) olur. Hessian’ın özdeğerleri \(1\) ve \(3\) pozitif olduğundan \((2, -1)\) kesin yerel minimumdur; \(f\) konveks bir karesel fonksiyon olduğundan bu aynı zamanda mutlak minimumdur. BFGS aynı noktayı 5 iterasyonda, varsayılan gradyan toleransının izin verdiği \(10^{-6}\) doğrulukla buldu. \(\blacksquare\)
Alıştırma 10.10 (Himmelblau Fonksiyonunun Minimumları) Himmelblau fonksiyonu adıyla bilinen
\[ f(x, y) = (x^2 + y - 11)^2 + (x + y^2 - 7)^2 \]
fonksiyonunu minimize ile \((1, 1)\), \((-1, 1)\), \((-1, -1)\) ve \((1, -1)\) başlangıç noktalarından minimize ediniz ve bulunan minimumları karşılaştırınız.
Çözüm
Gradyan. \(a = x^2 + y - 11\) ve \(b = x + y^2 - 7\) dersek \(f = a^2 + b^2\) ve zincir kuralıyla
\[ \frac{\partial f}{\partial x} = 4xa + 2b, \qquad \frac{\partial f}{\partial y} = 2a + 4yb \]
bulunur.
Kod.
import numpy as np
from scipy.optimize import minimize
def f(v):
x, y = v
return (x**2 + y - 11)**2 + (x + y**2 - 7)**2
def grad(v):
x, y = v
a, b = x**2 + y - 11, x + y**2 - 7
return np.array([4*x*a + 2*b, 2*a + 4*y*b])
for start in ([1, 1], [-1, 1], [-1, -1], [1, -1]):
res = minimize(f, start, jac=grad, method="BFGS")
print(f"başlangıç {str(start):8s} -> x = {res.x[0]: .6f}, "
f"y = {res.x[1]: .6f}, f = {res.fun:.1e}")Çıktı:
başlangıç [1, 1] -> x = 3.000000, y = 2.000000, f = 2.1e-17
başlangıç [-1, 1] -> x = -2.805118, y = 3.131312, f = 7.8e-14
başlangıç [-1, -1] -> x = -3.779310, y = -3.283186, f = 4.4e-13
başlangıç [1, -1] -> x = 3.584428, y = -1.848127, f = 7.2e-16Karşılaştırma. Dört başlangıç noktası dört farklı minimum verdi. Her birinde \(f \approx 0\)’dır. \(f\) iki karenin toplamı olduğundan \(f \ge 0\)’dır; dolayısıyla dört noktanın dördü de mutlak minimumdur. Örneğin \((3, 2)\) için \(a = 9 + 2 - 11 = 0\) ve \(b = 3 + 4 - 7 = 0\)’dır. minimize hangi minimumun “havzasında” başlarsa onu bulur; başlangıç noktasını değiştirmek farklı minimumları keşfetmenin basit bir yoludur. \(\blacksquare\)
Alıştırma 10.11 (Düzleme En Yakın Nokta) \(x + 2y + 3z = 14\) düzleminin orijine en yakın noktasını kısıtlı optimizasyonla bulunuz.
Çözüm
Problem. Uzaklığın karesini, \(f(\mathbf{v}) = x^2 + y^2 + z^2 = \mathbf{v} \cdot \mathbf{v}\)’yi, \(\mathbf{n} \cdot \mathbf{v} - 14 = 0\) eşitlik kısıtı altında en küçük yapıyoruz; burada \(\mathbf{n} = (1, 2, 3)\) düzlemin normalidir. Gradyan \(\nabla f = 2\mathbf{v}\)’dir.
import numpy as np
from scipy.optimize import minimize
def f(v):
return v @ v # x² + y² + z²
n = np.array([1.0, 2.0, 3.0]) # düzlemin normal vektörü
cons = [{"type": "eq", "fun": lambda v: n @ v - 14}]
res = minimize(f, x0=np.zeros(3), jac=lambda v: 2 * v,
method="SLSQP", constraints=cons)
print("nokta:", res.x)
print(f"uzaklık = {np.sqrt(res.fun):.10f}, √14 = {np.sqrt(14):.10f}")Çıktı:
nokta: [1.00000013 2.00000006 2.99999989]
uzaklık = 3.7416573674, √14 = 3.7416573868Kesin sonuç. Lagrange koşulu \(2\mathbf{v} = \lambda \mathbf{n}\), en yakın noktanın normal doğrultusunda olduğunu söyler: \(\mathbf{v} = t\,\mathbf{n}\). Kısıttan \(t\,(\mathbf{n} \cdot \mathbf{n}) = 14t = 14\), yani \(t = 1\) bulunur. En yakın nokta \((1, 2, 3)\), uzaklık \(\sqrt{14} \approx 3{,}7416573868\)’dir. SLSQP bunu varsayılan toleransının izin verdiği yaklaşık \(10^{-7}\) doğrulukla verdi. \(\blacksquare\)
Bu bölümde gradyanları ya elle yazdık ya da SciPy’ın sonlu farklarla yaklaşık hesaplamasına bıraktık. Sonlu farkların ne kadar doğru olduğu ve integrallerin sayısal olarak nasıl hesaplandığı bir sonraki bölümün, Sayısal Türev ve İntegral bölümünün konusudur.