381 - A6 - Muller's Method

461 days ago by Professor381

# f = x^2 - 16 # exact = 4 f = tan(x/4)-1 exact = pi x0 = -1 x1 = 2 x2 = 4 xs = [x0,x1,x2] points = plot(f,(x,min(xs),max(xs))) points += point((x0,f(x=x0)),color='red',size=40)+point((x1,f(x=x1)),color='red',size=40)+point((x2,f(x=x2)),color='red',size=40) show(points) 
       
p2 = x2 p1 = x1 p0 = x0 fx2 = f(x=p2) fx1 = f(x=p1) fx0 = f(x=p0) errors = [] for k in range(5): p1mp2 = p1-p2 p0mp1 = p0-p1 p0mp2 = p0-p2 c = n(fx2) b = n(p0mp2*(fx1-fx2)/(p1mp2*p0mp1)- p1mp2*(fx0-fx2)/(p0mp2*p0mp1)) a = n((fx0-fx2)/(p0mp2*p0mp1)-(fx1-fx2)/(p1mp2*p0mp1)) d = sqrt(b^2-4*a*c) if (abs(b-d) <= abs(b+d)): delx = 2*c/(b+d) else: delx = 2*c/(b-d) end p0 = p1 fx0 = fx1 p1 = p2 fx1 = fx2 p2 = p2 - delx # set the new estimate fx2 = f(x=p2) points += point((p2,f(x=p2)),color='orange',size=40) errors.append(n(abs(p2-exact)/exact)) pretty_print(n(p2)) end show(points) pretty_print((p2,fx2)) pretty_print(' or ',(n(p2),n(fx2))) pretty_print(errors) 
       

                                
                            

                                
# Experiment with r to see if there is a (relatively) constant value for the following ratio of successive error terms r = 1.84 for k in range(len(errors)-1): errors[k+1]/errors[k]^r 
       
0.942601562340136
1.59555369831210
0.552756307812208
1.05055477603422
0.942601562340136
1.59555369831210
0.552756307812208
1.05055477603422