CoCalc Public Filesassignments / assignment14 / assignment14_montecarlo.htmlOpen with one click!
Authors: Gabriela Dorfman Furman, Or Golan, guy miller
assignment14_montecarlo

שמות:

In [6]:
from IPython.display import HTML
import matplotlib.pyplot as plt

<font color = blue>
חישוב שטחים בשיטת מונטה קרלו
</font>

<font color = blue>

בתרגיל זה ניישם שיטה חישובית המכונה מונטה קרלו. בעזרת שיטה זו נחשב בהמשך שטח של אליפסה. אבל בשלב ראשון ניישם את השיטה המתוארת בהמשך לחישוב π\pi.
</font>

<font color = green>
קראו בעיון את הקטע שבהמשך

In [7]:
HTML('''<iframe src="https://he.wikipedia.org/wiki/%D7%A9%D7%99%D7%98%D7%AA_%D7%9E%D7%95%D7%A0%D7%98%D7%94_%D7%A7%D7%A8%D7%9C%D7%95"
     width=900 height=700"></iframe>''')
Out[7]:

<font color = green>

בהמשך קטע קוד המיישם את האלגוריתם המתואר למעלה למציאת π\pi. הוסיפו ליד כל סימן סולמית הערה המסבירה את תפקיד משפט הקוד שמשמאלה.
</font>

In [5]:
from random import uniform #הבאת פונקציה שתביא מספר רנדומאלי לא שלם בתחום שנכתב
N=10000#מספר הנסיונות\נקודות
count=0#מספר הקליעות
for i in range(N):
    x=uniform(-1,1)#הגרלה של שיעור ה-איקס בין 1- ל1
    y=uniform(-1,1)#הגרלה של שיעור הוואי בין 1- ל1
    if x**2+y**2<=1:# בדיקה אם הוא נמצא בתוך או על המעגל
        count=count+1
float(count)/N*4
Out[5]:
3.1408

<font color = blue>

חישוב π\pi בשיטה המתוארת כאן מכונה חישוב בדגימה ישירה
</font>

<font color = green>

הריצו את הפונקציה תוך כדי הגדלת ערכי N מ-10 ועד ל- 10000. מה המסקנה? הסבירו במה במה תלוי הדיוק בחישוב π\pi?
</font>

In [7]:
from random import uniform
def pi_Evaluation(N):
    count=0
    for i in range(N):
        x=uniform(0,1)
        y=uniform(0,1)
        if x**2+y**2<=1:
            count=count+1
    return float(count)/N*4
pi_Evaluation(100000)
#הדיוק בחישוב פאי תלוי במספר ההרצות
Out[7]:
3.13716

<font color =green>

מידת הפיזור של תוצאות "ניסוי" מכונה סטיית התקן בסידרת התרגילים שבהמשך תכתבו פונקציות שמטרתן לחשב את סטיית התקן המתקבלת בחישוב π\pi בדגימה ישירה. בשלב ראשון עליכם לכתוב פונקציה המריצה עבור אותו N את הפונקציה pi_Evaluation מספר פעמים המועבר אליה כפרמטר. הפונקציה צריכה להחזיר רשימה של כל ערכי ה-π\pi שהתקבלו.

השלימו את קטע הקוד שבהמשך.

</div><

In [9]:
def get_trails(k,N):
    pi_list=[pi_Evaluation(N) for i in range(k)]
    return pi_list
get_trails(5,10000)
Out[9]:
[3.1468, 3.1424, 3.13, 3.12, 3.1432]

<font color = green>

שורת הפקודות שלמטה משרטט תרשים שכיחויות של תוצאות של הרצות שונות עבור אותו גודל של N.הריצו את הקוד עבור ערכי N שונים ובחנו כיצד משתנה פיזור התוצאות.

In [12]:
import matplotlib.pyplot as plt
N=5000
trails=10000
o=plt.hist(get_trails(trails,N),bins=30)
Out[12]:

<font color = green>

העזרו בחבילה אותה פיתחתם לחישוב פונקציות סטיסטיות כדי לענות על השאלות שבהמשך.

<font color = green>

בעזרת הפונקציה המחשבת ממוצע חשבו ערך ממוצע של π\pi המתקבל לאחר ביצוע של k ניסויים בהם בוחרים N נקודות.

In [14]:
import py_compile
py_compile.compile('my_statistics.py')
In [17]:
from my_statistics import mean
def mean_for_pi(k,N):
    return mean(get_trails(k,N))
mean_for_pi(5,10000)
Out[17]:
3.14032

<font color =green>

הקישור שלמטה הוא קישור לדף ווקיפדיה של המונח סטיית תקן. קראו על המונח . העזרו בחבילה אותה פיתחתם כדי לשרטט גרף של סטיית התקן כתלות ב- N. שנו את N מ- 10 ל-100 במרווחים של 10. קבעו את k ל-1000. השתמשו בפונקציה זו כדי לשרטט גרף של סטיית התקן כתלות ב-N. מה המסקנה?
</font>

In [2]:
from my_statistics import std
import matplotlib.pyplot as plt

get_trails(5,1000)
 10001

<font color = green>

שרטטו גרף של סטיית התקן כתלות ב: NN\sqrt{N}\over{N}. האם קיבלתם בקרוב יחס ישר? אם כן כיצד סטיית התקן תלויה ב-N?

In [0]:
 

<div dir =rtl> <font color = blue>אליפסה היא המקום הגיאומטרי של כל הנקודות שסכום מרחקיהן משתי נקודות קבועות המכונות מוקדי האליפסה הוא קבוע ושווה לקוטר הגדול של האליפסה, 2a, הקוטר הקטן של האליפסה שווה ל-2b. </font>
<font color =green> כתבו פונקציה בשם ellipse_area המקבל שלושה פרמטרים a, b ו- tol המחשבת שטח של אליפסה שחצי הקוטר הגדול שלה הוא a וחצי הקוטר הקטן הוא b בדיוק של tol. השגיאה בחישוב השטח היא מסדר גודל של סטיית התקן.</font></div>

In [45]:
import random
from math import pi
def ellipse_area(a,b,tol):
    rec_area = 2*a*2*b
    con1=0
    con2=0
    ellipse_area1=pi*a*b#חישוב לפי נוסחא
    for i in range(tol):
        x=uniform(-a,a)
        y=uniform(-b,b)
        if ((x*x)/(a*a) + (y*y)/(b*b) <= 1):
            count=count+1
    ellipse_area2=(float(count)/tol)*rec_area#חישוב לפי שיטת מונטה קרלו
    return ellipse_area1, ellipse_area2#השוואה בין השטחים
ellipse_area(6,3,10000)
Out[45]:
(56.548667764616276, 56.455200000000005)
In [0]: