MONTE CARLO YÖNTEMLERİ VE KISMİ DİFFERANSİYEL DENKLEM ÇÖZÜMÜNE UYGULANMALARI

Monte Carlo Yöntemleri ve Kısmî Differansiyel Denklem Çözümüne Uygulanmaları

Kaya YAZGAN TCDD

ÖZET

Bu yazıda Monte Carlo yöntemleri ve rasgele sayı üretimi hakkında temel bilgiler verilmiştir. Eliptik bir kısmî differansiyel denklem veren, kare şeklindeki bir levha üzerinde ısı dağılımı probleminin çözümünde kullanılabilecek iki Monte Carlo yöntemi tanıtılmıştır.

SUMMARY An ıntroductory Information on-Monte Carlo methods and random number generation. is 'presented in this article. A problem of steady-state heat distribution on a square plate which gives a partial differential equation of elliptic type is approached by two different Monte Carlo methods.

1. GÎRİŞ tik bakışta farklı olan birçok fiziksel olayın matematiksel olarak ifade edildiklerinde eşdeğer denklemler verdikleri uzun süredir bilinen bir gerçek. Yüzyılımızda bir yandan sosyal bilimlerin gelişimi ile istatistik de bilimsel bir nitelik kazanmış, diğer yandan başta atom fiziği olmak üzere bazı doğa bilimi alanlarında istatistiksel yöntemler uygulanmaya başlamıştır. Kökü eskilere dayanmasına karşın Monte Carlo yöntemlerinin geniş uygulama alanı bulması da bu gelişim içindedir.

Monte Carlo yöntemlerinin kesin bir tanımını vermek yerine, ilişkili diğer kavramlarla birlikte Monte Carlo'nun da ne olduğunu belirtmeye çalışalım. Bir fiziksel olayın matematiksel olarak formüle edilmesi «uygulamalı matematik» alanına girer. İstatistiksel bir olaym «matematiksel ifadesi» matematiksel istatistik olur. İstatistiksel bir olaym bir başka istatistiksel olaya dönüşümü «benzeşim» (simulasyon) yöntemidir. İstatistiksel bir olaym kendi içinden alınan örneklerle incelenmesi «örnekleme» (sampling) diye adlandırılır. «Monte Carlo» ise önce matematiksel bir sorunun istatistiksel yoldan çözülmesi olarak görüldü, örneğin 1777 de Buffon, bir çubuk atma deneyinin istatistiksel sonucundan yararlanarak ^ sayışma yakınsayan bir deney önerdi. Yöntemin sistemleştirilip ünlü kumar şehrinden esinlenen adını alması ise 1949 da Metropolis ve Ulam'ın yazısı iledir [1].

Monte Carlo yöntemlerinin bugünkü uygulanış alanları ve şekilleri üç gurupta toplanabilir.

1. istatistiksel bir olaym önce matematiksel bir yapıya dönüştürülmesi, sonra da bu matematiksel olayın bir başka istatistiksel olay aracılığıyla incelenmesi,

2. İstatistiksel özellikleri oldukça belirgin fiziksel bir olayın doğrudan doğruya istatistiksel bir yöntemle çözülmesi, örneğin atom fiziğinde parçacıkların hareketleri ile ilgili problemlerin çözümü,

3. İlk bakışta istatistiksel özelliği göze çarpmayan fiziksel bir olayın önce matematiksel olarak ifade edilmesi, sonra da bu matematiksel yapının istatistiksel bir yöntemle çözülmesi.

Bu incelemede özellikle üçüncü tip uygulamalar üzerinde durulmakta, örnek olarak ikinci aşamadan kısmî differansiyel denklemlerin çözümü ele alınmaktadır.

İkinci aşamadan kısmî differansiyel denklemler, örneğin Laplace, Poisson, ısı iletimi gibi denklemler, fizik ve mühendislik dallarında pek çok karşılaşılan denklemlerdir. Bunların analitik çözümleri ancak belirli özel durumlarda, oldukça uzun ve çapraşık yöntemler izlenerek elde edilebilmektedir. Bilgisayarların gelişimine paralel olarak gelişen sayısal yöntemler ise, uygulama kolaylığı, hızı ve analitik çözümü olmayan problemleri bile çözülebilir kılması ile günümüzde bu alanda kesin bir egemenlik sağlamıştır.

