29 Aralık 2017 Cuma

Monty Hall problemi ve Monte Carlo çözümü

Fizik ve matematikte bir problemle karşılaştığımızda problemi kağıt-kalem yardımıyla 'analitik' olarak çözemediğimizde yada bu şekilde çözmeye girişmeden önce çözüme dair bir fikir edinmek adına ufak tefek deneyler yapmak için 'simulasyon' dediğimiz, bilgisayar yardımıyla bir takım sayısal yöntemler kullanırız. Bu yöntemlere dair fikir vermek adına ünlü bir olasılık problemini 'Monte Carlo' yöntemi olarak adlandırılan ve rastgele üretilmiş sayıları kullanan bir simulasyon ile çözmeye çalışacağız. Problemimiz ünlü 'Monty Hall' problemi!

Bir yarışmadasınız ve önünüzde bir tanesinin arkasında ödül, ikisinin ardında 'keçi' olan üç kapı bulunuyor. Kapılardan arkasında ödül olanı seçtiğinizde ödülü alıyor, diğer kapıları seçtiğinizde hiçbir şey kazanamıyorsunuz. Sunucu bir kapıyı seçmenizi istiyor ve siz de seçiyorsunuz. Ardından ödülün hangi kapı ardında olduğunu bilen sunucu kapılardan birini açıyor ve arkasından keçi çıkıyor. Size dönüp açtığı kapının arkasında keçi olduğunu gördüğünüz halde seçtiğiniz kapıyı değiştirip değiştirmeyeceğinizi soruyor. Sizce kazanmak için kapıyı değiştirmeli misiniz?

Problem 70'li yıllarda ortaya atılıyor ve değme matematikçilerin çözüm konusunda birbirine girmesine neden oluyor (hikayesi için tıklayınız). Çözümünün, birçok basit olasılık probleminin çözümünde olduğu gibi, 'sağduyuya karşı' görünen yapısı, problemi epey ilginç ve popüler kılmış durumda. Basit olasılık argümanları ile seçtiğimiz kapıyı değiştirmenin mi yoksa değiştirmemenin mi avantajlı çıkacağı gösterilebilir fakat biz bunu birbirinin aynısı birçok deneme ile benzetim (simulasyon) yaparak bulmaya çalışacağız.

Simulasyon için Python programlama dili ve Numpy kütüphanesini kullanacağız. Öncelikle ilgili kütüphaneleri çağıralım.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

Yapacağımız şey belirli sayıda örnek oyunlar oynayıp, her seferinde rastgele seçimler yaparak seçimimizi değiştirdiğimizde mi yoksa seçimimiz aynı kaldığında mı daha sık kazanıyoruz, bunu gözlemek!

Öncelikle ödülün arkasında olduğu kapıyı simule edelim. nsim tane örnek kapı üreteceğiz. Kapılar $0$,$1$ ve $2$ değerleri alabiliyor. Bunun için numpyın random.randint fonksiyonunu kullanacağız. Bu fonksiyon verdiğiniz bir aralıkta (bizim için $0$ ve $3$) size ile belirleyeceğiniz sayıda düzgün dağılmış (eşit olasılıklı) rastgele tam sayılar üretiyor.

In [2]:
"""
Fonksiyon: simulasyon_odul

Ardında ödül olan kapıyı rastgele belirleyen fonksiyon.

Örnek
-------
>>> print simulasyon_odul(3)
array([0, 0, 2])
"""

def simulasyon_odul(nsim):
    odul = np.random.randint(0,3, size = nsim)
    return odul

Şimdi de tahminlerimizi üreteceğimiz fonksiyonu tanımlayalım. Yine aynı şekilde nsim tane $0$,$1$ ve $2$ olan rastgele kapı numaraları üretiyoruz. Bizim durumumuzda her seferinde başlangıçta 'rastgele'bir kapı seçiyor; fakat tamamen farklı bir yönemle, örneğin hep aynı kapıyı da seçtirebiliriz.

In [3]:
"""
Foksiyon: simulasyon_tahmin

Yarışmacının seçeceği kapıyı rastgele seçen fonksiyon. 

Örnek
-------
>>> print simulasyon_tahmin(5)
array([0, 0, 0, 0, 0])
"""

def simulasyon_tahmin(nsim):
    tahmin = np.random.randint(0,3, size = nsim)
    return tahmin

Sırada sunucunun açacağı ve arkasında 'keçi' olan kapıyı belirlemek var. Bunun için sunucunun her zaman arkasında 'keçi' olan kapıyı açacağını bildiğimizden, bizim seçimimize göre eğer seçtiğimizin arkasında keçi varsa diğer keçili kapıyı, ödül varsa diğer iki keçili kapıdan birini belirliyoruz.

In [4]:
"""
Fonksiyon: keci_kapi

Arkasında ödül olmayan (keçi olan) ve kişinin seçmediği kapılardan birini seçen fonksiyon

Örnek
--------
>>> print keci_kapi(np.array([0, 1, 2]), np.array([1, 1, 1]))
>>> array([2, 2, 0])
"""
def keci_kapi(odul, tahmin):
    acik = []
    for i in range(len(odul)):
        kapilar = [0,1,2]
        if(odul[i] == tahmin[i]):
            kapilar.remove(odul[i])
            acik.append(np.random.choice(kapilar))
        else:
            kapilar.remove(odul[i])
            kapilar.remove(tahmin[i])
            acik.append(kapilar[0])
    return np.array(acik)

