9 SymPy ile Sembolik Hesap
Şimdiye kadar bilgisayara hep sayılarla hesap yaptırdık. Bu sayılar yaklaşıktır: Python ile İlk Adımlar bölümünde \(0{,}1 + 0{,}2\) toplamının tam olarak \(0{,}3\) çıkmadığını gördük ve NumPy de aynı 64 bitlik sayılarla çalışır. Oysa matematik yaparken çoğu zaman bir sayı değil, bir formül isteriz. \(x^x\) fonksiyonunun türevini, \(\int x e^x\,dx\) integralini ya da \(\sin x / x\) oranının sıfırdaki limitini kâğıt üzerinde harflerle ve kesin olarak hesaplarız.
Bu tür hesapları yapan programlara bilgisayar cebir sistemi (computer algebra system) denir. SymPy, tamamen Python ile yazılmış, açık kaynaklı bir bilgisayar cebir sistemidir. Ayrı bir program değil, import ile yüklenen sıradan bir kütüphanedir. Bu bölümde SymPy ile sembol tanımlamayı, ifadeleri açıp çarpanlarına ayırmayı, limit, türev, integral ve seri hesaplamayı, denklem ve diferansiyel denklem çözmeyi ve matrislerle kesin hesap yapmayı öğreneceğiz. En sonda da sembolik sonuçları lambdify ile NumPy’ye aktarıp sayısal hesapta kullanacağız.
SymPy pip install sympy komutuyla kurulur; Colab gibi çevrim içi Jupyter ortamlarında hazır gelir. Bütün kodlarda kütüphaneyi import sympy as sp kısaltmasıyla yükleyeceğiz. Analiz ve Lineer Cebir derslerinde elle yaptığınız hesapları SymPy ile denetleyebilirsiniz. Yine de SymPy’nin yanıtını doğru yorumlamak, arkasındaki matematiği bilmeyi gerektirir. Bu yüzden bölüm boyunca her hesabı ilgili ders notuna bağlayacağız.
9.1 Sayısal ve Sembolik Hesap
İlk iş, sembolik hesabın sayısal hesaptan nerede ayrıldığını görmek.
Tanım 9.1 (Sembolik Hesap) Sayıların ve değişkenlerin yaklaşık değerleriyle değil, kesin matematiksel ifadeler olarak işlendiği hesaba sembolik hesap denir. Sembolik hesabın sonucu bir ifadedir: \(\sqrt{8}\) yerine \(2\sqrt{2}\), \(\frac{1}{3} + \frac{1}{6}\) yerine \(\frac{1}{2}\), \((\sin x)'\) yerine \(\cos x\).
Yani sembolik hesap, kâğıt kalemle yaptığımız hesabın bilgisayardaki karşılığıdır. Sayısal hesapta \(\sqrt{2}\), bellekte sonlu bir ikili açılımla tutulur (ekranda \(1{,}4142135623730951\) olarak görünür) ve her işlemde küçük bir yuvarlama hatası birikir (bkz. Nümerik Analiz). Sembolik hesapta ise \(\sqrt{2}\), karesi tam olarak \(2\) olan sayı olarak kalır.
Farkı aynı hesabı iki kütüphaneyle yaparak görelim:
import math
import sympy as sp
print(math.sqrt(8))
print(sp.sqrt(8))
print(math.sqrt(2) ** 2)
print(sp.sqrt(2) ** 2)Çıktı:
2.8284271247461903
2*sqrt(2)
2.0000000000000004
2math.sqrt bir float döndürür. sp.sqrt ise \(\sqrt{8} = 2\sqrt{2}\) sadeleştirmesini yapar ve sonucu kesin biçimde saklar; bu yüzden \(\sqrt{2}\)’nin karesi tam olarak \(2\) çıkar. Kesin bir ifadenin ondalık açılımını istediğimiz anda, istediğimiz basamak sayısıyla evalf metodu ya da sp.N fonksiyonu verir. Kesirler için sp.Rational kullanılır:
import sympy as sp
a = sp.Rational(1, 3) + sp.Rational(1, 6)
print(a)
print(sp.pi)
print(sp.pi.evalf(30))
print(sp.N(sp.sqrt(2), 50))
print(float(sp.sqrt(2)))Çıktı:
1/2
pi
3.14159265358979323846264338328
1.4142135623730950488016887242096980785696718753769
1.4142135623730951sp.pi sembolik \(\pi\) sayısıdır; evalf(30) onun 30 anlamlı basamağa doğru yuvarlanmış değerini verir. \(\pi\)’nin açılımı \(\ldots 3832795\ldots\) diye sürdüğünden otuzuncu basamaktaki \(7\) yukarı yuvarlanıp \(8\) olmuştur; öncekilerin hepsi \(\pi\)’nin basamaklarıyla aynıdır. Bu doğruluk tesadüf değildir: SymPy ondalık açılımı istenen duyarlıkta yeniden hesaplar. float(...) ise ifadeyi sıradan bir Python float’una çevirir. Sembolik dünyadan sayısal dünyaya geçmenin en basit yolu budur; daha güçlüsünü bölümün sonunda lambdify ile göreceğiz.
9.2 Semboller ve İfadeler
Sembolik hesabın ilk adımı, harflerle gösterilen değişkenleri SymPy’ye tanıtmaktır.
Tanım 9.2 (Sembol) SymPy’de bir matematiksel değişkeni temsil eden nesneye sembol denir. Semboller sp.symbols fonksiyonuyla oluşturulur: x, y = sp.symbols("x y") satırı adları x ve y olan iki sembol üretir ve onları aynı adlı iki Python değişkenine bağlar.
Yani x artık bir sayı değil, değeri belirsiz bir harftir. Semboller ve sayılar +, -, *, /, ** işlemleriyle ve sp.sin, sp.exp, sp.log gibi SymPy fonksiyonlarıyla birleştirilerek ifadeler (expression) kurulur. Bir ifadedeki sembolün yerine bir sayı ya da başka bir ifade koymak için subs metodu kullanılır:
import sympy as sp
x, y = sp.symbols("x y")
f = x**2 + 2*x*y + sp.sin(x)
print(f)
print(f.subs(x, 1))
print(f.subs({x: 2, y: 3}))
print(f.subs(x, sp.pi))
print(f.subs(y, x**2))Çıktı:
x**2 + 2*x*y + sin(x)
2*y + sin(1) + 1
sin(2) + 16
2*pi*y + pi**2
2*x**3 + x**2 + sin(x)subs sonucu kendiliğinden sadeleştirir: \(x = \pi\) konunca \(\sin \pi = 0\) terimi düşer. \(\sin 1\) ve \(\sin 2\) ise kesin sayılar olarak kalır, ondalığa çevrilmez. Burada sp.sin yerine math.sin yazılamaz, çünkü math.sin yalnız sayı kabul eder ve sembolle çağrılınca hata verir. NumPy’nin np.sin fonksiyonu da sembollerle çalışmaz; sembolik ifadelerde daima sp. önekli fonksiyonları kullanırız.
x + 1/3 yazıldığında Python önce kendi kuralıyla 1/3 bölmesini yapar ve bir float üretir; SymPy yalnız bu yaklaşık sayıyı görür. Kesin kesir için sp.Rational(1, 3) yazılmalıdır. Kuvvet işaretine de dikkat: Python’da üs alma ** ile yapılır. ^ başka bir işlemdir (bitlerde dışlayan veya, XOR) ve 2 ^ 3 ifadesinin değeri 1’dir.
import sympy as sp
x = sp.symbols("x")
print(x + 1/3)
print(x + sp.Rational(1, 3))
print(2 ^ 3, 2 ** 3)Çıktı:
x + 0.333333333333333
x + 1/3
1 8SymPy bir ifadeyi bir metin olarak değil, bir ağaç olarak saklar.
Tanım 9.3 (İfade Ağacı) Bir SymPy ifadesinin ifade ağacı, iç düğümleri işlemler ve fonksiyonlar (Add, Mul, Pow, sin, …), yaprakları semboller ve sayılar olan ağaçtır. Kökteki düğümün türü func niteliğiyle, alt ağaçları args niteliğiyle okunur.
Yani \(x^2 + 2xy + \sin x\) ifadesi SymPy için “üç terimin toplamı”dır: kökte bir Add düğümü, onun altında bir kuvvet (Pow), bir çarpım (Mul) ve bir sin düğümü vardır. Çıkarma ve bölme için ayrı düğüm yoktur: \(a - b\) ifadesi Add(a, Mul(-1, b)), \(a/b\) ifadesi de Mul(a, Pow(b, -1)) olarak saklanır. Her alt ağacın ham yapısını sp.srepr gösterir:
import sympy as sp
x, y = sp.symbols("x y")
f = x**2 + 2*x*y + sp.sin(x)
print(f.func)
print(f.args)
for term in f.args:
print(sp.srepr(term))Çıktı:
<class 'sympy.core.add.Add'>
(x**2, 2*x*y, sin(x))
Pow(Symbol('x'), Integer(2))
Mul(Integer(2), Symbol('x'), Symbol('y'))
sin(Symbol('x'))Çıktıdan ağacın tamamı okunabilir:
Add düğümünün üç alt ağacı vardır: Pow, Mul ve sin. Yapraklar semboller (x, y) ve tam sayılardır (2). Şekil, yukarıdaki sp.srepr çıktılarının çizimidir.subs gibi işlemler bu ağacı dolaşarak çalışır: f.subs(x, 1) her x yaprağını 1 ile değiştirir, sonra düğümleri aşağıdan yukarıya yeniden hesaplar. SymPy’nin birçok davranışı bu yapıdan gelir. Örneğin iki ifadenin “aynı” olup olmadığı sorusu, ileride göreceğimiz gibi, ağaçların aynı olup olmadığı sorusudur.
9.3 Varsayımlar
Bir sembol hakkında bildiklerimizi SymPy’ye söylemezsek, SymPy en genel durumu düşünür.
Tanım 9.4 (Varsayım) Bir sembol oluşturulurken real=True, positive=True, integer=True, nonnegative=True gibi anahtar sözcüklerle verilen bilgilere o sembolün varsayımları (assumptions) denir. Varsayımı verilmemiş bir sembol, keyfi bir karmaşık sayıyı temsil eder.
Yani SymPy bir dönüşümü ancak sembolün alabileceği her değer için doğruysa yapar. \(\sqrt{x^2} = x\) eşitliği \(x = -1\) için yanlıştır. \(\sqrt{x^2} = |x|\) eşitliği de \(x = i\) için yanlıştır: \(\sqrt{i^2} = \sqrt{-1} = i\) iken \(|i| = 1\)’dir. Bu yüzden varsayımsız bir sembolde \(\sqrt{x^2}\) olduğu gibi kalır; sembolün reel ya da pozitif olduğunu söyleyince sadeleşir:
import sympy as sp
z = sp.symbols("z")
x = sp.symbols("x", real=True)
p = sp.symbols("p", positive=True)
n = sp.symbols("n", integer=True)
print(sp.sqrt(z**2), sp.sqrt(x**2), sp.sqrt(p**2))
print(sp.sin(n * sp.pi), sp.cos(n * sp.pi))
print(p.is_positive, x.is_positive, n.is_integer)Çıktı:
sqrt(z**2) Abs(x) p
0 (-1)**n
True None Truen tam sayı olduğundan SymPy \(\sin n\pi = 0\) ve \(\cos n\pi = (-1)^n\) dönüşümlerini de yapar. Son satırdaki is_positive gibi nitelikler üç değerden birini alır: True, False ya da None. None “bilinmiyor” demektir: reel bir \(x\) pozitif de olabilir, olmayabilir de. Bir if koşulunda None yanlış sayılır. Bu yüzden if x.is_positive: bloğunun çalışmaması “x pozitif değil” demek değildir; yalnız “pozitif olduğu bilinmiyor” demektir.
Örnek 9.1 (Logaritmanın Çarpım Kuralı ve Varsayımlar) \(\ln(ab) = \ln a + \ln b\) dönüşümünü sp.expand_log fonksiyonuyla önce varsayımsız, sonra pozitif \(a\), \(b\) sembolleriyle yaptırınız. Sonuçların farkını açıklayınız.
Çözüm
İki durumu aynı kodda deneyelim. İkinci satırdaki sp.symbols çağrısı a ve b adlarını yeni, pozitif sembollere bağlar:
import sympy as sp
a, b = sp.symbols("a b")
print(sp.expand_log(sp.log(a * b)))
a, b = sp.symbols("a b", positive=True)
print(sp.expand_log(sp.log(a * b)))
print(sp.expand_log(sp.log(a**3 / b)))Çıktı:
log(a*b)
log(a) + log(b)
3*log(a) - log(b)Varsayımsız sembollerde ifade değişmez, çünkü kural karmaşık sayılarda genel olarak yanlıştır. Örneğin \(a = b = -1\) için \(\ln(ab) = \ln 1 = 0\)’dır; oysa karmaşık logaritmanın esas değeriyle \(\ln(-1) = i\pi\) olduğundan \(\ln(-1) + \ln(-1) = 2\pi i\) olur. \(a, b > 0\) olduğu bilinince kural her değer için doğrudur ve SymPy onu uygular; \(\ln(a^3/b) = 3\ln a - \ln b\) dönüşümü de aynı nedenle yapılır. Reel pozitif sayılarla çalışılacaksa sembolleri baştan positive=True ile tanımlamak, force=True seçeneğiyle SymPy’yi kuralı uygulamaya zorlamaktan daha güvenlidir.
\(\blacksquare\)
9.4 İfadeleri Dönüştürmek
Aynı matematiksel ifadenin birçok yazılışı vardır ve SymPy’de her yazılışa götüren ayrı bir fonksiyon bulunur. En sık kullanılanlar şunlardır:
sp.expand: çarpımları ve kuvvetleri açar.sp.factor: bir polinomu rasyonel katsayılı indirgenemez çarpanlarına ayırır.sp.collect: bir sembolün kuvvetlerine göre terimleri toplar.sp.cancel: bir rasyonel fonksiyonu ortak çarpanları sadeleşmiş pay/payda biçimine getirir.sp.apart: bir rasyonel fonksiyonu basit kesirlere ayırır.sp.together: kesirleri ortak paydada toplar.sp.expand_trig,sp.trigsimp: trigonometrik ifadeleri açar ya da sadeleştirir.
Önce polinomlarla başlayalım:
import sympy as sp
x, y = sp.symbols("x y")
print(sp.expand((x + 1)**5))
print(sp.factor(x**4 - 1))
print(sp.factor(x**3 - x**2*y - x*y**2 + y**3))
print(sp.factor(x**2 - 2))
print(sp.factor(x**2 - 2, extension=sp.sqrt(2)))
print(sp.collect(sp.expand((x + y + 1)**2), x))Çıktı:
x**5 + 5*x**4 + 10*x**3 + 10*x**2 + 5*x + 1
(x - 1)*(x + 1)*(x**2 + 1)
(x - y)**2*(x + y)
x**2 - 2
(x - sqrt(2))*(x + sqrt(2))
x**2 + x*(2*y + 2) + y**2 + 2*y + 1factor, çarpanlara ayırmayı rasyonel sayılar üzerinde yapar. \(x^2 + 1\) ve \(x^2 - 2\) rasyonel katsayılı çarpanlara ayrılmaz; bu yüzden ilki bir çarpan olarak kalır, ikincisi de değişmeden geri döner. extension=sp.sqrt(2) seçeneği sayı kümesine \(\sqrt{2}\)’yi ekler; o zaman \(x^2 - 2 = (x - \sqrt{2})(x + \sqrt{2})\) ayrışması bulunur.
Rasyonel fonksiyonlarda üç fonksiyon birbirini tamamlar:
import sympy as sp
x = sp.symbols("x")
r = (x**3 - 1) / (x**2 - 1)
print(r)
print(sp.cancel(r))
print(sp.apart((x**2 + 1) / (x**3 - x)))
print(sp.together(1/x + 1/(x + 1)))Çıktı:
(x**3 - 1)/(x**2 - 1)
(x**2 + x + 1)/(x + 1)
1/(x + 1) + 1/(x - 1) - 1/x
(2*x + 1)/(x*(x + 1))cancel, pay ve paydadaki ortak \(x - 1\) çarpanını sadeleştirdi. Bu, \(x = 1\) dışında geçerli bir eşitliktir; özgün kesir \(x = 1\)’de tanımsızdır, sadeleşmiş kesir değildir. apart, Analiz 2’deki basit kesirlere ayırma algoritmasını uygular (bkz. Analiz 2). together ise bunun tersini yapar.
Trigonometrik ifadeler ve genel amaçlı sp.simplify için:
import sympy as sp
x = sp.symbols("x")
print(sp.expand_trig(sp.sin(2*x)))
print(sp.expand_trig(sp.cos(3*x)))
print(sp.trigsimp(sp.sin(x)**4 - sp.cos(x)**4))
print(sp.simplify(sp.sin(x)**2 + sp.cos(x)**2))
print(sp.simplify((x**2 - 1) / (x - 1)))Çıktı:
2*sin(x)*cos(x)
4*cos(x)**3 - 3*cos(x)
-cos(2*x)
1
x + 1simplify, yukarıdaki dönüşümlerin birçoğunu sırayla dener ve bulduğu en kısa ifadeyi döndürür. Hangi dönüşümün işe yarayacağını bilmediğimizde kullanışlıdır. Ama yavaştır ve “en sade biçim” kesin tanımlı bir kavram olmadığından sonucun biçimi önceden kestirilemez. Hedef biçim belliyse (çarpanlara ayrılmış, açılmış, basit kesirlere ayrılmış) ona götüren özel fonksiyonu kullanmak daha iyidir.
İki ifadeyi karşılaştırırken bir inceliğe dikkat etmek gerekir.
Tanım 9.5 (Yapısal Eşitlik) SymPy’de a == b karşılaştırması, iki ifadenin ağaçlarının aynı olup olmadığını sınar. Buna yapısal eşitlik denir. İki ifadenin matematiksel olarak eşit olması, yani değişkenlerin her değeri için aynı sayıyı vermesi ayrı bir sorudur.
Yani \((x+1)^2\) ile \(x^2 + 2x + 1\) matematiksel olarak eşit, yapısal olarak farklıdır: birinin kökünde Pow, ötekinin kökünde Add vardır. Matematiksel eşitliği sınamak için farkı sadeleştirip sıfırla karşılaştırırız ya da equals metodunu kullanırız:
import sympy as sp
x = sp.symbols("x")
a = (x + 1)**2
b = x**2 + 2*x + 1
print(a == b)
print(sp.expand(a) == b)
print(sp.simplify(a - b) == 0)
print(a.equals(b))Çıktı:
False
True
True
True- Eşitliğin iki tarafını ayrı ifadeler olarak yazın:
lhsverhs. - Farkı sadeleştirin:
sp.simplify(lhs - rhs). - Sonuç
0ise özdeşlik doğrudur. Sıfırdan farklı bir ifade çıkarsa farkı birkaç sayıdasubsile hesaplayın: sıfırdan farklı bir değer görülürse iddia yanlıştır. Bu kontrol gereklidir, çünküsimplifyher sıfırı tanıyamayabilir.
Örnek 9.2 (Bir Trigonometrik Özdeşlik) \((\sin x + \cos x)^2 = 1 + \sin 2x\) özdeşliğini SymPy ile doğrulayınız.
Çözüm
Tarifin üç adımını uygulayalım:
import sympy as sp
x = sp.symbols("x")
lhs = (sp.sin(x) + sp.cos(x))**2
rhs = 1 + sp.sin(2*x)
print(sp.simplify(lhs - rhs))
print(lhs == rhs)Çıktı:
0
FalseFark sıfıra sadeleştiğinden özdeşlik her \(x\) için doğrudur. İkinci satır, iki tarafın yapısal olarak farklı olduğunu hatırlatır (Tanım 9.5). Elle de kolayca görülür:
\[(\sin x + \cos x)^2 = \sin^2 x + \cos^2 x + 2\sin x \cos x = 1 + \sin 2x.\]
\(\blacksquare\)
9.5 Limit
Analiz 1’de ε–δ ile tanımladığımız limit (bkz. Analiz 1) SymPy’de tek bir fonksiyonla hesaplanır: sp.limit(f, x, a). Sonsuz, iki küçük o harfiyle sp.oo olarak yazılır.
import sympy as sp
x = sp.symbols("x")
n = sp.symbols("n", positive=True, integer=True)
print(sp.limit(sp.sin(x) / x, x, 0))
print(sp.limit((1 - sp.cos(x)) / x**2, x, 0))
print(sp.limit((1 + 1/n)**n, n, sp.oo))
print(sp.limit(x * sp.log(x), x, 0, dir="+"))
print(sp.limit(x**2 / sp.exp(x), x, sp.oo))Çıktı:
1
1/2
E
0
0Üçüncü satırdaki E, SymPy’nin \(e\) sayısıdır (bkz. Analiz 1). SymPy limitleri, değişkene küçük sayılar koyarak tahmin etmez; ifadenin seri açılımları üzerinde çalışan bir algoritmayla kesin olarak hesaplar. Bu yüzden \(0 \cdot \infty\) ya da \(\infty/\infty\) gibi belirsiz biçimler de sorun çıkarmaz.
Tek yönlü limitler dir seçeneğiyle istenir: dir="+" sağdan, dir="-" soldan, dir="+-" iki yönlü limit demektir. \(\arctan(1/x)\) fonksiyonunun sıfırdaki davranışına bakalım:
import sympy as sp
x = sp.symbols("x")
f = sp.atan(1 / x)
print(sp.limit(f, x, 0, dir="+"))
print(sp.limit(f, x, 0, dir="-"))
try:
print(sp.limit(f, x, 0, dir="+-"))
except ValueError:
print("İki yönlü limit yok: sol ve sağ limit farklı.")Çıktı:
pi/2
-pi/2
İki yönlü limit yok: sol ve sağ limit farklı.\(x\) sağdan sıfıra yaklaşırken \(1/x \to +\infty\) olduğundan \(\arctan(1/x) \to \pi/2\); soldan yaklaşırken \(1/x \to -\infty\) olduğundan \(\arctan(1/x) \to -\pi/2\) olur. Sağ ve sol limit farklı olduğundan iki yönlü limit yoktur (bkz. Analiz 1) ve SymPy ValueError hatası verir. Hatayı try/except ile yakaladık.
dir verilmezse sp.limit sağ limiti hesaplar. Bu yüzden sp.limit(1/x, x, 0) sonucu oo olur, oysa \(1/x\) fonksiyonunun sıfırda iki yönlü limiti yoktur. İki yönlü limit isteniyorsa dir="+-" yazılmalı ya da iki yön ayrı ayrı hesaplanmalıdır:
import sympy as sp
x = sp.symbols("x")
print(sp.limit(1 / x, x, 0))
print(sp.limit(1 / x, x, 0, dir="-"))
print(sp.limit(1 / x**2, x, 0, dir="+-"))Çıktı:
oo
-oo
oo9.6 Türev
Türev sp.diff fonksiyonuyla alınır: sp.diff(f, x) birinci türevi, sp.diff(f, x, n) \(n\)’inci türevi, sp.diff(f, x, y) önce \(x\)’e, sonra \(y\)’ye göre kısmi türevi verir. Türevin tanımı ve kuralları için bkz. Analiz 2.
import sympy as sp
x, y = sp.symbols("x y")
print(sp.diff(x**x, x))
print(sp.diff(sp.sin(x)**2, x, 2))
print(sp.diff(sp.exp(x * y**2), x, y))
d = sp.Derivative(sp.atan(x), x)
print(d, "=", d.doit())Çıktı:
x**x*(log(x) + 1)
2*(-sin(x)**2 + cos(x)**2)
2*y*(x*y**2 + 1)*exp(x*y**2)
Derivative(atan(x), x) = 1/(x**2 + 1)İlk satır, logaritmik türevle bulunan \((x^x)' = x^x(\ln x + 1)\) formülüdür. İkinci satır \((\sin^2 x)'' = 2(\cos^2 x - \sin^2 x) = 2\cos 2x\) sonucunu açılmış biçimde verir. sp.Derivative türevi hesaplamadan, yalnız gösterir; doit metodu hesabı yapar. Aynı ayrım integral, limit ve toplam için de vardır: sp.Integral, sp.Limit, sp.Sum.
Türevi bir noktada hesaplamak için önce türev ifadesini bulur, sonra subs ile noktayı koyarız. \(f(x) = x^3 - 2x\) eğrisinin \(x = 1\) noktasındaki teğetini bulalım. Teğet doğrusu \(y = f(a) + f'(a)(x - a)\) denklemiyle verilir (bkz. Analiz 2):
import sympy as sp
x = sp.symbols("x")
f = x**3 - 2*x
a = 1
slope = sp.diff(f, x).subs(x, a)
tangent = f.subs(x, a) + slope * (x - a)
print("eğim:", slope)
print("teğet: y =", sp.expand(tangent))
print(sp.factor(f - tangent))Çıktı:
eğim: 1
teğet: y = x - 2
(x - 1)**2*(x + 2)Son satır ilginç bir şey söyler: eğri ile teğet arasındaki fark \((x - 1)^2(x + 2)\)’dir. \(x = 1\)’deki çift kök teğetliği gösterir; eğri doğruya bu noktada değer ama onu kesmez. Basit kök \(x = -2\) ise teğetin eğriyi ikinci kez kestiği \((-2, -4)\) noktasıdır.
9.7 İntegral
İntegral için tek fonksiyon vardır: sp.integrate(f, x) belirsiz integrali, sp.integrate(f, (x, a, b)) \([a, b]\) üzerindeki belirli integrali hesaplar. Sınırlar sp.oo olabilir; birden çok sınır üçlüsü verilirse katlı integral içten dışa doğru hesaplanır.
import sympy as sp
x, y = sp.symbols("x y")
print(sp.integrate(x * sp.exp(x), x))
print(sp.integrate(1 / (1 + x**2), x))
print(sp.integrate(sp.sin(x)**2, (x, 0, sp.pi)))
print(sp.integrate(sp.exp(-x**2), (x, -sp.oo, sp.oo)))
print(sp.integrate(1 / x**2, (x, 1, sp.oo)))
print(sp.integrate(x * y, (y, 0, x), (x, 0, 1)))Çıktı:
(x - 1)*exp(x)
atan(x)
pi/2
sqrt(pi)
1
1/8İlk satır kısmi integrasyonla bulunan \(\int x e^x\,dx = (x - 1)e^x\), ikinci satır da temel \(\int \frac{dx}{1 + x^2} = \arctan x\) integralidir. Ardından \(\int_0^\pi \sin^2 x\,dx = \pi/2\), Gauss integrali \(\int_{-\infty}^{\infty} e^{-x^2}\,dx = \sqrt{\pi}\) ve yakınsak bir genelleştirilmiş integral gelir. Son satır, önce \(y\)’ye göre \(0\)’dan \(x\)’e, sonra \(x\)’e göre \(0\)’dan \(1\)’e alınan
\[\int_0^1 \int_0^x xy\,dy\,dx = \int_0^1 \frac{x^3}{2}\,dx = \frac{1}{8}\]
katlı integralidir.
sp.integrate(f, x) bir ilkel fonksiyon döndürür, \(+C\) sabitini yazmaz. Ayrıca \(\int dx/x\) için log(x) verir, mutlak değer koymaz. Bu, \(x > 0\) aralığında doğru bir ilkeldir; \(x < 0\) için ilkel \(\ln|x|\)’tir. Bulunan ilkeli türev alarak denetlemek iyi bir alışkanlıktır:
import sympy as sp
x = sp.symbols("x")
F = sp.integrate(x * sp.exp(x), x)
print(sp.simplify(sp.diff(F, x) - x * sp.exp(x)))
print(sp.integrate(1 / x, x))Çıktı:
0
log(x)Her integralin ilkeli temel fonksiyonlarla yazılamaz. SymPy böyle durumlarda ya integralle tanımlanan özel bir fonksiyon kullanır ya da integrali hesaplanmamış olarak geri verir:
import sympy as sp
x = sp.symbols("x")
print(sp.integrate(sp.exp(-x**2), x))
print(sp.integrate(sp.sin(x) / x, x))
I = sp.Integral(x**x, (x, 0, 1))
print(I.doit())
print(I.evalf(20))Çıktı:
sqrt(pi)*erf(x)/2
Si(x)
Integral(x**x, (x, 0, 1))
0.78343051071213440706erf hata fonksiyonu, \(\operatorname{erf}(x) = \frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2}\,dt\) ile; Si sinüs integrali de \(\operatorname{Si}(x) = \int_0^x \frac{\sin t}{t}\,dt\) ile tanımlanır. \(\int_0^1 x^x\,dx\) için kapalı bir biçim bulunamadı ve doit integrali olduğu gibi döndürdü. evalf(20) ise integrali sayısal olarak 20 anlamlı basamakla hesapladı. Sayısal integrali Sayısal Türev ve İntegral bölümünde ayrıntısıyla göreceğiz.
Örnek 9.3 (İki Eğri Arasındaki Alan) \(y = x^2\) parabolü ile \(y = x + 2\) doğrusu arasında kalan sınırlı bölgenin alanını SymPy ile hesaplayınız.
Çözüm
Önce iki eğrinin kesiştiği noktaları buluruz. Bunun için sp.solve fonksiyonunu kullanacağız. Fonksiyonu Denklem Çözmek kısmında ayrıntısıyla göreceğiz; burada denklemin çözümlerini bir liste olarak verir. sp.Eq(top, bottom) de iki tarafın eşitliğini, yani \(x + 2 = x^2\) denklemini yazar. Sonra kesişimler arasında üstteki eğriden alttakini çıkarıp integral alırız (bkz. Analiz 2):
import sympy as sp
x = sp.symbols("x")
top = x + 2
bottom = x**2
a, b = sp.solve(sp.Eq(top, bottom), x)
print("kesişimler:", a, b)
print("orta noktada fark:", (top - bottom).subs(x, (a + b) / 2))
area = sp.integrate(top - bottom, (x, a, b))
print("alan:", area, "=", area.evalf())Çıktı:
kesişimler: -1 2
orta noktada fark: 9/4
alan: 9/2 = 4.50000000000000Eğriler \(x = -1\) ve \(x = 2\)’de kesişir. Aradaki bir noktada farkın pozitif çıkması, bu aralıkta doğrunun parabolün üstünde kaldığını doğrular. Elle hesap aynı sonucu verir:
\[\int_{-1}^{2} (x + 2 - x^2)\,dx = \left[\frac{x^2}{2} + 2x - \frac{x^3}{3}\right]_{-1}^{2} = \frac{10}{3} + \frac{7}{6} = \frac{9}{2}.\]
\(\blacksquare\)
9.8 Seriler ve Toplamlar
Bir fonksiyonun bir nokta etrafındaki Taylor açılımı sp.series(f, x, a, n) ile bulunur. Sonuç \((x - a)^{n-1}\) terimine kadar olan açılım ile bir \(O\big((x - a)^n\big)\) kalan teriminden oluşur. removeO metodu kalan terimini atar ve geriye Taylor polinomu kalır (bkz. Analiz 2).
import sympy as sp
x = sp.symbols("x")
s = sp.series(sp.exp(x), x, 0, 6)
print(s)
print(s.removeO())
print(sp.series(sp.tan(x), x, 0, 8))
print(sp.series(sp.sqrt(x), x, 4, 3))Çıktı:
1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + O(x**6)
x**5/120 + x**4/24 + x**3/6 + x**2/2 + x + 1
x + x**3/3 + 2*x**5/15 + 17*x**7/315 + O(x**8)
1 - (x - 4)**2/64 + x/4 + O((x - 4)**3, (x, 4))n parametresi terim sayısını değil, kalanın kuvvetini belirler: sp.series(sp.exp(x), x, 0, 6) en çok \(x^5\) terimini içerir. Son satır \(\sqrt{x}\)’in \(4\) etrafındaki açılımıdır. SymPy \(2 + \frac{x - 4}{4}\) toplamını \(1 + \frac{x}{4}\) olarak yazdığından açılım ilk bakışta tanıdık görünmez; ama \(\sqrt{x} \approx 2 + \frac{x-4}{4} - \frac{(x-4)^2}{64}\) ile aynıdır.
Sonlu ve sonsuz toplamlar sp.summation ile hesaplanır:
import sympy as sp
k, n = sp.symbols("k n", integer=True, positive=True)
s = sp.summation(k**2, (k, 1, n))
print(s)
print(sp.factor(s))
print(sp.summation(1 / k**2, (k, 1, sp.oo)))
print(sp.summation(sp.Rational(1, 2)**k, (k, 0, sp.oo)))Çıktı:
n**3/3 + n**2/2 + n/6
n*(n + 1)*(2*n + 1)/6
pi**2/6
2İlk iki satır \(1^2 + 2^2 + \cdots + n^2 = \frac{n(n+1)(2n+1)}{6}\) formülüdür. Ardından Basel problemi olarak bilinen \(\sum 1/k^2 = \pi^2/6\) eşitliği ve \(1/2\) oranlı geometrik serinin toplamı gelir.
Örnek 9.4 (e Sayısına Taylor Polinomlarıyla Yaklaşmak) \(e^x\) fonksiyonunun \(0\) etrafındaki \(n\)’inci Taylor polinomu \(T_n\) olsun. \(|e - T_n(1)| < 10^{-4}\) eşitsizliğini sağlayan en küçük \(n\)’yi SymPy ile bulunuz.
Çözüm
\(T_n(1) = \sum_{k=0}^{n} 1/k!\) değerini kesir olarak hesaplayıp hatayı evalf ile ölçelim:
import sympy as sp
x = sp.symbols("x")
for n in range(4, 9):
T = sp.series(sp.exp(x), x, 0, n + 1).removeO()
value = T.subs(x, 1)
error = (sp.E - value).evalf(5)
print(n, value, error)Çıktı:
4 65/24 0.0099485
5 163/60 0.0016152
6 1957/720 0.00022627
7 685/252 2.7860e-5
8 109601/40320 3.0586e-6Hata \(n = 6\) için \(2{,}26 \cdot 10^{-4}\), \(n = 7\) için \(2{,}79 \cdot 10^{-5}\)’tir. İstenen en küçük değer \(n = 7\)’dir ve \(e \approx T_7(1) = \frac{685}{252}\) olur. Hata kesin bir kesir olan \(T_n(1)\) ile sembolik \(e\) arasındaki farktan hesaplandığından, kendisi bir yuvarlama hatası taşımaz. Benzer bir soruyu \(5 \cdot 10^{-5}\) duyarlıkla, Lagrange kalanıyla ve hesap yapmadan önce çözen bir örnek için bkz. Analiz 2. Oradaki \(e - T_n(1) < \frac{3}{(n+1)!}\) sınırı bizim \(10^{-4}\) duyarlığımızda da \(n = 7\) verir: sınırın \(10^{-4}\)’ten küçük olması için \((n+1)! > 30\,000\) gerekir ve \(7! = 5040 < 30\,000 < 40\,320 = 8!\) olduğundan en küçük seçim \(n + 1 = 8\)’dir.
\(\blacksquare\)
9.9 Denklem Çözmek
Alan örneğinde kısaca kullandığımız sp.solve, SymPy’nin en çok başvurulan fonksiyonlarından biridir. sp.solve(expr, x), \(\text{expr} = 0\) denkleminin çözümlerini bir liste olarak verir. Denklemi iki taraflı yazmak için sp.Eq(lhs, rhs) kullanılır. Python’da = atama, == de yapısal eşitlik olduğundan denklem için ayrı bir nesne gerekir.
import sympy as sp
x, a, b, c = sp.symbols("x a b c")
print(sp.solve(x**2 - 5*x + 6, x))
print(sp.solve(sp.Eq(x**3, 8), x))
print(sp.solve(x**2 + 1, x))
print(sp.solve(a*x**2 + b*x + c, x))Çıktı:
[2, 3]
[2, -1 - sqrt(3)*I, -1 + sqrt(3)*I]
[-I, I]
[(-b - sqrt(-4*a*c + b**2))/(2*a), (-b + sqrt(-4*a*c + b**2))/(2*a)]I, SymPy’nin sanal birimidir: \(x^3 = 8\) denkleminin bir reel, iki karmaşık kökü vardır. Son satır ikinci derece denklemin kök formülüdür. Ancak dikkat: SymPy burada \(a \ne 0\) olduğunu sessizce varsaydı. solve “genel” durumu çözer; parametrelerin özel değerlerini ayrıca incelemek bize düşer.
Çözümleri belli bir kümede aramak için sp.solveset daha uygundur. domain seçeneği çözüm kümesini belirler ve sonuç her zaman bir küme olur; eşitsizlikler de çözülebilir:
import sympy as sp
x = sp.symbols("x")
print(sp.solveset(x**2 + 1, x, domain=sp.S.Reals))
print(sp.solveset(sp.sin(x), x, domain=sp.Interval(0, 10)))
print(sp.solveset(x**2 - 4 > 0, x, domain=sp.S.Reals))Çıktı:
EmptySet
{0, pi, 2*pi, 3*pi}
Union(Interval.open(-oo, -2), Interval.open(2, oo))\(x^2 + 1 = 0\) denkleminin reel çözümü yoktur (EmptySet, boş küme). \([0, 10]\) aralığında \(\sin x = 0\) denkleminin çözümleri \(0, \pi, 2\pi, 3\pi\)’dir. \(x^2 > 4\) eşitsizliğinin çözümü de \((-\infty, -2) \cup (2, \infty)\) kümesidir.
Denklem sistemleri için denklemler bir liste, bilinmeyenler ikinci bir liste olarak verilir. Bir çember ile bir doğrunun kesişimini bulalım: \(x^2 + y^2 = 5\), \(y = x + 1\).
import sympy as sp
x, y = sp.symbols("x y")
eqs = [sp.Eq(x**2 + y**2, 5), sp.Eq(y, x + 1)]
sols = sp.solve(eqs, [x, y], dict=True)
print(sols)
for s in sols:
print(s[x], s[y], [e.subs(s) for e in eqs])Çıktı:
[{x: -2, y: -1}, {x: 1, y: 2}]
-2 -1 [True, True]
1 2 [True, True]dict=True seçeneği her çözümü {x: ..., y: ...} biçiminde bir sözlük olarak verir. Böylece çözüm doğrudan subs’a verilebilir. Son iki satırdaki [True, True] listeleri, iki çözümün de iki denklemi sağladığını gösterir.
sp.solve fonksiyonunun bulduğu iki çözümdür.Doğrusal sistemlerde sp.linsolve çözüm kümesini verir; sonsuz çözüm varsa onu serbest değişkenlere bağlı olarak yazar:
import sympy as sp
x, y, z = sp.symbols("x y z")
eqs = [x + y + z - 6, x - y + 2*z - 5, 2*x + y - z - 1]
print(sp.linsolve(eqs, [x, y, z]))
print(sp.linsolve([x + y - 2, 2*x + 2*y - 4], [x, y]))Çıktı:
{(1, 2, 3)}
{(2 - y, y)}İlk sistemin tek çözümü \((1, 2, 3)\)’tür. İkinci sistemin iki denklemi aynı doğruyu verdiğinden sonsuz çözüm vardır: her \(y\) için \((2 - y, y)\) bir çözümdür.
Her denklemin kapalı biçimde bir çözümü yoktur. \(\cos x = x\) denkleminin kökü temel fonksiyonlarla yazılamaz ve sp.solve bunu NotImplementedError hatasıyla bildirir. Böyle durumlarda sp.nsolve, verilen bir başlangıç noktasından sayısal bir yöntemle kökü bulur:
import sympy as sp
x = sp.symbols("x")
try:
sp.solve(sp.cos(x) - x, x)
except NotImplementedError:
print("solve kapalı biçimde bir çözüm bulamadı.")
print(sp.nsolve(sp.cos(x) - x, x, 0.7))Çıktı:
solve kapalı biçimde bir çözüm bulamadı.
0.739085133215161Aynı kökü Newton-Raphson yöntemiyle elle hesaplayan örnek için bkz. Nümerik Analiz. Sayısal kök bulmayı bir sonraki bölümde SciPy ile ayrıntılı olarak ele alacağız.
9.10 Diferansiyel Denklem Çözmek
SymPy bazı diferansiyel denklem türlerini de kapalı biçimde çözer. Bilinmeyen fonksiyon sp.Function("y") ile tanımlanır ve y(x).diff(x) onun türevini gösterir. Denklem sp.Eq ile kurulur, sp.dsolve ile çözülür. Başlangıç koşulları ics sözlüğüyle verilir.
Birinci mertebeden doğrusal \(y' + y = x\) denklemini önce genel olarak, sonra \(y(0) = 1\) koşuluyla çözelim:
import sympy as sp
x = sp.symbols("x")
y = sp.Function("y")
ode = sp.Eq(y(x).diff(x) + y(x), x)
print(sp.dsolve(ode, y(x)))
sol = sp.dsolve(ode, y(x), ics={y(0): 1})
print(sol)
print(sp.checkodesol(ode, sol))Çıktı:
Eq(y(x), C1*exp(-x) + x - 1)
Eq(y(x), x - 1 + 2*exp(-x))
(True, 0)Genel çözüm \(y = C_1 e^{-x} + x - 1\) ailesidir; C1 keyfi sabittir. \(y(0) = 1\) koşulu \(C_1 - 1 = 1\), yani \(C_1 = 2\) verir. sp.checkodesol, çözümü denkleme koyar ve artığın sıfır olduğunu (True, 0) ile bildirir. Aynı sonuca elle de ulaşılır. Denklemi \(e^x\) integrasyon çarpanıyla çarparsak (bkz. Diferansiyel Denklemler)
\[(e^x y)' = x e^x \quad\Longrightarrow\quad e^x y = (x - 1)e^x + C_1\]
olur. Aşağıdaki şekil denklemin yön alanını ve aileden birkaç çözümü gösterir:
\(C_1 = 0\) için çözüm \(y = x - 1\) doğrusudur. Diğer bütün çözümler \(C_1 e^{-x} \to 0\) olduğundan \(x \to \infty\) iken bu doğruya yaklaşır. Kapalı biçimde çözülemeyen denklemleri Diferansiyel Denklemlerin Sayısal Çözümü bölümünde sayısal yöntemlerle çözeceğiz.
9.11 Matrisler
Lineer cebirin hesapları da sembolik olarak, yani kesin kesirlerle yapılabilir. sp.Matrix ile kurulan matrislerin girdileri kesin sayılar ya da ifadelerdir. Matris çarpımı * ya da @ ile yapılır; .T transpozu, .det() determinantı, .inv() tersi, .rank() rankı verir. NumPy ile Lineer Cebir bölümündeki np.linalg fonksiyonlarından farkı, sonuçların yaklaşık ondalık sayılar değil, kesin kesirler olmasıdır.
import sympy as sp
A = sp.Matrix([[2, 1, 0], [1, 3, 1], [0, 1, 2]])
print(A.det())
print(A.inv())
print(A * A.inv() == sp.eye(3))Çıktı:
8
Matrix([[5/8, -1/4, 1/8], [-1/4, 1/2, -1/4], [1/8, -1/4, 5/8]])
Truesp.eye(3) \(3 \times 3\) birim matristir. Ters matrisin girdileri kesin kesirler olduğundan \(AA^{-1} = I\) eşitliği yapısal olarak da tam tutar.
Satır indirgenmiş merdiven biçim (bkz. Lineer Cebir) rref metoduyla, sıfır uzayının bir tabanı da nullspace metoduyla bulunur:
import sympy as sp
B = sp.Matrix([[1, 2, 1, 1],
[2, 4, 0, 6],
[1, 2, 3, -1]])
R, pivots = B.rref()
print(R)
print(pivots)
print([v.T for v in B.nullspace()])Çıktı:
Matrix([[1, 2, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
(0, 2, 3)
[Matrix([[-2, 1, 0, 0]])]rref iki şey döndürür: indirgenmiş matrisi ve pivot sütunlarının indislerini. İndisler Python alışkanlığıyla \(0\)’dan başlar; yani pivotlar 1., 3. ve 4. sütunlardadır. Pivotsuz tek sütun ikincisi olduğundan sıfır uzayı tek boyutludur. Tabanı \((-2, 1, 0, 0)\) vektörüdür: \(x_1 + 2x_2 = 0\), \(x_3 = x_4 = 0\). Vektörleri tek satırda yazdırmak için transpozlarını (v.T) yazdırdık.
Özdeğer ve özvektörler (bkz. Lineer Cebir) için charpoly, eigenvals, eigenvects ve diagonalize metodları vardır. Karakteristik polinomun değişkeni olan \(\lambda\) sembolünü lam adlı değişkene bağlıyoruz, çünkü lambda Python’da ayrılmış bir sözcüktür (bkz. Tanım 2.12):
import sympy as sp
lam = sp.symbols("lambda")
A = sp.Matrix([[4, 1], [2, 3]])
p = A.charpoly(lam).as_expr()
print(p, "=", sp.factor(p))
print(A.eigenvals())
for value, mult, vectors in A.eigenvects():
print(value, mult, [v.T for v in vectors])
P, D = A.diagonalize()
print(P, D)
print(P * D * P.inv() == A)Çıktı:
lambda**2 - 7*lambda + 10 = (lambda - 5)*(lambda - 2)
{5: 1, 2: 1}
2 1 [Matrix([[-1/2, 1]])]
5 1 [Matrix([[1, 1]])]
Matrix([[-1, 1], [2, 1]]) Matrix([[2, 0], [0, 5]])
TrueKarakteristik polinom \(\lambda^2 - 7\lambda + 10 = (\lambda - 2)(\lambda - 5)\)’tir. eigenvals her özdeğeri cebirsel katıyla birlikte bir sözlük olarak verir. eigenvects her özdeğer için katı ve öz uzayın bir tabanını listeler. diagonalize, \(A = PDP^{-1}\) olacak biçimde \(P\) ve köşegen \(D\) matrislerini döndürür. \(P\)’nin sütunları özvektörlerdir; SymPy \((-1/2, 1)\) özvektörünü \(2\) ile çarpıp \((-1, 2)\) olarak almıştır (bkz. Lineer Cebir).
Girdileri sembol olan matrislerle de çalışılabilir. \(\theta\) ve \(\varphi\) açılı iki dönme matrisinin çarpımını hesaplayalım:
import sympy as sp
theta, phi = sp.symbols("theta phi", real=True)
def rot(a):
return sp.Matrix([[sp.cos(a), -sp.sin(a)],
[sp.sin(a), sp.cos(a)]])
M = sp.simplify(rot(theta) * rot(phi))
for row in M.tolist():
print(row)
print(M == rot(theta + phi))
print(sp.simplify(rot(theta).det()))Çıktı:
[cos(phi + theta), -sin(phi + theta)]
[sin(phi + theta), cos(phi + theta)]
True
1İki dönmenin bileşkesi, açıları toplanmış dönmedir. simplify bunu görmek için \(\cos(\theta + \varphi) = \cos\theta\cos\varphi - \sin\theta\sin\varphi\) toplam formüllerini kullandı. Dönme matrisinin determinantı da \(\cos^2\theta + \sin^2\theta = 1\)’dir.
9.12 Sembolikten Sayısala: lambdify
Sembolik bir sonuç çoğu zaman yolun yarısıdır. Bulduğumuz bir türevi binlerce noktada hesaplamak, grafiğini çizmek ya da bir sayısal yöntemde kullanmak isteriz. subs ile evalf bunu yapabilir, ama her noktada ifade ağacını baştan dolaştığı için yavaştır.
Tanım 9.6 (Lambdify) sp.lambdify(x, expr, "numpy") çağrısı, expr ifadesini, x yerine bir sayı ya da bir NumPy dizisi alan sıradan bir Python fonksiyonuna çevirir. Üretilen fonksiyonda sp.sin, sp.exp gibi SymPy fonksiyonlarının yerini np.sin, np.exp gibi NumPy karşılıkları alır.
Yani lambdify sembolik dünyadan sayısal dünyaya geçen köprüdür: formülü SymPy ile bir kez türetiriz, sonra NumPy hızıyla bir dizinin bütün elemanları üzerinde hesaplarız (bkz. Vektörizasyon ve Broadcasting).
import numpy as np
import sympy as sp
x = sp.symbols("x")
expr = sp.diff(x**2 * sp.exp(-x), x)
print(expr)
f = sp.lambdify(x, expr, "numpy")
xs = np.linspace(0, 4, 5)
print(f(xs))
print(f(1.0), expr.subs(x, 1).evalf())Çıktı:
-x**2*exp(-x) + 2*x*exp(-x)
[ 0. 0.36787944 0. -0.14936121 -0.14652511]
0.36787944117144233 0.367879441171442f, \(x^2 e^{-x}\) fonksiyonunun \((2x - x^2)e^{-x}\) türevini hesaplayan bir NumPy fonksiyonudur. Bir dizi alınca bütün elemanlarda birden hesaplar; tek bir sayıda ise subs ve evalf ile aynı değeri verir. Hız farkını ölçelim:
import time
import numpy as np
import sympy as sp
x = sp.symbols("x")
expr = sp.exp(-x**2) * sp.cos(3 * x)
xs = np.linspace(0, 2, 2000)
t0 = time.perf_counter()
slow = np.array([float(expr.subs(x, v)) for v in xs])
t1 = time.perf_counter()
f = sp.lambdify(x, expr, "numpy")
t2 = time.perf_counter()
fast = f(xs)
t3 = time.perf_counter()
print(f"subs ile döngü : {1000 * (t1 - t0):.0f} ms")
print(f"lambdify (üretim) : {1000 * (t2 - t1):.0f} ms")
print(f"lambdify (hesap) : {1000 * (t3 - t2):.2f} ms")
print("en büyük fark :", np.max(np.abs(slow - fast)))Çıktı:
subs ile döngü : 962 ms
lambdify (üretim) : 110 ms
lambdify (hesap) : 0.05 ms
en büyük fark : 2.220446049250313e-16Süreler bu notların yazıldığı bilgisayarda ölçülmüştür; başka bir bilgisayarda, hatta aynı bilgisayarda her çalıştırmada biraz farklı çıkar. Oranlar ise kabaca aynı kalır. subs ile 2000 noktalık döngü yaklaşık bir saniye sürdü. lambdify ile üretilen fonksiyon aynı hesabı milisaniyenin çok altında, binlerce kat hızlı yaptı. Fonksiyonu üretmenin maliyeti bir kerelik bir bedeldir: fonksiyon bir kez üretilir, sonra istenen kadar çağrılır. Sonuçlar ise neredeyse tıpatıp aynıdır; en büyük fark son birkaç ikili basamaktaki yuvarlamadan gelir.
- Formülü SymPy ile türetin (
diff,integrate,series,solve, …). sp.lambdify(x, expr, "numpy")ile ifadeyi bir Python fonksiyonuna çevirin.- Fonksiyonu NumPy dizileriyle çağırıp sonuçları hesaplayın ya da çizin.
Tarifi sinüsün Taylor polinomlarına uygulayalım. \(T_1, T_3, T_5, T_7\) polinomlarını SymPy ile bulup NumPy ile hesaplayacak ve Matplotlib ile Grafik Çizimi bölümündeki gibi \(\sin x\) ile aynı eksende çizeceğiz:
import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
x = sp.symbols("x")
xs = np.linspace(-2 * np.pi, 2 * np.pi, 400)
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(xs, np.sin(xs), color="black", lw=2.5, label="sin x")
for n in [1, 3, 5, 7]:
T = sp.series(sp.sin(x), x, 0, n + 1).removeO() # 1. adım
T_num = sp.lambdify(x, T, "numpy") # 2. adım
ax.plot(xs, T_num(xs), label=f"T{n}") # 3. adım
print(f"T{n}(x) =", T)
ax.set_ylim(-3, 3)
ax.axhline(0, color="gray", lw=0.5)
ax.legend()
fig.savefig("taylor-sin.png", dpi=150)Çıktı:
T1(x) = x
T3(x) = -x**3/6 + x
T5(x) = x**5/120 - x**3/6 + x
T7(x) = -x**7/5040 + x**5/120 - x**3/6 + xKod, polinomları yazdırır ve grafiği taylor-sin.png dosyasına kaydeder. Aşağıdaki şekil aynı grafiğin aynı verilerle çizilmiş hâlidir:
Bütün polinomlar \(0\) yakınında \(\sin x\) ile çakışır; derece arttıkça çakışmanın sürdüğü aralık genişler. Uzakta ise her polinom en yüksek dereceli teriminin işaretine göre \(\pm\infty\)’a gider. Taylor polinomu yerel bir yaklaşımdır ve hatası Lagrange kalanıyla kestirilir (bkz. Analiz 2).
9.13 Alıştırmalar
Aşağıdaki alıştırmaları önce kendiniz çözmeye çalışın; her çözümde kod ve çıktısı verilmiştir.
Alıştırma 9.1 (Kosinüsün Dördüncü Kuvveti) \(\cos^4 x - \sin^4 x = \cos 2x\) özdeşliğini SymPy ile doğrulayınız.
Çözüm
Adım 1. İki tarafı ayrı ifadeler olarak yazıp farkı sadeleştiririz. Bu, özdeşlik doğrulama tarifidir (bkz. Örnek 9.2).
Adım 2. Sol tarafı çarpanlarına ayırarak özdeşliğin neden doğru olduğunu da görelim:
import sympy as sp
x = sp.symbols("x")
lhs = sp.cos(x)**4 - sp.sin(x)**4
rhs = sp.cos(2*x)
print(sp.simplify(lhs - rhs))
print(sp.factor(lhs))Çıktı:
0
(-sin(x) + cos(x))*(sin(x) + cos(x))*(sin(x)**2 + cos(x)**2)Adım 3. Fark \(0\) olduğundan özdeşlik doğrudur. Çarpanlara ayrılmış biçim elle ispatı da verir. factor, \(a^4 - b^4 = (a - b)(a + b)(a^2 + b^2)\) ayrışımını yaptı; ilk iki çarpanın çarpımı \(\cos^2 x - \sin^2 x = \cos 2x\), üçüncü çarpan da \(1\)’dir:
\[\cos^4 x - \sin^4 x = (\cos^2 x - \sin^2 x)(\cos^2 x + \sin^2 x) = \cos 2x \cdot 1.\]
\(\blacksquare\)
Alıştırma 9.2 (Köklü Bir Limit) \(\displaystyle\lim_{x \to 0} \frac{\sqrt{x + 4} - 2}{x}\) limitini SymPy ile hesaplayınız ve sonucu eşlenikle çarparak doğrulayınız.
Çözüm
Adım 1. Limiti doğrudan hesaplarız.
Adım 2. Payı eşleniği \(\sqrt{x + 4} + 2\) ile çarpınca paydaki köklerin gittiğini görürüz. Böylece kesir \(1/(\sqrt{x + 4} + 2)\) olur ve bunun limiti doğrudan alınır:
import sympy as sp
x = sp.symbols("x")
f = (sp.sqrt(x + 4) - 2) / x
print(sp.limit(f, x, 0))
conj = sp.sqrt(x + 4) + 2
print(sp.expand((sp.sqrt(x + 4) - 2) * conj))
print(sp.limit(1 / conj, x, 0))Çıktı:
1/4
x
1/4Adım 3. \((\sqrt{x + 4} - 2)(\sqrt{x + 4} + 2) = x\) olduğundan \(x \ne 0\) için
\[\frac{\sqrt{x + 4} - 2}{x} = \frac{1}{\sqrt{x + 4} + 2} \longrightarrow \frac{1}{2 + 2} = \frac{1}{4}\]
bulunur. İki yol aynı sonucu verir.
\(\blacksquare\)
Alıştırma 9.3 (Bir Kübik Polinomun Yerel Ekstremumları) \(f(x) = x^3 - 3x + 1\) fonksiyonunun yerel ekstremum noktalarını ve değerlerini SymPy ile bulunuz.
Çözüm
Adım 1. Kritik noktalar \(f'(x) = 0\) denkleminin çözümleridir.
Adım 2. Her kritik noktada ikinci türevin işaretine bakarız: \(f''(c) < 0\) ise yerel maksimum, \(f''(c) > 0\) ise yerel minimum vardır (bkz. Analiz 2).
import sympy as sp
x = sp.symbols("x")
f = x**3 - 3*x + 1
f1 = sp.diff(f, x)
f2 = sp.diff(f, x, 2)
critical = sp.solve(f1, x)
print("f'(x) =", f1, " kritik noktalar:", critical)
for c in critical:
print(c, "f'' =", f2.subs(x, c), "f =", f.subs(x, c))Çıktı:
f'(x) = 3*x**2 - 3 kritik noktalar: [-1, 1]
-1 f'' = -6 f = 3
1 f'' = 6 f = -1Adım 3. \(f'(x) = 3x^2 - 3 = 0\) denkleminin kökleri \(\pm 1\)’dir. \(f''(-1) = -6 < 0\) olduğundan \(x = -1\)’de yerel maksimum vardır ve değeri \(f(-1) = 3\)’tür. \(f''(1) = 6 > 0\) olduğundan \(x = 1\)’de yerel minimum vardır ve değeri \(f(1) = -1\)’dir.
\(\blacksquare\)
Alıştırma 9.4 (Kısmi İntegrasyonlu Bir Belirli İntegral) \(\displaystyle\int_0^1 x^2 e^x\,dx\) integralini kesin olarak hesaplayınız ve ondalık değerini veriniz.
Çözüm
Adım 1. Önce ilkeli, sonra belirli integrali hesaplarız.
Adım 2. Kesin sonucu evalf ile ondalığa çeviririz:
import sympy as sp
x = sp.symbols("x")
F = sp.integrate(x**2 * sp.exp(x), x)
print(F)
I = sp.integrate(x**2 * sp.exp(x), (x, 0, 1))
print(I, "=", I.evalf())Çıktı:
(x**2 - 2*x + 2)*exp(x)
-2 + E = 0.718281828459045Adım 3. İlkel \((x^2 - 2x + 2)e^x\)’tir; iki kez kısmi integrasyonla elle de bulunur. Sınırları koyarsak
\[\int_0^1 x^2 e^x\,dx = (1 - 2 + 2)e - 2 = e - 2 \approx 0{,}718281828\]
olur.
\(\blacksquare\)
Alıştırma 9.5 (Basit Kesirlerle İntegral) \(\displaystyle\int \frac{3x + 5}{(x - 1)(x + 2)}\,dx\) integralini önce apart ile basit kesirlere ayırarak, sonra integrate ile hesaplayınız.
Çözüm
Adım 1. apart ile kesri basit kesirlere ayırırız.
Adım 2. İntegrali alır ve bulunan ilkeli türev alarak denetleriz:
import sympy as sp
x = sp.symbols("x")
r = (3*x + 5) / ((x - 1) * (x + 2))
print(sp.apart(r))
F = sp.integrate(r, x)
print(F)
print(sp.simplify(sp.diff(F, x) - r))Çıktı:
1/(3*(x + 2)) + 8/(3*(x - 1))
8*log(x - 1)/3 + log(x + 2)/3
0Adım 3. Ayrışım
\[\frac{3x + 5}{(x - 1)(x + 2)} = \frac{8}{3(x - 1)} + \frac{1}{3(x + 2)}\]
biçimindedir. Gerçekten \(x = 1\) koyunca \(A = 8/3\), \(x = -2\) koyunca \(B = (-1)/(-3) = 1/3\) bulunur. Buradan
\[\int \frac{3x + 5}{(x - 1)(x + 2)}\,dx = \frac{8}{3}\ln|x - 1| + \frac{1}{3}\ln|x + 2| + C\]
olur. SymPy mutlak değerleri ve sabiti yazmaz; bu farkı İntegral kısmındaki uyarıda görmüştük (bkz. Analiz 2).
\(\blacksquare\)
Alıştırma 9.6 (Logaritmanın Taylor Polinomu) \(\ln(1 + x)\) fonksiyonunun \(0\) etrafındaki beşinci Taylor polinomu \(T_5\)’i bulunuz. \(T_5(0{,}1)\) değerinin \(\ln 1{,}1\) sayısından farkını hesaplayıp \(0{,}1^6/6\) sınırıyla karşılaştırınız.
Çözüm
Adım 1. series ile \(x^5\) terimine kadar açılımı alırız.
Adım 2. \(0{,}1\)’i sp.Rational(1, 10) olarak koyarız. Böylece \(T_5(0{,}1)\) kesin bir kesir olur; hesaplanan fark yuvarlamadan değil, yalnız açılımın kesilmesinden gelir:
import sympy as sp
x = sp.symbols("x")
T5 = sp.series(sp.log(1 + x), x, 0, 6).removeO()
print(T5)
h = sp.Rational(1, 10)
approx = T5.subs(x, h)
exact = sp.log(1 + h)
print(approx, "=", approx.evalf(12))
print("ln(1.1) =", exact.evalf(12))
print("hata =", (approx - exact).evalf(5))
print("sınır =", (h**6 / 6).evalf(5))Çıktı:
x**5/5 - x**4/4 + x**3/3 - x**2/2 + x
285931/3000000 = 0.0953103333333
ln(1.1) = 0.0953101798043
hata = 1.5353e-7
sınır = 1.6667e-7Adım 3. \(T_5(x) = x - \frac{x^2}{2} + \frac{x^3}{3} - \frac{x^4}{4} + \frac{x^5}{5}\)’tir. Hata yaklaşık \(1{,}535 \cdot 10^{-7}\)’dir ve \(0{,}1^6/6 \approx 1{,}667 \cdot 10^{-7}\) sınırının altında kalır. \(x = 0{,}1\) için açılım alterne bir seridir ve terimlerin mutlak değeri azalır. Bu yüzden hata, atılan ilk terimin mutlak değerinden küçüktür; atılan ilk terim de \(-x^6/6\)’dır (bkz. Analiz 2).
\(\blacksquare\)
Alıştırma 9.7 (Teleskopik Bir Toplam) \(\displaystyle\sum_{k=1}^{n} \frac{1}{k(k+1)}\) toplamının kapalı biçimini ve \(n \to \infty\) iken limitini SymPy ile bulunuz.
Çözüm
Adım 1. Genel terimi apart ile basit kesirlere ayırırız; ikinci argüman k, ayrıştırmanın hangi değişkene göre yapılacağını söyler.
Adım 2. Kısmi toplamı summation ile hesaplar, sonra limitini alır ve sonsuz toplamla karşılaştırırız:
import sympy as sp
k, n = sp.symbols("k n", integer=True, positive=True)
a = 1 / (k * (k + 1))
print(sp.apart(a, k))
S = sp.simplify(sp.summation(a, (k, 1, n)))
print(S)
print(sp.limit(S, n, sp.oo), sp.summation(a, (k, 1, sp.oo)))Çıktı:
-1/(k + 1) + 1/k
n/(n + 1)
1 1Adım 3. \(\frac{1}{k(k+1)} = \frac{1}{k} - \frac{1}{k+1}\) olduğundan toplamda ardışık terimler birbirini götürür:
\[\sum_{k=1}^{n} \left(\frac{1}{k} - \frac{1}{k+1}\right) = 1 - \frac{1}{n+1} = \frac{n}{n+1}.\]
Limiti \(1\)’dir; serinin toplamı da \(1\)’dir (bkz. Analiz 2).
\(\blacksquare\)
Alıştırma 9.8 (Doğrusal Olmayan Bir Sistem) \(x^2 + y^2 = 25\), \(xy = 12\) sisteminin bütün çözümlerini SymPy ile bulunuz ve doğrulayınız.
Çözüm
Adım 1. Denklemleri sp.Eq ile yazıp solve’a veririz. Sözlük seçeneği olmadan solve çözümleri \((x, y)\) demetleri olarak döndürür.
Adım 2. Her çözümü iki denkleme koyup hepsinin sağlandığını all ile denetleriz:
import sympy as sp
x, y = sp.symbols("x y")
eqs = [sp.Eq(x**2 + y**2, 25), sp.Eq(x * y, 12)]
sols = sp.solve(eqs, [x, y])
print(sols)
print(all(e.subs({x: a, y: b}) for a, b in sols for e in eqs))Çıktı:
[(-4, -3), (-3, -4), (3, 4), (4, 3)]
TrueAdım 3. Dört çözüm vardır: \((3, 4)\), \((4, 3)\), \((-3, -4)\), \((-4, -3)\). Elle de bulunur: \((x + y)^2 = 25 + 24 = 49\) ve \((x - y)^2 = 25 - 24 = 1\) olduğundan \(x + y = \pm 7\), \(x - y = \pm 1\)’dir. Dört işaret seçimi dört çözümü verir.
\(\blacksquare\)
Alıştırma 9.9 (Köşegenleştirmeyle Matris Kuvveti) \(A = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}\) matrisini köşegenleştirerek \(A^n\) için genel bir formül bulunuz ve \(n = 10\) için doğrulayınız.
Çözüm
Adım 1. diagonalize ile \(A = PDP^{-1}\) ayrışımını buluruz.
Adım 2. \(A^n = PD^nP^{-1}\)’dir ve köşegen matrisin kuvveti köşegen girdilerin kuvvetidir. \(n\)’yi negatif olmayan bir tam sayı sembolü olarak tanımlarsak SymPy formülü sembolik olarak hesaplar:
import sympy as sp
n = sp.symbols("n", integer=True, nonnegative=True)
A = sp.Matrix([[2, 1], [1, 2]])
P, D = A.diagonalize()
print(P, D)
An = P * D**n * P.inv()
print(An)
print(An.subs(n, 10) == A**10, A**10)Çıktı:
Matrix([[-1, 1], [1, 1]]) Matrix([[1, 0], [0, 3]])
Matrix([[3**n/2 + 1/2, 3**n/2 - 1/2], [3**n/2 - 1/2, 3**n/2 + 1/2]])
True Matrix([[29525, 29524], [29524, 29525]])Adım 3. Özdeğerler \(1\) ve \(3\)’tür, özvektörler \((-1, 1)\) ve \((1, 1)\)’dir. Buradan
\[A^n = \frac{1}{2}\begin{pmatrix} 3^n + 1 & 3^n - 1 \\ 3^n - 1 & 3^n + 1 \end{pmatrix}\]
bulunur. \(n = 10\) için \(3^{10} = 59049\) olduğundan köşegen girdiler \(29525\), diğerleri \(29524\)’tür ve doğrudan hesaplanan \(A^{10}\) ile aynıdır.
\(\blacksquare\)
Alıştırma 9.10 (Parametreli Bir Matrisin Tersinirliği) \(K = \begin{pmatrix} 1 & k & 0 \\ k & 1 & k \\ 0 & k & 1 \end{pmatrix}\) matrisinin hangi \(k\) değerleri için tersinir olmadığını SymPy ile bulunuz.
Çözüm
Adım 1. Bir kare matris ancak ve ancak determinantı sıfırdan farklıysa tersinirdir (bkz. Lineer Cebir). Determinantı \(k\)’ye bağlı bir ifade olarak hesaplarız.
Adım 2. Determinantı sıfır yapan \(k\) değerlerini solve ile bulur, bu değerlerde rankı da denetleriz:
import sympy as sp
k = sp.symbols("k")
K = sp.Matrix([[1, k, 0], [k, 1, k], [0, k, 1]])
d = K.det()
print(d)
roots = sp.solve(d, k)
print(roots)
for r in roots:
print(r, K.subs(k, r).rank())Çıktı:
1 - 2*k**2
[-sqrt(2)/2, sqrt(2)/2]
-sqrt(2)/2 2
sqrt(2)/2 2Adım 3. \(\det K = 1 - 2k^2\)’dir ve \(k = \pm\frac{\sqrt{2}}{2}\) için sıfır olur. Bu iki değerde rank \(2\)’ye düşer, yani \(K\) tersinir değildir. Diğer bütün \(k\) değerlerinde \(K\) tersinirdir.
\(\blacksquare\)
Alıştırma 9.11 (Sembolik Türevle Newton Yöntemi) \(x^3 - 2x - 5 = 0\) denkleminin kökünü Newton yöntemiyle \(x_0 = 2\)’den başlayarak beş adımda hesaplayınız. Türevi sp.diff ile bulunuz, fonksiyonları lambdify ile sayısallaştırınız ve sonucu sp.nsolve ile karşılaştırınız.
Çözüm
Adım 1. Newton yinelemesi \(x_{k+1} = x_k - f(x_k)/f'(x_k)\)’dir (bkz. Nümerik Analiz). Türevi elle yazmak yerine SymPy’ye bulduruyoruz.
Adım 2. \(f\) ve \(f'\)’yi lambdify ile Python fonksiyonlarına çevirip yinelemeyi sıradan bir döngüyle yaparız:
import sympy as sp
x = sp.symbols("x")
expr = x**3 - 2*x - 5
f = sp.lambdify(x, expr, "numpy")
fp = sp.lambdify(x, sp.diff(expr, x), "numpy")
xk = 2.0
for k in range(1, 6):
xk = xk - f(xk) / fp(xk)
print(k, xk)
print(sp.nsolve(expr, x, 2))Çıktı:
1 2.1
2 2.094568121104185
3 2.094551481698199
4 2.0945514815423265
5 2.0945514815423265
2.09455148154233Adım 3. \(f'(x) = 3x^2 - 2\)’dir. İlk adım \(x_1 = 2 - \frac{-1}{10} = 2{,}1\)’dir. Yineleme dördüncü adımda \(2{,}0945514815423265\) değerine oturur ve nsolve aynı kökü verir. Doğru basamak sayısının her adımda kabaca ikiye katlanması, Newton yönteminin karesel yakınsamasıdır.
\(\blacksquare\)
Bu bölümde SymPy ile kesin hesap yapmayı öğrendik: semboller ve varsayımlar, ifade dönüşümleri, limit, türev, integral, seri, denklem, diferansiyel denklem ve matris hesapları. Sembolik sonuçları da lambdify ile NumPy’ye taşıdık. Ancak birçok problemin kapalı biçimde çözümü yoktur. \(\cos x = x\) denkleminde gördüğümüz gibi o zaman sayısal yöntemlere başvururuz. Bir sonraki bölüm Kök Bulma ve Optimizasyon, denklemlerin köklerini ve fonksiyonların minimumlarını SciPy ile sayısal olarak bulmayı anlatıyor.