This repository has been archived by the owner on Mar 16, 2020. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 0
/
bertrand.py
186 lines (152 loc) · 6.05 KB
/
bertrand.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
from random import random
from matplotlib import pyplot as plt
import numpy as np
import math
class Bertrand:
def __init__(self):
# creates a figure and an Axes object(the plotting area)
self.fig, self.ax = plt.subplots()
# side of the equilateral triangle
self.a = 0
def distance(self, A, B):
return math.sqrt((A[0] - B[0])**2 + (A[1] - B[1])**2)
def draw_circle(self, r):
self.ax.set_xlim((-r - 1, r + 1)), self.ax.set_ylim((-r - 1, r + 1))
self.ax.set_xlabel('x'), self.ax.set_ylabel('y')
circle = plt.Circle((0, 0), r, color='black', fill=False)
# Artist: represent a standard graphical object, knows how to use the renderer to paint it
self.ax.add_artist(circle)
def draw_triangle(self, P):
x, y = P
# stack arrays in sequence horizontally (column wise)
points = np.hstack([x, y])
self.ax.plot(x, y, color='cyan', marker='o', label='P')
# calculate the other 2 points of the equilateral triangle
# rotate P by 60, 120 degrees => A, B points
angle = np.deg2rad(120)
x = P[0] * np.cos(angle) - (P[1] * np.sin(angle))
y = P[0] * np.sin(angle) + (P[1] * np.cos(angle))
A = (x, y)
points = np.vstack([points, [x, y]])
self.ax.plot(x, y, color='yellow', marker='*')
angle = np.deg2rad(240)
x = P[0] * np.cos(angle) - (P[1] * np.sin(angle))
y = P[0] * np.sin(angle) + (P[1] * np.cos(angle))
B = (x, y)
self.a = self.distance(A, B)
points = np.vstack([points, [x, y]])
self.ax.plot(x, y, color='orange', marker='*')
points = np.vstack([points, [P[0], P[1]]])
self.ax.plot(points[:, 0], points[:, 1], linewidth=3, color='blue', label='equilateral triangle')
def first_method(self):
r = float(input('r = '))
n = int(input('n = '))
self.fig.suptitle('Bertrand paradox, 1st method')
self.draw_circle(r)
# get a fixed P point on the circumference
alpha = random() * (2 * np.pi)
x, y = r * np.cos(alpha), r * np.sin(alpha)
P = (x, y)
self.draw_triangle(P)
favorable = 0
# generate random points on the circumference
for i in range(n):
alpha = random() * (2 * np.pi)
x, y = r * np.cos(alpha), r * np.sin(alpha)
# the new point
M = (x, y)
self.ax.plot(x, y, color='black', marker='.')
if self.distance(P, M) > self.a:
self.ax.plot([P[0], x], [P[1], y], color='green')
favorable += 1
else:
self.ax.plot([P[0], x], [P[1], y], color='red')
print('{0}/{1} = {2}'.format(favorable, n, favorable/n))
plt.legend()
plt.grid()
plt.show()
def second_method(self):
r = float(input('r = '))
n = int(input('n = '))
self.fig.suptitle('Bertrand paradox, 2nd method')
self.draw_circle(r)
# draw a random radius
alpha = random() * (2 * np.pi)
x, y = r * np.cos(alpha), r * np.sin(alpha)
plt.plot([0, x], [0, y], color='black', linewidth=3, label='radius')
# perpendicular angle to alpha
beta = abs(alpha - np.pi/2)
x, y = r * np.cos(alpha), r * np.sin(alpha)
# get a fixed P point on the circumference
P = (x, y)
self.draw_triangle(P)
A = r_middlex, r_middley = r / 2 * np.cos(alpha), r / 2 * np.sin(alpha)
plt.plot(r_middlex, r_middley, color='darkblue', marker='o', label='middle of the radius')
favorable = 0
for i in range(n):
# random point on the radius
m = random() * r
x, y = m * np.cos(alpha), m * np.sin(alpha)
M = (x, y)
plt.plot(x, y, color='orange', marker='.')
# perpendicular line to the radius in the generated point
x, y = M
# cord length
length = math.sqrt(pow(r, 2) - pow(self.distance(M, (0, 0)), 2))
x1, y1 = x + length * np.cos(beta), y + length * np.sin(beta)
x2, y2 = x - length * np.cos(beta), y - length * np.sin(beta)
if self.distance(A, (0, 0)) < self.distance(M, (0, 0)):
plt.plot([x1, x2], [y1, y2], color='red')
else:
plt.plot([x1, x2], [y1, y2], color='green')
favorable += 1
print('{0}/{1} = {2}'.format(favorable, n, favorable / n))
plt.legend()
plt.grid()
plt.show()
def third_method(self):
r = float(input('r = '))
n = int(input('n = '))
self.fig.suptitle('Bertrand paradox, 3rd method')
self.draw_circle(r)
# get a fixed P point on the circumference
alpha = random() * (2 * np.pi)
P = x, y = r * np.cos(alpha), r * np.sin(alpha)
self.draw_triangle(P)
# draw a smaller circle
circle = plt.Circle((0, 0), r / 2, color='grey', fill=False)
self.ax.add_artist(circle)
favorable = 0
for i in range(n):
# generate a new point within the circle
alpha = random() * (2 * np.pi)
rand_r = math.sqrt(random())
x, y = rand_r * r * np.cos(alpha), rand_r * r * np.sin(alpha)
# the new point
M = (x, y)
if self.distance(M, (0, 0)) <= r / 2:
plt.plot(x, y, color='green', marker='.')
favorable += 1
else:
plt.plot(x, y, color='red', marker='.')
print('{0}/{1} = {2}'.format(favorable, n, favorable / n))
plt.legend()
plt.grid()
plt.show()
def main():
print('choose from the 3 methods(1|2|3):\n\t1\n\t2\n\t3')
choice = 0
while 0 <= choice <= 3:
choice = int(input('choice = '))
bert = Bertrand()
if choice == 1:
bert.first_method()
break
elif choice == 2:
bert.second_method()
break
else:
bert.third_method()
break
if __name__ == '__main__':
main()