Son olarak da tahmin değiştirme fonksiyonumuzu (tahmin_degistir) yazıyoruz. Bu fonksiyon tahminimizi ve sunucunun açtığı, arkasında keçi olan kapıyı girdi olarak alıyor ve tahmini değiştirdiği kapıyı döndürüyor.

In [5]:
"""
Fonksiyon: tahmin_degistir

Arkasında keçi olan kapı açıldığında her zaman seçtiği kapıyı değiştiren strateji

Örnek
--------
>>> print tahmin_degistir(np.array([0, 1, 2]), np.array([1, 2, 1]))
>>> array([2, 0, 0])
"""
def tahmin_degistir(tahmin, kecikapilar):
    degistir = []
    for i in range(len(odul)):
        kapilar = [0,1,2]
        if(kecikapilar[i] == tahmin[i]):
            kapilar.remove(kecikapilar[i])
            degistir.append(np.random.choice(kapilar))
        else:
            kapilar.remove(tahmin[i])
            kapilar.remove(kecikapilar[i])
            degistir.append(kapilar[0])
    return np.array(degistir)

Son olarak oynadığımız oyunlardan kapıyı değiştirdiklerimizin ve değiştirmediklerimizin kaçında kazandığımızı hesaplamak için bir fonksiyon belirliyoruz.

In [6]:
"""
Fonksiyon: kazanma_yuzdesi
Kazanma yüzdesini hesaplayan fonksiyon

Örnek
---------
>>> print kazanma_yuzdesi(np.array([0, 1, 2]), np.array([0, 0, 0]))
33.333
"""
def kazanma_yuzdesi(tahmin, odul):
    return np.sum(tahmin == odul)/float(len(tahmin))

Fonksiyonlarımızı belirledikten sonra şimdi simulasyona başlayalım. nsim$=1000$ tane oyun oynayacağız. Tek tek ardından ödül olan ve tahmin ettiğimiz kapı değerlerini rastgele oluşturuyoruz. Sonrasında sunucunun seçeceği ve arkasında keçi çıkan kapıyı belirliyoruz. Sonrasında kapımızı değiştirip değiştirmeyeceğimize karar veriyoruz. Ardından elde ettiğimiz sonucu yazdırıyoruz.

In [7]:
nsim = 1000
odul = simulasyon_odul(nsim)
tahmin = simulasyon_tahmin(nsim)
kecikapilar = keci_kapi(odul,tahmin)
degistir = tahmin_degistir(tahmin,kecikapilar)

print("Eğer ilk seçtiği kapıyı değiştirmezse ", kazanma_yuzdesi(tahmin, odul))
print("Eğer ilk seçtiği kapıyı değiştirirse: ", kazanma_yuzdesi(degistir, odul))
Eğer ilk seçtiği kapıyı değiştirmezse  0.305
Eğer ilk seçtiği kapıyı değiştirirse:  0.695

Görüldüğü üzere kapıyı değiştirdiğinde %67'e yakın yani $\frac{2}{3}$ oranında kazanırken, değiştirmediğinde $\frac{1}{3}$ yakın bir oranda kazanıyor. Yani seçtiğiniz kapıyı değiştirmek en makul seçenek gibi görünüyor!

Bu oranın oynadığımız oyun sayısına bağlı olup olmadığını görmek için $1$'den başlayıp $10000$'e kadar $20$'şer artan sayılarda oyun oynayıp değiştirip değiştirmediği durumlarda kazanma yüzdesini çizdirelim, bakalım nasıl bir şey çıkacak ortaya.

In [8]:
degistirirse = []
degistirmezse = []
for nsim in range(1,10000,20):
    odul = simulasyon_odul(nsim)
    tahmin = simulasyon_tahmin(nsim)
    kecikapilar = keci_kapi(odul,tahmin)
    degistir = tahmin_degistir(tahmin,kecikapilar)
    degistirirse.append(kazanma_yuzdesi(degistir, odul))
    degistirmezse.append(kazanma_yuzdesi(tahmin, odul))

plt.figure(figsize=(14,10))    
plt.plot(range(1,10000,20),degistirirse, 'r.', label='Degistirirse')
plt.plot(range(1,10000,20),degistirmezse, 'b.', label='Degistirmezse')
plt.xlabel('Simulasyon sayısı')
plt.ylabel('Olasılık')
plt.legend()
plt.show()
    

Gördüğünüz gibi ilk başta düşük sayıda oyun oynadığımızda olasılıklar değişkenlik gösterirken, 10-15 oyundan sonra sabit olasılıklara yani 0.33 ve 0.67 civarında sabitleniyor. Elbette rastgele oyunlar oynadığımız için bu değerler her seferinde tam olarak 0.33 ve 0.66 çıkmıyor ve bu durum Monte Carlo yöntemlerinin genel özelliği. Her ne kadar 'teorik' değerlerle birebir aynı olmasa da ona çok çok yakın olduğunu görüyoruz.

