## Look at the log-product function to solve Backwards Euler
f = x -> x*exp(x) - (0.1)*(1.0+0.1)
## Find a region where f has a different sign
a = 0.0
b = 10.0
@show f(a)
@show f(b)
## Illustrate bisection
using Plots
pyplot()
a = 0.0
b = 10.0
anim = @animate for i=1:10
plot(f,linspace(a,b,50))
fa = f(a)
fb = f(b)
ab2 = 0.5*a + 0.5*b
fab2 = f(ab2)
if sign(fab2*fb) <= 0
a = ab2
else
b = ab2
end
yl = max(abs(fa),abs(fb)) # show symmetric y
ylims!(-yl,yl)
title!(@sprintf("iter = %i",i))
end
gif(anim, "bisection.gif", fps = 1)
"""
`bisection`
===========
-`bisection(f, a, b)`
-`bisection(f, a, b, delta)`
Find a root of the scalar function f, given that f(a) and f(b) have
different signs and f is continuous. The value x will be within
delta of a root. The method is repeated interval bisection.
"""
function bisection(f,a,b,delta)
fa = f(a)
fb = f(b)
@assert(sign(fa*fb) <= 0)
maxit = 150 # see more on why so large
for i=1:maxit
ab2 = 0.5*a + 0.5*b
fab2 = f(ab2)
if abs(fab2) <= eps(1.0)
break
end
if sign(fab2*fb) <= 0
a = ab2
fa = fab2
else
b = ab2
fb= fab2
end
if abs(b-a) <= delta
break
end
end
return 0.5*a + 0.5*b
end
# set default
bisection(f,a,b) = bisection(f,a,b,eps(1.0))
@show x = bisection(f,0.0,10.0)
@show f(x)
@show f(x+eps(1.0))