Sayısal yöntemler kendi içinde iki gurupta incelenebilir. «Deterministik» sayısal yöntemlerde kısmî differansiyel denklem bir sonlu farklar denklemine dönüştürülür, ilerde değinileceği gi-

bi, katsayıların oluşturduğu matrisin tersi bulunmaya çalışılır (*) [2]. Elde edilen sonlu farklar denklemleri «olasılıksal» yöntemler, rasgele sayılar dizisi ve istatistiksel değerlendirmeler kullanılarak çözülürse asıl konumuz olan Monte Carlo yöntemleri alanına girilmiş olur. önce yukarda adı geçen rasgele sayılar dizisi üzerinde duralım.

2. RASGELE SAYILAR DİZİSİ ÜRETİMİ

Birbiriyle hiçbir ilişkisi olmayan sayılardan oluşan bir dizi elde etmek yüzyılımızın başından beri üzerinde çalışılan bir konudur [3], [4]. 1927 yılında Tippett tarafından yayınlanan 40000 sayılık diziyi, 1939 da 100000, 1955 de 1000000 sayılık dizilerin yayınlanması izledi. Günümüzde önceden basılmış sayı dizilerinin kullanılmasından çok algorıtmik bir yolla rasgele sayıların teker teker elde edilmesi ve hesaplamalarda kullanılması yoluna gidilmekte.

Algoritmik bir yöntemle rasgele sayılar elde etme fikri ilk bakışta olanaksız görünüyorsa da geliştirilen birçok yöntemle «rasgeleymiş gibi davranan» sayılardan oluşan dizilerin elde edilmesi günümüz bilgisayarlannda çok görülen bir uygulamadır. Yazımızda rasgele sayı dizisi denildiğinde böyle bir dizi amaçlanmaktadır.

John von Neumann tarafından 1946 da ortaya atılan «orta kareler yöntemi» birbirinden oldukça bağımsız sayılar dizisi veren bir yöntemdir. 2n basamaklı bir sayı, serinin ilk sayısı olarak seçilir, bu sayının orta n basamağının karesi olan sayı serinin ikinci sayısıdır. Her sayı bir sonrakinin bulunmasında kulanılarak seri sürdürülür, örneğin 2425 sayısıyla başladığımızı düşünelim. 42 nin karesi olan 1764 ikinci sayımızdır, 76 nm karesi olan 5776 üçüncü sayımızdır... Bu yöntemin sakıncası kısa süre sonra kendini yineleyen serilere dönüşmesi, orta basamaklarda sıfır elde edince ise hep sıfır vermesidir, örneğin yukarda başlanan dizi 16 tane rasgele sayı verdikten sonra sıfıra ulaşır (•*). Daha iyi bir başlangıç sayısı seçimi ile çok daha uzun diziler elde edilebilir, örneğin Metropolis 38 bitlik 75000 sayılık bir dizi elde etti. Çok kullanılan diğer bir yöntemde rasgele sayı, kendinden bir önceki sayıya dayanarak ve elden geldiği kadar büyük bir periyotla kendini yinelemesi için belirli koşullara uyularak hesaplanır. X rasgele sayıyı, n ise adım sayısını gösterirse, Xn+1;

Xn+1 = (aXn + c) mod m (*) Bu genellemenin dışında kalan bir istisna, hiperbolik denklemlere uygulana*^ Karakteristik Eğriler Yöntemidir.

(*«)2425, 1764, 5776, 5929, 8464, 2116, 0121, 0144, 0196, 0361 1296, 0841, 7056, 0025, 0004, 0000. 0000,

formülünden hesaplanabilir. Burada,

Xo a

0 serinin ilk sayısı 2 çarpan

c ^ 0 artım

