Monte Carlo

 1 import random
 2 import matplotlib.pyplot as plt
 3 import numpy as np
 4 
 5  
 6 
 7 
 8 # 函式模擬點的隨機掉落，並分為兩組
 9 def check(tmp):
10     global ans, ansX, ansY
11     global inX, inY
12     global outX, outY
13 
14     tmpX = random.random()
15     tmpY = random.random()
 
16     if (tmpX ** 2 + tmpY ** 2) <= 1:
17         ans = ans + 1
18         inX = np.append(inX, tmpX)
19         inY = np.append(inY, tmpY)
20     else:
21         outX = np.append(outX, tmpX)
22         outY = np.append(outY, tmpY)
23     ansX = np.append(ansX, tmp)
24     ansY = np.append(ansY, float(ans / tmp * 4))
 
25 
26 
27 # 變數宣告
28 ans = 0;
29 inX = np.array([])
30 inY = np.array([])
31 outX = np.array([])
32 outY = np.array([])
33 ansX = np.array([])
34 ansY = np.array([])
35 circleX = np.linspace(0, 1, 10000)
36 circleY = (1 - circleX ** 2) ** (0.5)
37 
38 # 主體
39 N = int(input("迴圈次數: "))
40 # N = 1000
 

41 for i in range(N):
42     check(i + 1)
43 
44 fig = plt.figure()
45 p1 = fig.add_subplot(121)
46 p2 = fig.add_subplot(122)
47 # i行j列，一維順序下的第k個
48 
49 p1.axis("square")
50 p1.axis([0.0, 1.0, 0.0, 1.0])
51 p1.scatter(inX, inY, c="r", marker=".")
52 p1.scatter(outX, outY, c="b", marker=".")
53 p1.plot(circleX, circleY, 'r')
54 
55 p2.plot(ansX, ansY)
56 
57 # 結果輸出
58 plt.show()
59 print(format(float(ans / N * 4), "0.6f"))