Monte Carlo yöntemleri fizikte analitik çözümünü yapamadığımız birçok sistem hakkında bilgi edinmek ve birçok zor integrali hesaplamak için sık sık kullanılan yöntemlerden biri. Bu dönem bilgisayar mühendisliği bölümünden aldığım yüksek lisans Monte Carlo dersine konuya dair birçok şey öğrenme fırsatım oldu. Bu yazı yeni yılla birlikte bu dersin ardından öğrendiklerim üzerine hazırlayacağım Monte Carlo dizisi için ufak, eğlenceli bir açılış yazısı olsun istedim. İlerleyen günlerdeki yazılarda konuyu derinlemesine inceleyip aşağıdaki konu başlıklarına değinmeyi planlıyorum:

  • Monte Carlo yöntemlerine giriş
  • Markov Zincirleri ve Monte Carlo (MCMC)
  • Önem Örneklemesi ve Metropolis-Hastings Algoritması
  • Gibbs Örneklemesi
  • Hamilton Monte Carlo yöntemi
  • Kesin Örnekleme, Propp-Wilson yöntemi ve Ising Modeli

Yeni yılda görüşmek üzere!

*Yazıda kullanılan kodlar Boğaziçi Ünv. Fizik Bölümünde bu dönem aldığım 'Machine Learning in Physics' dersinde yaptığımız bir ödevden uyarlanmıştır.

28 Aralık 2017 Perşembe

(Aslında) Zamandan bağımsız Schrödinger Denklemi

Fizikte elimizdeki sistemi tanımlamak için 'durum' (state) dediğimiz objeler kullanırız. Sistemlerin davranışını ifade etmek için bu durumların sayıları (istatistiksel özellikleri), nasıl değiştikleri (dinamik özellikleri), bu durumlarla ilişkili büyüklüklerin korunup korunmadığı (simetri yasaları) gibi sorularla uğraşırız. Bu nedenle 'durum' kavramının fiziğin en temelinde yatan şeylerden biri diyebiliriz. Örneğin bir ideal gazı tanımlamak için onun durumunu, basınç ve hacim parametrelerini kullanabiliriz. Yani bana bahsi geçen gazın $P$ ve $V$ değerlerini verildiğinde ben o sistemi tanımlamak için (örneğin sıcaklık ve toplam enerjisini hesaplamak için) yeterli bilgiye sahibim diyebilirim. Gazın parametrelerinin değişimi için ise termodinamiğin dönüşüm yasalarını kullanıruz. Benzer şekilde klasik mekanikte de bir parçacığın durumunu belirtmek için o parçacığın konumunu ve momentumu belirtmem gerekir. Örneğin parçacık 3 boyutta hareket ediyorsa, her üç boyut için $x$, $y$ ve $z$, aynı zamanda her üç boyuttaki momentumu $p_x$, $p_y$ ve $p_z$'yi söylersek o parçacığı tanımlamış oluruz. Bu durumların zamanla değişimini ise lisede 'Newton Yasaları' ($F = ma$), biraz büyünce de 'Hamilton Denklemleri' adını verdiğimiz kurallar veriyor. Kuantum mekaniğine geldiğimizde ise, elimizdeki sistemi tanımlamak için 'dalga fonksiyonu' adını verdiğimiz durumlarımız bulunuyor. Dalga fonksiyonunun kendisi matematiksel bir yapı olsa da kendisi sistemin erişilebilir durumlarının olasılıklarını hesaplamak için kullanılıyor. Kuantum mekaniğinde durumların zamanla değişimini de Schrödinger Denklemi adını verdiğimiz denklem ifade ediyor.

Kuantum mekaniğinde sistemi tanımlamak için kullandığımız 'durum'lar klasik mekanikteki parametrelerden biraz farklı davranıyorlar; daha farklı bir matematiksel yapıya sahipler aslında. Durumlar, klasik mekanikteki konum ve momentum değerleri gibi bir koordinat sisteminde ifade ettiğimiz 'nokta'lardan farklı olarak 'vektör uzayı' dediğimiz bir matematiksel yapının elemanları, dolayısıyla birer 'vektör'ler. Durumları zamanla değiştirmek için de vektörleri birbirine dönüştürmemize yani bir 'dönüşüm operatörüne (işlemci)' ihtiyacımız var. Operatör'leri kabaca matrislerin genel halleri gibi düşünebiliriz; bir vektörü matrisle çarparak farklı bir vektör elde ediyoruz gibi... Yani sistemin zamanla evrilmesi demek durum vektörlerimizin 'zamanda evriltme (öteleme) operatörü' dediğimiz bir bir operatör aracılığıyla dönüşmesi anlamına geliyor.

Kuantum mekaniği kitaplarında sistemin durumunun değişimini formüle etmek için öncelikle 'gökten inmiş' gibi davranılan Schrödinger denklemi ile başlanıp, sonrasında buradan yola çıkıp yukarıda bahsettiğim 'zamanda ötemele operatörü' tanımlanıyor. Şu gökten inme olayını ortadan kaldırmak adına olaya tersten baktığımızda aslında Schrödinger Denklemi'nin nereden geldiği ve ne olduğunu daha iyi anlama fırsatımız olacak. Hadi başlayalım!