m > Xo, a, c modülusdur. Mümkün olan en büyük periyotlu diziyi elde etmek için bu sabit sayılar birkaç teorem yardımıyla belirlenir [5], [6], [7].

Bu iki yöntemin dışmda, daha önce elde edilen

sayıların belirli katsayılarla çarpımları toplanarak, çembersel kaymalar ve VEYA mantık işlemleri kullanarak algoritmik rasgele sayı üretim yöntemleri de vardır.

-,-

3. MONTE CARLO YÖNTEMLERİYLE İNTEGRAL ALMAK

Rasgele sayılar dizisinin nasıl elde edildiğini gördükten sonra,

I = f f (x) dx

(D

o

eşitliğindeki I yi iki Monte Carlo yöntemi ile bulmaya çalışalım (Şekil 1).

f(x)

M

0 Şekil 1.

X

a. tA aradığımız I değerine, Şekil 1 deki taralı alana ulaşan bir ortalama olsun. (0,1) aralığı içinde düzgün dağılmış bağımsız rasgele sayılardan i'inci deneyde elde edilenleri u; ve Vj ile gösterelim.

tAı = M, eğer (Uj.v,) noktası Şekil 1 in taralı kısmında ise, f (u,) ^ v, M t ;Ai = 0, eğer f (u,) < v, M •••

N deney sonunda elde edilen sonuçtur. B beklenen değeri gösterirsek,

B | t. != I yazılabilir.

b. iki ayn rasgele sayı kullanmak yerine daha basit bir yöntemle de (1) integralinin sonucuna yaklaşabiliriz.

t« = Yine,

2 f(Ul)

1=1

tB J = îlk bakışta umulanm tersine bu iki yöntemden daha basit olanı, b yöntemi, diğerine göre daha küçük bir varyans vermektedir, sonuca daha büyük bir hızla yaklaşmak mümkündür [8].

Monte Carlo yöntemleri hakkında bir fikir edindikten sonra ikinci aşamadan eliptik bir kısmî differansiyel denklemin kesikli ve sürekli Monte Carlo yöntemleriyle çözülmesini inceleyelim.

Problem Bir kenarı 100 °C da, diğer üç yanı 0°C da tutulan kare şeklinde bir levhadaki durgun durum ısı dağılımını inceleyelim (Şekil 2).

t=0

Şekil 2.

«t» ile gösterilen ısının Laplace denklemine göre dağılmış olduğu bilinir:

3* t

=o

(2)

Karenin kenarları a uzunlukta ise, denklemle birlikte aşağıdaki smır şartlan problemin tanımlanmasını tamamlar :

t(x,y) =0 x = a, 0 < y < a y = 0, 0 < x < a

y = a, 0 < x < a için 't (x, y) = 100 x = a, 0 < y < a için.

Bu problemin analitik çözümü bilinmektedir [9] ve bu nedenle sayısal çözümlerimizdeki yanılmayı bulmakta kullanılabilir.

4. KESİKLİ MONTE CARLO

Bu yöntemde ısının yayılması bir ızgara üzerindeki rasgele yürüyüş yardımıyla elde edilir (Şekil 3). Deterministik sayısal yöntemlerde olduğu gibi, verilen (2) denklemi sonlu farklar

J+l

j

1>

V(t E)

Şekil 3.

denklemine çevrilir. Eğer tjj , (ij) noktasmdaki ısıyı gösterirse ve üçüncü aşamadan daha yukarı türevli terimler ihmal edilirse aşağıdaki iki Taylor seri açılışı elde edilir.

t

= ti., - A* 3t .

3x2

(AxP

yt

3!

(3)

s tlı(

( A * ) 2 a21

d31 3!

(4)

(3) ve (4) denklemlerinin taraf tarafa toplamı (32t/9x2) için bir yaklaşım verir:

t,.,., + t ı + I ı l s 2 tM + (A*)2

a21 3x2

*,_,., — 2t, LJ + 1 )

(5)

(A*)*

Benzer işlem (9zt/3y2) için de yapılırsa,

(6)

3y2

(Ay)2

elde edilir.

Laplace denklemi (2), içine (5) ve (6) da elde edilen sonlu fark yaklaşımları yerleştirilirse ve Ax = Ay = h alınırsa,

M=

(7)

bulunur.

Deterministik sayısal yöntemlerle Laplace denkleminin çözümü, (7) denkleminin her nokta için yazılması ile elde edilen sonlu farklar denklemleri sisteminin çözümüne dönüşür.

Şimdi de olasılıksal bir yönde problemi ele alalım. Isının levha üstünde yayılmasını, Şekil 3 deki noktalara uğrayarak yapılan bir rasgele yürüyüşe benzetelim. P (i,j), (i,j) noktasında olmanın olasılığını göstersin. Yürüyüşü yapan cisim bir adımda küçük karelerin ancak bir kenarı boyunca ilerliyorsa ve bir noktadayken eşit olasılıklarla komşu dört noktadan birine yöneliyorsa (i,j) noktasında olmanın olasılığı, bir adım önce komşu dört noktada olma olasılıklarının toplamına eşittir.

(i+lj)

kış yapıldığını gösteren sayıyı) bu noktanın sınır değeri (bizim problemimizde 0 veya 100) ile çarpın. Bu çarpımları toplayın ve Y sayısının içindeki sayıya (toplam deney sayısına) bölün. Büyük deney sayılan için bu bölüm No başlangıç noktasındaki ısıya yaklaşır (*).

Dikkat edilirse yukardaki algoritma ile levhanın çevresinde çeşitli sınır noktalarında uygulanan ismin No başlangıç noktasına olan etkileri, bu noktaya yakmlıklan ve sınır değerinin büyüklüğü ile artmaktadır. Bu da fiziksel olayın temel niteliğidir.

5. SÜREKLİ MONTE CARLO

Sürekli Monte Carlo yönteminde ise, Şekil 3 de görülen ızgara üzerindeki kesikli hareket yerine Şekil 4 de görüldüğü gibi bir başlangıç noktasından başlayıp bir çember çevresi üzerindeki noktalara eşit olasılıkla gidilir. Çemberler çözüm alanı içinde kalmak koşuluyla her büyüklükte olabilir, hattâ çemberden başka şekiller aracılığıyla da rasgele yürüyüş yapılabilir [10].

P (i,j-l

(8)

(7) ve (8) denklemlerinin aynı fiziksel olayı tanımladıkları açıktır.

Sözü edilen rasgele yürüyüşü aşağıdaki algoritmaya göre bilgisayar için programlamak çok kolaydır.

4.1. Kesikli Monte Carlo Algoritması

ADIM 1. Yolculuk sayılarını saymak için kullanılacak Y sayıcısını sıfırlayın.

ADIM 2. Koordinat parametrelerini ( I J ) , başlangıç noktasının (No) koordinatlarına eşitleyin.

ADIM 3. 0 ile 1 arasında dağılmış bağımsız rasgele sayı dizisinden bir R sayısı alnı. R nin değerine göre eşit olasılıklarla komşu noktalardan birine gidin, (örneğin 0<R<0,25 ise sola ilerleyin, yani I parametresini bir arttınn).

Eğer yeni ulaşılan nokta bir sınır noktası ise, ADIM 4 e aksi halde ADIM 3 e gidin.

ADIM 4. Ulaşılan sımr noktasına ait sayıcıyı bir arttınn. Yolculuk sayısını veren Y sayıcısını bir arttınn. Y nin değeri önceden belirlenen bir üst limite ulaştıysa ADIM 5 e, aksi halde ADIM 2 ye gidin.

ADIM 5. Her sınır noktası için ayrılan sayıcının içindeki sayıyı (bu noktadan toplamda kaç çı-

Şekil 4.

5.1. Sürekli Monte Carlo Algoritması

ADIM 1. PAY ve PAYDA adlı iki sayıcıyı Gifırlayın.

ADIM 2. Koordinat parametrelerini (X, Y) başlangıç noktası Nom koordinatlarına eşitleyin. ADIM 3. (X, Y) noktasını merkez alan, çözüm alanının içinde kalan en büyük çemberi çizin.

(*) Burada sunulan algoritma yöntemi kolay açıklamak amacıyla verilmiştir. Bilgisayar için daha verimli bir şekilde uygulanabilir, örneğin her sınır noktası için bellekte bir yer ayrılmasına gerek yoktur.

•BunüHiJçfevre'siiıdev"e'şit öfâsıliklaFİânbakğele 'bir iîSlcta 'seçin v^'^Y^f^âfâm-etreîerifîi'Jbü? yeni noktanınkilerei eşîtleyiîîAl'ğeT > '(X,Y<J riöktası çöZürm'Jalanıriın''Şevresiîferişizilenl belirli' bir raafjih (mi),-jçin'âe ise 'AÜIM '4 e «idin/ aksi' halde ADIM 3*ü 'yineleyin.' > '.:'>.'•' ''•• . •'-'•' -ADIM 4İPAY; sayİGisına'nulâşılan sınırdaki sı•nırv değerini': ekleyin? ı PAYıDA":sayısına bir-- ekle-yin'ijEğer- PAYDA/dakİJisapl?öriceden belirlenen ıbiari'üsbdimite ulaştıysa'/itADIM ;5e; aksi'Jhalde ADIM 2yel gidin.'< nhzM: v ' j ( ' r •>><•( ADIM 5. PAYı PAYDA ya bölün, sonuç No noktasındaki ısı için aranan yaklaşımı verir.

6. SOttUÇ

•1r•/

 ' j , . - ı l'J .

•Bu.yazıda MonteiCaclo iyöm'temleri'-hakkındâ geneJL bilgi verilmişi yen basitidir fisi? probleminin

iki Monte Carlo yöntemi'iile nasıl /çözülebileceği• fgösferilmiştiF.y'Bürarâdâ pek: çok önemli konu, )amacınr"dışın'ai''taşmamak için 'ele iSİMiina-

HMŞıtırs. örneğin; dahasikarmaşık •denklemler.'verildiğinde çeşitli yönlerde ilerlemenin .LoTa'sıliğı [11]; deney birçok kez yinelendiğinde elde

edilen sonuçların dağılımının varyansının yolculuk sayısı (bîzinf"aigo^fiimlarımızdaki {Y veya ÇAYDA) ile olanvilişkislr[12],-bin yolculuk

içinf beklenen aaım^sayısı,1 içerdeki ijki 'nokta

arâsm'daki ^Hâreketnı olasılığı, sonuçtaki doğru

basamakların'"" sayısı ["1"3] gibi konular tartışılmadı.

\I

-->UTU

KesikliJVIp'rite Carloıiçin verilen algoritma bellek vöriüiidenfcgeliştİFİlelbîîecel'i^gibi, çakışmamız

sıi|asıhda daha_ç.oto^ellek kullanma, pahasma bir

yolcululej sırasmdaAîğranan noktalar heık|kındaki, bilgiler, de-saklanarak «geçmişlerin üstüste

— ,--o—

-jj

\ ıs

I' i

bindirilmesi» ycmtemi geliştirildi [14']. I

Monte ^C'arlcT^ yöntemleFİnin diğer/ sayısal yöntemlerle' kVrşılaştıraöıas'ından aşağıdaki sonuçlar-çıkmak-tadÎF." ' * ^ - ^ * ^ »»«<-- - - - ~ -J

1. Çözüm, bütün çözüm alanında değil de alanın bir veya birkaç noktasında istenirse, ancak Monte Casio. nyöa|ejalejudiğer«n'ö'ktâlardâki ç&zümü bulmaksızın istenen nokta veya noktalardaki' çözümü ve'rir'.^ «">'*• " 'YA •' ;;

2. Monte^ Carlo yöntemleri • çok. işlem gerektirmelerine ' karşılık ' determini^tik ' yontemlerdeki bûyıik'matrisİerin; saflân-masınİ' gerekfiffmemektgdi,r:-,Bu,d.a günümüzün».küçük belleklifve çok hızlı, bilgisayarları ^çinvönemli' bir özelliktir.

3. Monte Carlo yöntemlerinde "doğruluk", "çözüm üalanirun boyutlarından bağımsızdır,, halbuki, determmis,tik sayısal .yöntemlerde boyutların artması ile'doğruluk-azalar [15].

â) yaklaşık ^-sonuçların) yeterli' ı olduğu, t <mühendislıki lu.y'gulamafaıınnd'af/. küçük bellekli- nı'zk bilgisayarlar için Monte' Garlo 'yöntemleri alışılagelmiş deterministik sayısal yöntemlerden daha iyi sonuçlar verebilirler. ' '

KAYNAKLAR

1. Metrppolis N. ve Ulam S.: «The Monte Car^n'Jio-;MeÎHöd5;ıİr!rÂmer.'Statisı:''Askoc.i Völ. 44, İ1J ' 33S-341İ '1949." ''•"•'

2. Forsythe, G? E.,: «Solving Linean Algebraic Equations Can be Interesting», Bull. Amer.

'' Mathr.'JSoc!,CİVol.'57,i299-329; 1953: Oı ' '

3. . Knuth, D., E.: • «The Art of Computer Erog» ' ramming», ;Mass., Addison-Wesley, 1969.n •

4. Nance, jR. £. veCldude, O. Jr.: ,«A Bibliography on Rândom Nuniper Generatiön»,

' Comput. Rey.'.'iyoİ. 18,' 495-5q8,' 1972. "ı( : ' '

T5J Greenbergerı:'iîHotes 'on a' N'ew Pseudo-Raiı-

'•; ' dom' Nümber-i'Generâitdr»i'?JP'AlCMr(VölJb8-,

383-389, 1961.

/ " ' I . , •• ' " ' ' ' •''- •

6. Hull ve Dobell: «Random Number Gene», İlAM(. Rev;.,. VoL 4, 230-254, 1962.'.

7. Carmichael: Bull. Amer. Math. Soc, Vol.

^ } 16, 232-238, 1910.1

(

8. Morton, K. W.: «On the Treatment of Mont , 'te|

Textı Bookş», MTAÇ,

'" ' Vol. l6,"223-224, 1956.

_$.n,Garslow,,,H. S.^ejaeger, J. Ç.: Operational I Me,thods in, Applied Mathematics, 1963.

10. Muller, M. E.: «Some Continuous Monte Carlo Methods for the Dirichlet Problem», Ann. Ma«Hî rSıtâlis.,'Völ. 27j^5'69-589,^ 1956: ••

•Ui<.yEkr.lİch\:.L<:?W.:< «Monte Garlo-'Sölutioris< of Boundary Value- Prbblems vilnvol'ving- the

A D^fference Analog, fI? ı ^ ^ u ^ + K y - ^ u ^ 0», '• , Jr ACM, Vjol,,^ 204-218,495?,. îu,, ... , f ,

12., Şfarfiider, ,Yw. A. ,(ed.): The ..Monte iCaiilo n n ı Methpd, Pergamon, Press,, ,1966 ; (1962^'tani;j^rihli Rusça, baskıdan'(çeyiid-). , •• • ,

13: "Cuftİss, J; H::f«Samplingl;Metİıf6ldş Aİppli'ed •("to 'b'ifferential and ^Dîffereh'ce' E'quations»,

T' P-roe. ,ıof IBMıSeminar onJ ıScientific Computation',,;1949/, Nfevv,AYofk-.r .-.' i- ' • <

14/ ' Yazgan.;' K.: «A Survey o'rij'th'e NumerıHal

ne' 'Solutions'öf- Partial'DifferenferEquatibhs»;

''t'' Ms Thesis, ODTÜ; 1973. >'> ' ' > • • •

•

t' ' ' .

15. Kac, M.: «Applications of Statistical Methods to Differential and Integral Equa-

T tiöns»',cNotes' oh lectures rdelivere'd ât'Mî'T,

Elektrik Mühendisliği' 204