|
4 | 4 | # DEFINE FUNCTION IN THE CODE for now
|
5 | 5 | # the polynomial is defined here
|
6 | 6 | def poly(x):
|
7 |
| - return (5*x*x*x+3*x+4) |
8 |
| - |
9 |
| - #from math import exp |
10 |
| - #return (exp(x)+x) |
11 |
| - |
12 |
| - #return (3*x*x-2*x-8) |
| 7 | + return (x*x*x)-(x)- 1 |
13 | 8 |
|
14 | 9 | found = stop = False
|
15 | 10 | digits = 4
|
16 | 11 | pi = 3.1415926535897932384
|
17 | 12 |
|
18 | 13 | # 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
|
19 | 14 | def chord_xintercept(a,b):
|
20 |
| - slope = (poly(b)-poly(a))/(b-a) |
21 |
| - |
22 |
| - x_a = a + (-1*poly(a))/slope |
23 |
| - x_b = b - poly(b)/slope |
24 |
| - if(x_a == x_b): |
25 |
| - return x_a |
26 |
| - else: |
27 |
| - # debugging |
28 |
| - print "\nx_a = "+str(x_a)+" "+str(x_a==x_b) |
29 |
| - print "x_b = "+str(x_b)+"\n" |
30 |
| - return "done" |
| 15 | + return a + ( -poly(a) * (b-a) )/( poly(b) - poly(a) ) |
31 | 16 |
|
32 | 17 | def regulafalsi(start,end):
|
| 18 | + var = "junk" |
| 19 | + counter = 0 # to catch var staying on the same side for more than 2 iterations - Illinois algorithm |
33 | 20 | while(True):
|
| 21 | + temp = var |
34 | 22 | # the coreof this function,
|
35 | 23 | var = chord_xintercept(start,end)
|
36 | 24 |
|
| 25 | + # we need to check if the same endpoint persists for more than 2 iterations ie Illinois algorithm |
| 26 | + if(counter >= 2): |
| 27 | + # we need to bias the next iteration towards start, we have had the last 2 iterations giving us the var towards end's side |
| 28 | + var = ( poly(end)*start - 0.5*poly(start)*end )/( poly(end) - 0.5*poly(start) ) |
| 29 | + if(counter <= -2): |
| 30 | + # we need to bias the next iteration towards end, we have had the last 2 iterations giving us the var towards start's side |
| 31 | + var = ( 0.5*poly(start)*end - poly(end)*start )/( 0.5*poly(start) - poly(end) ) |
| 32 | + |
| 33 | + # if the maximum precision possible has been reached |
| 34 | + if(var == temp): |
| 35 | + return var |
| 36 | + |
37 | 37 | # if we are lucky enough to be looking for a root
|
38 | 38 | # that can be described using a finte number of digits, and find it
|
39 | 39 | if(poly(start) == 0):
|
40 | 40 | return start
|
41 | 41 | if(poly(end) == 0):
|
42 | 42 | return end
|
| 43 | + |
43 | 44 | # if we have reached the last level of precision possible
|
44 | 45 | if(var == "done"):
|
45 |
| - if(abs(poly(start))<abs(poly(end))): |
| 46 | + if( (start - end) <= 10**(-1*digits) ): |
| 47 | + print (start - end) |
46 | 48 | return start
|
47 |
| - else: |
48 |
| - return end |
49 |
| - |
50 |
| - # Let us establish an alternate base condition, if the required level of precision has been reached, |
51 |
| - #if(poly(var) <= 10**(-1*digits) and poly(var) >= (-1*(10**(-1*digits)))): |
52 |
| - #return var |
53 |
| - |
| 49 | + |
54 | 50 | # The root lies in between start and end, never play outside the boundaries :p
|
55 | 51 | assert (poly(start)*poly(end) <= 0)
|
56 |
| - |
| 52 | + |
57 | 53 | # check on which side our var is on towards start or end?
|
58 |
| - if( (poly(var)<0 and poly(start)<0) or (poly(var)>0 and poly(start)>0) ): |
59 |
| - # its on start's side, doesnt matter wheather thats positive or negative |
60 |
| - # we know that the root is between var and end, the new start is |
61 |
| - start = var |
| 54 | + if( poly(var) * poly(end) > 0): |
| 55 | + # its on end's side, thus the new end is |
| 56 | + end = var |
| 57 | + counter = counter + 1 |
62 | 58 | else:
|
63 | 59 | # then var has to lie on the other side
|
64 |
| - end = var |
| 60 | + start = var |
| 61 | + counter = counter - 1 |
| 62 | + |
65 | 63 | # for debuggind and letting the user know whats happening
|
66 | 64 | #print "\nThe start,end,var is "+str(start)+" "+str(end)+" "+str(var)
|
67 | 65 | #print "The start,end,var is "+str(poly(start))+" "+str(poly(end))+" "+str(poly(var))
|
68 | 66 | print "\nThe root lies between "+str(start)+" AND "+str(end)+" \n working... shortening the range"
|
| 67 | + print "\n counter is "+str(counter) |
69 | 68 |
|
70 | 69 | # the recursive function to find the root using Bisection Method-
|
71 | 70 | # what we do is using the range the user helped us find, we start by
|
|
0 commit comments