Elimizde bir 'zamanda ötemele operatörü' $U(t)$ olsun ve bu operatörü herhangi bir $t_0$ anında durum vektörü $\left|\psi(t_0)\right\rangle$'ne uyguladığımızda bize o durumun $t$ kadar sonraki halini versin, yani sistemi zamanda $t$ kadar ötelesin/evriltsin.
$$U(t) \left|\psi(t_0)\right\rangle =\left|\psi(t_0 + t)\right\rangle$$
Bu dönüşümün öncelikle durumlar arasındaki "ilişkiyi" korumasını istiyoruz. Yani durumlar başlangıçta birbirinden tamamen farklıysa, örneğin birbirine dikse (orthogonal), dönüşüm sonrasında da yani $t$ kadar süre sonra da durumların birbirine dik kalmalarını istiyoruz. Daha genel olarak durumların birbirlerine 'benzerliklerinin' bir ölçüsü olarak ifade edebileceğimiz iç çarpımlarının (inner product) değişmemesi istiyoruz. Bu koşul bizim dönüşüm operatörümüz için ektra bir özellik gerektirecek, o da operatörün unitary (Türkçesi: birimcil - ben İngilizce olan 'unitary'yi kullanacağım) olması:
$$\left\langle\phi\right|\psi\rangle = \left\langle\phi\right|U^\dagger U | \psi\rangle \hspace{1cm} \Rightarrow \hspace{1cm} U^\dagger U = I \text{  (unitary)} $$
Yukarıda $\left|\phi\right\rangle$ ve $\left|\psi\right\rangle$ gibi iki vektörün ilk baştaki iç çarpımları (braket notasyonu için bknz.), $U$ ve $U^\dagger$ yani $U$'nun Hermitsel eşleniği (kompleks eşlenik) ile dönüşümü sonucunda elde edilen iki vektörün iç çarpımına eşit olması koşulunun $U$'nun unitary bir operatör olduğu sonucuna vardığını görüyoruz. Yani kısacası bir durumu zamanda $t$ kadar ilerletmenin yolu onu $U(t)$ operatörü ile işleme sokmak.
$$U(t) \left|\psi(0)\right\rangle =\left|\psi(t)\right\rangle$$
Bir durumu $t = 0$ kadar ilerletmek ona hiç bir şey yapmamakla eşdeğer, yani birim işlemci (identity) uygulamayla eşdeğer:
$$U(0) = I  \hspace{1cm} \Rightarrow \hspace{1cm} U^\dagger(0) = I$$
Elimimizdeki durumu zamanda $\epsilon$ kadar zamanda ötelemek için gereken operatörün şu formda olmasını bekliyoruz:
$$U(\epsilon) = I + \epsilon G \hspace{1cm} \Rightarrow \hspace{1cm} U^\dagger(\epsilon) = I + \epsilon G^\dagger$$
$G$ durumu $\epsilon$ kadar ilerleten operatör (kısa bir zaman dönüşümü ürettiği için üreteç (generator) olarak adlandırılıyor). Yeni oluşturduğumuz $U(\epsilon)$ operatörünün unitary olması gerekiyor:
$$(I + \epsilon G^\dagger)(I + \epsilon G) = I \\ \epsilon(G^\dagger + G) = 0 \\ G^\dagger = - G$$
Son satırda elimizdeki $G$ operatörünün anti-Hermitian özelliği olduğunu görüyoruz. Fakat kuantum mekaniğinde operatörlerin Hermitian özelliğe sahip olmasını istiyoruz (özdeğerleri real olması gibi özelliklerinden dolayı); dolayısıyla operatörümüzü alternatif olarak şöyle tanımlarsak istediğimiz koşul sağlanacak:
$$U(\epsilon) = I - \frac{i}{\hbar} \epsilon H \hspace{1cm} \Rightarrow \hspace{1cm} U^\dagger(\epsilon) = I + \frac{i}{\hbar} \epsilon H^\dagger \\ \boxed{H = H^\dagger}$$
$G$ operatörünü yeniden adlandırarak $H$ ismini verdik çünkü göreceğiz ki zamanda ötemenin üreteci sistemin enerjisiyle ilişkili Hamiltonyen operatörü.$\hbar$ Planck sabiti ve bunu denklemdeki boyutların tutması için gerekli bir sabit gibi düşünebilirsiniz.

Şu ana kadar yaptığımız şey, elimizdeki sistemi $\epsilon$ kadar küçük bir zaman için evriltmek ve sonucunu görmek. Buradan yola çıkıp sistemi daha uzun zamanlar evriltebilir miyiz peki? Elimizdeki kısa süreli evriltme işlemini üst üste uygulayarak daha uzun sürelerdeki değişimi elde edebiliriz. O halde yapalım; $n \rightarrow \infty $ ve $\epsilon$  çok küçük olup, çarpımları $\epsilon n = t$ olduğunda:
$$(I- \frac{i}{\hbar} \epsilon H)^n \hspace{0.5cm} \Rightarrow \hspace{0.5cm} e^{-\frac{i}{\hbar}Ht}$$
Eksponensiyel fonksiyonun tanımı gereği, $\epsilon$ kadar küçük zaman ötelemesini $n$ kere uygulayıp sistemi toplamda $t$ kadar ötelemek için $e^{-\frac{i}{\hbar}Ht}$ gibi bir operatör kullanmam gerekiyor. (Bir operatörün eksponansiyeli için bknz.)

