-
Notifications
You must be signed in to change notification settings - Fork 1
/
110borwein
executable file
·115 lines (105 loc) · 3.29 KB
/
110borwein
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
#!/usr/bin/env python3
#
# Simple program to compute Borwein integrals, using the rectangle method,
# the trapezoidal rule and the Simpson rule,
# and print both the In and the difference between this value and π/2
# made in python 3, you can launch it with:
# ./110borwein constant_defining_the_integral
# (See subject for more details)
#
import sys
import math
import os
import random
def man_help():
print("USAGE")
print("\t\t./110borwein n\n")
print("DESCRIPTION")
print("\t\tn\tconstant defining the integral to be computed");
def valid_string(string):
i = 0
while (i < len(string)):
if (string[i] != '-' and (string[i] < '0' or string[i] > '9')):
sys.stdout.write("Invalid argument\n")
exit(84)
i = i + 1
def rectangle_method(x):
result = 0.0
b = 5000
a = 0
n = 10000
h = (b - a) / n
while (a < 10000):
result = result + get_sum(x, a * h)
a = a + 1
result = result * h
return (result)
def trapezoids_method(x):
result = 0.0
b = 5000
n = 10000
h = b / n
a = 1
while (a < 10000):
result = result + get_sum(x, a * h)
a = a + 1
result = result * 2
result = result + get_sum(x, 5000)
result = result + get_sum(x, 0)
result = result * (h / 2)
return (result)
def simpson_method(x):
res = 0.0
ult = 0.0
b = 5000
n = 10000
h = b / n
a = 1
while (a < 10000):
res = res + get_sum(x, a * h)
a = a + 1
a = 0
while (a < 10000):
ult = ult + get_sum(x, (a * h) + (h / 2))
a = a + 1
result = get_sum(x, 5000) + get_sum(x, 0) + (2 * res) + (4 * ult)
result = result * (b / (6 * n))
return (result)
def get_sum(x, i):
result = 1.0
a = 0
while (a <= x):
if (i != 0):
result = result * (math.sin(i / (2 * a + 1)) / (i / (2 * a + 1)))
a = a + 1
return (result)
if (len(sys.argv) == 2):
if (sys.argv[1] == "-h"):
man_help()
exit(0)
if (sys.argv[1].isdigit() == True):
x = float(sys.argv[1])
result = rectangle_method(x)
print(str("Rectangles:\nI" + str(int(x)) + " = " + "%.10f" % result))
if (str("%.10f" % (result - (math.pi / 2))) == "-0.0000000000"):
print(str("diff = " + "%.10f" % (-1 * (result - (math.pi / 2))) + "\n"))
else:
print(str("diff = " + "%.10f" % (result - (math.pi / 2)) + "\n"))
result = trapezoids_method(x)
print(str("Trapezoids:\nI" + str(int(x)) + " = " + "%.10f" % result))
if (str("%.10f" % (result - (math.pi / 2))) == "-0.0000000000"):
print(str("diff = " + "%.10f" % (-1 * (result - (math.pi / 2))) + "\n"))
else:
print(str("diff = " + "%.10f" % (result - (math.pi / 2)) + "\n"))
result = simpson_method(x)
print(str("Simpson:\nI" + str(int(x)) + " = " + "%.10f" % result))
if (str("%.10f" % (result - (math.pi / 2))) == "-0.0000000000"):
print(str("diff = " + "%.10f" % (-1 * (result - (math.pi / 2)))))
else:
print(str("diff = " + "%.10f" % (result - (math.pi / 2))))
else:
sys.stdout.write("Not an integer number\n")
exit(84)
else:
sys.stdout.write("Invalid number of argument\n")
exit(84)