# 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 |
|
|