Şimdi vurucu kısma geldik. Elimizde en başta bahsettiğimiz $\left|\psi(t)\right\rangle$ durumumuz var ve bunu $\Delta t$ kadar zamanda evriltmek istiyoruz. Yapacağımız işlem şu şekilde olacak:
$$\left|\psi(t + \Delta t)\right\rangle = e^{-\frac{i}{\hbar}H \Delta t} \left|\psi(t)\right\rangle$$
Eğer $\Delta t$ çok büyük değilse eksponansiyel fonksiyonu $(I - \frac{i}{\hbar}H \Delta t)$ olarak açabiliriz:
$$\left|\psi(t + \Delta t)\right\rangle = (I - \frac{i}{\hbar}H \Delta t) \left|\psi(t)\right\rangle$$
$$\left|\psi(t + \Delta t)\right\rangle - \left|\psi(t)\right\rangle = -\frac{i}{\hbar}H \Delta t \left|\psi(t)\right\rangle$$
$$\frac{\left|\psi(t + \Delta t)\right\rangle - \left|\psi(t)\right\rangle}{\Delta t} = -\frac{i}{\hbar}H \left|\psi(t)\right\rangle$$
Son denklem $\Delta t$ değeri $0$'a giderken doğrudan $t$'ye göre türeve eşit olacak, dolayısıyla şöyle bir denklem elde ediyoruz:
$$\boxed{\frac{\partial \left|\psi\right\rangle }{\partial t} = -\frac{i}{\hbar}H \left|\psi\right\rangle}$$

İşte meşhur (zamana bağlı) Schrödinger Denklemi! Klasik olarak kitaplarda ilk olarak bu denklemden başlanıp zamanda evrim operatorü $U$'nun eksponensiyel formda olduğu gösterilir. Fakat buradaki yolla ilerlediğimizde olayın mantığının ters, yani aslında Schrödinger Denklemi'nin mevcut formu zamanda evrim operatörünün eksponansiyel formda olmasından kaynaklandığı görülüyor. Aynı zamanda denklemin 'gökten inmediğini' de göstermiş olduk! Bu çözümlemeyi üstad Suskind'in 'Advanced Quantum Mechanics' serisinin ilk dersinde gördüm ve gördüğüm gibi yaşadığım 'aydınlanma' ile bunun tam bir 'Aslında Fizik' vakası olduğuna karar verdim. Umarım ilgilenen kişiler için de benzer etkiyi yaratır diyelim!

11 Aralık 2017 Pazartesi

Spin Sistemleri

