-
Notifications
You must be signed in to change notification settings - Fork 0
/
numerical_methods.py
186 lines (155 loc) · 8.54 KB
/
numerical_methods.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
# This program is to implement the various methods taught in MA202 by Dr T.N.Janakiraman
# DEFINE FUNCTION IN THE CODE for now
# the polynomial is defined here
def poly(x):
return (x*x*x)-(x)- 1
found = stop = False
digits = 4
pi = 3.1415926535897932384
# this function will give the x axis'x value at the point in the real line where the chord connecting the points (a,poly(a)) and (b,poly(b)) cross the xaxis
def chord_xintercept(a,b):
return a + ( -poly(a) * (b-a) )/( poly(b) - poly(a) )
def regulafalsi(start,end):
var = "junk"
counter = 0 # to catch var staying on the same side for more than 2 iterations - Illinois algorithm
while(True):
temp = var
# the coreof this function,
var = chord_xintercept(start,end)
# we need to check if the same endpoint persists for more than 2 iterations ie Illinois algorithm
if(counter >= 2):
# we need to bias the next iteration towards start, we have had the last 2 iterations giving us the var towards end's side
var = ( poly(end)*start - 0.5*poly(start)*end )/( poly(end) - 0.5*poly(start) )
if(counter <= -2):
# we need to bias the next iteration towards end, we have had the last 2 iterations giving us the var towards start's side
var = ( 0.5*poly(start)*end - poly(end)*start )/( 0.5*poly(start) - poly(end) )
# if the maximum precision possible has been reached
if(var == temp):
return var
# if we are lucky enough to be looking for a root
# that can be described using a finte number of digits, and find it
if(poly(start) == 0):
return start
if(poly(end) == 0):
return end
# if we have reached the last level of precision possible
if(var == "done"):
if( (start - end) <= 10**(-1*digits) ):
print (start - end)
return start
# The root lies in between start and end, never play outside the boundaries :p
assert (poly(start)*poly(end) <= 0)
# check on which side our var is on towards start or end?
if( poly(var) * poly(end) > 0):
# its on end's side, thus the new end is
end = var
counter = counter + 1
else:
# then var has to lie on the other side
start = var
counter = counter - 1
# for debuggind and letting the user know whats happening
#print "\nThe start,end,var is "+str(start)+" "+str(end)+" "+str(var)
#print "The start,end,var is "+str(poly(start))+" "+str(poly(end))+" "+str(poly(var))
print "\nThe root lies between "+str(start)+" AND "+str(end)+" \n working... shortening the range"
print "\n counter is "+str(counter)
# the recursive function to find the root using Bisection Method-
# what we do is using the range the user helped us find, we start by
# finding the mid of that range and then check where this mid lies, wrt the root
# using this info we shorten the range by removing the redundant side of the range.
# we continue doing this by now finding the mid of the newly formed range.
def bisection(start,end):
# recursion is not as fast I think. It will surely take more space as we go deeper,
# I experienced "RuntimeError: maximum recursion depth exceeded while calling a Python object" thus i am now using a loop
while(True):
# the coreof this function, we use 2.0 to get var as a float
var = (start+end)/2.0
# if we are lucky enough to be looking for a root
# that can be described using a finte number of digits, and find it
if(poly(start) == 0):
return start
if(poly(end) == 0):
return end
# if we have reached the last level of precision possible
if(start == end):
return start
# Let us establish an alternate base condition, if the required level of precision has been reached,
if(poly(var) <= 10**(-1*digits) and poly(var) >= (-1*(10**(-1*digits)))):
return var
# The root lies in between start and end, never play outside the boundaries :p
assert (poly(start)*poly(end) <= 0)
# check on which side our var is on towards start or end?
if( (poly(var)<0 and poly(start)<0) or (poly(var)>0 and poly(start)>0) ):
# its on start's side, doesnt matter wheather thats positive or negative
# we know that the root is between var and end, the new start is
start = var
else:
# then var has to lie on the other side
end = var
# for debuggind and letting the user know whats happening
#print "\nThe start,end,var is "+str(start)+" "+str(end)+" "+str(var)
#print "The start,end,var is "+str(poly(start))+" "+str(poly(end))+" "+str(poly(var))
print "\nThe root lies between "+str(start)+" AND "+str(end)+" \n working... shortening the range"
while(stop == False):
# At a root of odd multiplicity the graph of the function crosses the x-axis, whereas at a root of
# even multiplicity the graph touches the x-axis. the first 2 methods can only deal with roots of even mutiplicity
print "Enter 1 to find the root of a polynomial by the Bisection Method"
print "Enter 2 to find the root of a polynomial by the Regula-Falsi"
resp = raw_input()
if(resp == "1"):
# now to find the region where the function changes sign
print "\nTo start the approximation we need to establish a region (a,b) where the root lies \nthis method can only find roots which have even multiplicity ie, the curve crosses the graph at that point \n"
while(found == False):
# catch ValueError caused when user enters invalid a or b
try:
print "Enter a possible 'a'"
a = float(raw_input())
print "Enter a possible 'b'"
b = float(raw_input())
except ValueError:
print "You have entered a Invalid value for a or b, please try again."
continue
if((poly(a)*poly(b)) < 0):
print "Thank you There lies a root in this range! we shall find it."
found = True
else:
print "\nSorry Please try again. use Control + C to exit the program."
# Get the precision to which the root must be found NO NEED CAUSE WE BETTER AIM FOR SUPER PRECISION
print "Enter the number of digits after the decimal point to which the root should be found (maximum 11)"
digits = int(raw_input())
ans = bisection(a,b)
# we truncate the non-significant digits
print "\n\nThe root has been found to be approximately "+str(ans)[0:len(str(int(ans)))+1+digits]
elif(resp == "2"):
# now to find the region where the function changes sign
print "\nTo start the approximation we need to establish a region (a,b) where the root lies \nthis method can only find roots which have even multiplicity ie, the curve crosses the graph at that point \n"
while(found == False):
# catch ValueError caused when user enters invalid a or b
try:
print "Enter a possible 'a'"
a = float(raw_input())
print "Enter a possible 'b'"
b = float(raw_input())
except ValueError:
print "You have entered a Invalid value for a or b, please try again."
continue
if((poly(a)*poly(b)) < 0):
print "Thank you There lies a root in this range! we shall find it."
found = True
else:
print "\nSorry Please try again. use Control + C to exit the program."
# Get the precision to which the root must be found NO NEED CAUSE WE BETTER AIM FOR SUPER PRECISION
print "Enter the number of digits after the decimal point to which the root should be found (maximum 11)"
digits = int(raw_input())
ans = regulafalsi(a,b)
# we truncate the non-significant digits
print "\n\nThe root has been found to be approximately "+str(ans)[0:len(str(int(ans)))+1+digits]
else:
print "You have entered a Invalid choice"
print "Do you want to exit the program? (y/n)"
resp2 = raw_input()
if(resp2 == "y"):
stop = True
else:
# reset the flag so that the user can look for the root in a different range
found = False