Temel parçacıkların bir manyetik alan altında yön değiştirdikleri gözleminden yola çıkarak tıpkı birer küçük mıknatıs gibi davrandıkları söylenebilir. Bu özellik, belirli bir yüke sahip makroskopik cisimlerin dönüşleriyle ilişkili olarak kazandıkları manyetik özelliklerden yola çıkarak 'spin' olarak adlandırılmış fakat temel parçacıklar söz konusu olduğunda asıl olarak kuantum mekaniğinin temel çıkarımlarıyla ilişkili. Spin sistemleri istatistiksel fiziğin giriş aşamasında temel istatistiksel ilişkileri (ortalama/standart sapma/dağılımlar) ve zayıf etkileşimleri incelemek üzere oldukça kullanışlı 'oyuncak modeller'. Sahip oldukları manyetik moment ($\mu_0$) ile ilişkili olarak bulundukları noktadaki manyetik alanla ve tanımlanabilecek daha karmaşık yollarla etrafındaki spinlerle etkileşebiliyorlar. Bu örnek çalışmada spin $1/2$ olan parçacıkların davranışlarını, yani sadece 'yukarı' ve 'aşağı' şeklinde iki durumda yönlenebilen, dolayısıyla 2 farklı duruma sahip olan sistemleri inceleyeceğiz (örneğin elektron). (Spin'in 'aşağı' olması ortamdaki manyetik alanın tersi yönde, yukarı olması ise manyetik alanla ynı yönde yönlendiği anlamına geliyor.)

Bu uygulamada spin sistemlerinin rastgele dağılımlarından başlayarak ortalama manyetik momentlerini, bunların standart sapmalarını inceleyip ardından sistemin dışarda bir $T$ sıcaklığında bir ısı rezervi ile etkileşimi halinde nasıl dağılacağını modelleyip bu dağılımın Boltzmann dağılımı olduğunu göstereceğiz.

Çalışma için Python'un numpy, scipy ve matplotlib kütüphanelerini ağırlıklı olarak kullanacağız.

Öncelikle kullanacağımız kütüphaneleri yükleyerek başlayalım.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stat
import random, math, os, pylab
from matplotlib.patches import FancyArrowPatch

İlk örneğimizde $p$ olasılıkla yukarı, $1-p$ olasılıkla da aşağı yönlü $n$ tane spin oluşturalım. Spinlerimiz iki değer aldıklarından ve yalnızca bu değerlerin olasılıkları ile modellendiklerinden birer Bernoulli rassal değişkeni (random variable) şeklinde dağılacaklar. scipy.stats.bernoulli.rvs fonksiyonu bize $p$ olasılıklı $n$'tane $0$ ve $1$'ler üretiyor. Daha kolay çalışabilmek için (örneğin toplam spin sayısını doğrudan toplayarak bulmak için) bunları $1$ ve $-1$'e dönüştürmek için rassal değişkeni 2 ile çarpıp 1 çıkarıyoruz.

In [2]:
n = 36
p = 0.5
spin = stat.bernoulli.rvs(p, size=n)*2-1
print(spin)
print('Net spin = {} '.format(spin.sum()))
spin = list(spin)
[ 1  1 -1 -1 -1  1  1 -1  1 -1  1  1  1  1  1 -1  1  1  1  1  1 -1  1 -1 -1
 -1 -1  1  1  1 -1 -1  1  1 -1 -1]
Net spin = 6 

Oluşturduğumuz 36 spin'i $L \times L$ bir lattice'e yerleştirelim. Bunun için sayfanın sonunda ek olarak tanımladığımız görselleştirme fonksiyonumuz plot_configurations'u kullanacağız. Yeşil ile gösterilenler 'aşağı' yönlü spinlerken, kırmızı ile gösterilenler 'yukarı' yönlü spinler.

In [5]:
L = 6
plot_configurations([spin], L)

Üstteki konfigurasyon spin olasılıklarının eşit, yani yukarı gelme olasılığının da, aşağı gelme olasılığının da $0.5$ olduğu durumu gösteriyor. Aynı şeyi $p$ değerini değiştirip farklı olasılıklar için çizelim. ($p$ değerini spin'in 'yukarı' yönlü gelme ihtimali olarak düşünebilirsiniz.). Aşağıda $p = 0.1, 0.5, 0,8$ değerleri için konfigürasyonlar çizilmiş.

In [39]:
n = 16
L = 4
prob = np.arange(0.1,1.0,0.1)
for p in prob:
    spin = stat.bernoulli.rvs(p, size=n)*2-1
    if(p  == 0.1 or p == 0.5 or p == 0.8 or p == 1):
        print('Olasılık (p) = ', p)
        plot_configurations([spin], L)
Olasılık (p) =  0.1
Olasılık (p) =  0.5
<matplotlib.figure.Figure at 0x7fa2714f5390>
Olasılık (p) =  0.8
<matplotlib.figure.Figure at 0x7fa270f34c50>

Şimdi de net manyetik spin değerinin ortalamasını bulmaya çalışalım. Bunun için birbirinin 'benzeri', her birinde 36 spin yer alan 100 tane kopya sistem (ensemble) oluşturacağız. Sonrasında her birinin net spin değeri üzerinden ortalama ve ortalamadan sapmayı yani variance'ı hesaplayacağız.

In [7]:
n = 36
p = 0.25

N_trial = 100
ensemble = []
for i in range(N_trial):
    spin = stat.bernoulli.rvs(p, size=n)*2-1
    ensemble.append(spin)
net_spins = [np.sum(k) for k in ensemble]
print('Mean spin for prob = {} = {}'.format(p, np.mean(net_spins)))
print('Variance of spin for prob = {} = {}'.format(p, np.std(net_spins)/np.sqrt(n)))
Mean spin for prob = 0.25 = -17.78
Variance of spin for prob = 0.25 = 0.8457672387969531

Yukarıda yaptığımızı bir de farklı sayıda spin'le deneyelim. Aşağıda öncelikle number değişkeni ile 10'dan 100'e kadar 2'şer artan farklı spin sayıları oluşturup, ardından her spin sayısı değeri için yukarıdaki gibi 100'er örnekten oluşan ensemble'lar elde edip, her bir spin değeri için ortalamaları (u_t) ve göreli hataları (Std_t) kaydediyoruz. Göreli hata (relative error) $\sigma^2/N$ şekilde tanımlanıyor; $\sigma^2$ yani varyans $N$ ile gittiğinden, göreli hata $1/\sqrt{N}$ ile gidiyor.

In [17]:
number = np.arange(10,100,2)
u_t = []
Std_t = []
for k in number:
    p = 0.5
    N_trial = 100
    ensemble = []
    for i in range(N_trial):
        spin = stat.bernoulli.rvs(p, size=k)*2-1
        ensemble.append(spin)
    net_spins = [np.sum(m) for m in ensemble]
    u_t.append(np.mean(net_spins))
    Std_t.append(np.std(net_spins)/k)

The avarage value is $u = N \times (2p-1)$

In [19]:
plt.plot(number, u_t)
plt.title('Net spin ($p = 0.5$)')
plt.xlabel('Number of spins')
plt.show()

Ortalama değer, $p = 0.5$ olduğundan $0$ civarında dolanıyor fakat spin sayısı arttıkça varyansın arttığı da açıkça görülebiliyor. Farklı bir $p$ değeri için, örneğin $p = 0.3$ için deneyelim; sonuç şöyle olacak:

In [20]:
number = np.arange(10,100,2)
u_t = []
Std_t = []
for k in number:
    p = 0.3
    N_trial = 100
    ensemble = []
    for i in range(N_trial):
        spin = stat.bernoulli.rvs(p, size=k)*2-1
        ensemble.append(spin)
    net_spins = [np.sum(m) for m in ensemble]
    u_t.append(np.mean(net_spins))
    Std_t.append(np.std(net_spins)/k)
plt.plot(number, u_t)
plt.title('Net spin ($p = 0.3$)')
plt.show()

Göreli hata (relative error) ise $1/\sqrt{N}$ ile gittiğinden her iki durumda da azalıyor olacak.

In [22]:
plt.plot(number, Std_t)
plt.title('Relative Error')
plt.xlabel('Number of spins')
plt.show()

Isı Rezervi (Heat Bath) ve Boltzman Dağılımı

  • $T$ sıcaklığında bir ısı rezervi ile etkileşimde olan tek bir spin için 'aşağı' (down) ve 'yukarı' (up) durumunda olma ihtimalleri şu şekilde veriliyor ($\beta = \frac{1}{kT}$, $k$: Boltzmann sabiti): $$ P(up) = \frac{e^{-\beta E_{up}}}{Z}$$ $$P(down) = \frac{e^{-\beta E_{down}}}{Z}$$ $Z$ fonksiyonu 'bölüşüm fonksiyonu' (partition function) olarak adlandırılıyor ve olasılıkları normalize ediyor: $$ Z = \sum{e^{-\beta E_{up}} + e^{-\beta E_{down}}}$$
  • $\pm \mu_0$ manyetik spinine sahip iki spin durumunun $B$ manyetik alanında kazandıkları enerjiler: $E_{up} = - \mu_0 B$ and $E_{down} = \mu_0 B$
  • Sistemin ortalama enerjisini şu şekilde hesaplıyoruz: $$\bar{E} = -\frac{\partial \ln{Z}}{\partial \beta} = -\mu_0 B \tanh{\beta \mu_0 B}$$
  • Verilen $N$ spin için deneklemelerimiz $N$ ile ölçekleniyor; örneğin ortalama enerji şuna dönüşüyor $$\bar{E} = -N\mu_0 B \tanh{\beta \mu_0 B}$$
  • Limit sıcaklık değerleri için enerjinin limit değerleri: $$ \beta \rightarrow 0 \hspace{1cm} (T \rightarrow \infty), \hspace{1cm} \bar{E} \rightarrow -N \mu_0 B $$ $$ \beta \rightarrow \infty \hspace{1cm} (T \rightarrow 0), \hspace{1cm} \bar{E} \rightarrow 0 $$

Yukarıdaki değerleri simulasyon yaparak kestirmeye çalışalım. İlk olarak tek bir spin ile başlayalım. Kolaylık olsun diye enerjiyi $\mu_0B$ biriminde ölçelim (yani $\mu_0B = 1$) ve Boltzmann sabitini de $1$ olarak alalım. Verilen bir temp değerinde bir spin oluşturalım. Örneğin $T = 1.0$ için örnek bir spin üretelim:

In [40]:
temp = 1.0
beta = 1/temp
Z = np.exp(beta) + np.exp(-beta)
prob_up = np.exp(beta)/Z
prob_down = np.exp(-beta)/Z

a = np.random.uniform()
if a < prob_up:
    spin = +1
else:
    spin = -1

print('Spin is = ',spin)
Spin is =  1

Şimdi de 100 tane spin'den oluşan bir 'ensemble' oluşturup ortalama spin değeri ve varyansını hesaplayalım.

In [41]:
temp = 1.0
beta = 1/temp
Z = np.exp(beta) + np.exp(-beta)
prob_up = np.exp(beta)/Z
prob_down = np.exp(-beta)/Z
spin = []
for k in range(100):
    a = np.random.uniform()
    if a < prob_up:
        spin.append(+1)
    else:
        spin.append(-1)
print('Ortalama spin = ', np.mean(spin))
print('Varyans = ', np.var(spin))
Ortalama spin =  0.7
Varyans =  0.51

Bunu farklı sıcaklık değerleri için yapıp, ortalama enerjinin sıcaklık $T$ ile nasıl değiştiğine bakalım. Oluşturduğumuz random ensemble sonuçlarını teorik değer olan $\bar{E} = -N\mu_0 B \tanh{\beta \mu_0 B}$ ile karşılaştıralım ($\mu_0B = 1$).

In [42]:
temp_range = np.arange(0.01, 25.1, 0.1) # sıcaklık değerleri
beta = 1/temp_range
Z = np.exp(beta) + np.exp(-beta) #bölüşüm fonksiyonu
prob_up = np.exp(beta)/Z #'yukarı' olma ihtimali
prob_down = np.exp(-beta)/Z #'aşağı' olma ihtimali
diff_t_samples = []

for i in range(len(temp_range)):
    ensemble = []
    for k in range(100):
        a = np.random.uniform()
        if a < prob_up[i]:
            ensemble.append(+1)
        else:
            ensemble.append(-1)
    diff_t_samples.append(ensemble)
    
net_spins = [-np.sum(x) for x in diff_t_samples]

plt.plot(temp_range, net_spins,'b',label = 'data')
plt.plot(temp_range, -100*np.tanh(beta), 'r', label='N x tanh(beta)')
plt.xlabel('T - sıcaklık')
plt.ylabel('Net magnetic enerji')
plt.grid(True)
plt.legend()
plt.show()    

Benzer şekilde tek bir spin yerine N tane spin alıp, $T$ sıcaklığında bir ısı rezervi ile etkileştirdiğimizde, spin konfigürasyonlarının sıcaklıkla değişimi aşağıdaki gibi olacak. Öncelikle düşük sıcaklıkta tüm spinler manyetik alan yönünde yönlenme eğiliminde olacaklar yani 'yukarı' (up) durumda olacaklar. Sıcaklık arttıkça ortamla etkileşim exponensiyel olarak azalacağından 'aşağı' yönlü spinler de ortaya çıkmaya başlayacak ve sıcaklığı daha da arttırdığımızda ortalama yarısının 'yukarı', yarısının 'aşağı' olduğu maksimum düzensiz duruma doğru yaklaşılacak. Bu durumlar, yukarıda verilen limit değerler çerçevesinde tutarlı görünüyor.

In [37]:
temp_range = np.arange(0.01, 25.1, 0.1) # sıcaklık değerleri
beta = 1/temp_range
Z = np.exp(beta) + np.exp(-beta)
prob_up = np.exp(beta)/Z
prob_down = np.exp(-beta)/Z
diff_t_samples = []

for i in range(len(temp_range)):
    ensemble = []
    for k in range(100):
        a = np.random.uniform()
        if a < prob_up[i]:
            ensemble.append(+1)
        else:
            ensemble.append(-1)
    diff_t_samples.append(ensemble)
    if(i % 30 == 0): #Her üç derecede bir çiz
        print('Spin Dağılımı T= ', temp_range[i])
        plot_configurations([ensemble], 10)
        
    
net_spins = [-np.sum(x) for x in diff_t_samples]

plt.plot(temp_range, net_spins,'b',label = 'data')
plt.plot(temp_range,-100*np.tanh(beta), 'r', label='N x tanh(beta)')
plt.xlabel('Sıcalık')
plt.ylabel('Toplam net magnetik spin')
plt.grid(True)
plt.legend()
plt.show()   
Spin Dağılımı T=  0.01
Spin Dağılımı T=  3.01
<matplotlib.figure.Figure at 0x7fa270d3c9e8>
Spin Dağılımı T=  6.01
<matplotlib.figure.Figure at 0x7fa271599a20>
Spin Dağılımı T=  9.01
<matplotlib.figure.Figure at 0x7fa271629128>
Spin Dağılımı T=  12.01
<matplotlib.figure.Figure at 0x7fa2715c8080>
Spin Dağılımı T=  15.01
<matplotlib.figure.Figure at 0x7fa2715b0b00>
Spin Dağılımı T=  18.01
<matplotlib.figure.Figure at 0x7fa27125c588>
Spin Dağılımı T=  21.01
<matplotlib.figure.Figure at 0x7fa270fca048>
Spin Dağılımı T=  24.01
<matplotlib.figure.Figure at 0x7fa271426ba8>

Bir sonraki aşama spinleri sadece ısı rezervi ile değil, birbirleriyle etkileştirmek. Bunun için de istatistiksel fiziğin en temel modellerinden biri olan Ising Modeli'ni inceliyor olacağız.

In [43]:
##Spin konfigürasyonlarını görselleştirmek için kullanılan fonksiyon
def plot_configurations(conf, L):
    pylab.figure(figsize=(6 * len(conf), 6))
    s = 1.0 / L
    for i_c in range(len(conf)):
        c = conf[i_c]
        colors = {}
        for i in range(L*L):
            if(c[i] == 1):
                colors[str(i)] = 'r'
            else:
                colors[str(i)] = 'g'
        #print(colors)
        for l in range(L ** 2):
            x, y = ((l // L) + 0.5) * s, ((l - (l // L) * L) + 0.5) * s
            dy = c[l] * 0.85 / float(L)
            arrow = FancyArrowPatch((x, y - 0.5 * dy), (x, y + 0.5 * dy), \
                    fc= colors[str(l)], color='.2', lw=0, alpha=.8, \
                    arrowstyle="Simple, head_length=" + str(0.9 * 150 * s) \
                    + ", head_width=" + str(0.9 * 150 * s) + ", tail_width=" \
                    + str(0.9 * 40 * s))
            pylab.gca().add_patch(arrow)
        pylab.axis('scaled')
        pylab.axis([0, 1, 0, 1])
        pylab.gca().set_xticks([])
        pylab.gca().set_yticks([])
        [pylab.axhline(y=(i * s), ls='--', c='.2', lw=0.5) for i in range(L)]
        [pylab.axvline(x=(j * s), ls='--', c='.2', lw=0.5) for j in range(L)]
    
    pylab.tight_layout()
    #pylab.savefig(output_dir + '/' + filename)
    pylab.show()
    pylab.clf()
    
#Kod referansı: https://www.coursera.org/learn/statistical-mechanics 
#coursera statistical mechanics algorithms and computation - werner krauth