On Monday, February 24, 2014 7:48:50 PM UTC-8, Jason Merrill wrote:
>
> On Monday, February 24, 2014 3:34:22 AM UTC-8, Hans W Borchers wrote:
>
>> BTW It would be nice to have Ridders' algorithm available in Julia, too,
>> about which the "Numerical Recipes" say:
>>
>> "*In both reliability and speed, Ridders's method is generally
>> competitive with the more highly developed and better established (but more
>> complicated) method of Van Wijngaarden, Dekker, and Brent [...].*"
>>
>> My Julia version of Ridders' method would be as follows, perhaps it can
>> be improved with your proposed x1 < (x1 + x2)/2 < x2 :
>>
>
> I don't know whether the stopping criterion in your implementation is fine
> or not, but it makes me a little suspicious. There aren't any pairs of
> floats much larger than 1.0 that are separated by less than eps(), so it
> looks to me like if you try to isolate a root larger than 1.0, this
> algorithm will run either for maxiter iterations, or until the upper and
> lower bounds exactly coincide. Exact coincidence might be fine, or it might
> not. What I like about bisection is that it's so simple that I can analyze
> it confidently. The *only* arithmetic in bisect_root is a right shift in
> the _middle subroutine.
>
> I'd be curious to see ridder's method benchmarked against the bisection
> code for a few different objective functions. For simple objective
> functions evaluated over floats, it's harder to outperform bisection than a
> lot of people think because bisection has such low per-iteration overhead.
> The square root in here might hurt it some.
>
Because I can't leave well enough alone...
julia> function ridder_loop()
accum = 0.0
for i = 1:10000
accum += ridders(sin, 2.0, 4.0)
end
accum
end
julia> function bisect_loop()
accum = 0.0
for i = 1:10000
accum += bisect_root(sin, 2.0, 4.0)[1]
end
accum
end
julia> @time ridder_loop()
elapsed time: 10.430907316 seconds (1596537204 bytes allocated)
31415.92653589205
julia> @time ridder_loop()
elapsed time: 10.440722588 seconds (1596480048 bytes allocated)
31415.92653589205
julia> @time bisect_loop()
elapsed time: 0.152900673 seconds (25984456 bytes allocated)
31415.92653589205
julia> @time bisect_loop()
elapsed time: 0.144590689 seconds (25920064 bytes allocated)
31415.92653589205
So for this example, admittedly with an overly simplistic objective
function, bisection is winning by a factor of about 70.
> #-- --------------------------------------------------------------------
>> function ridders(f::Function, a::Real, b::Real;
>> maxiter::Integer = 1000, tol::Real = eps())
>>
>> x1 = a; x2 = b
>> f1 = f(x1); f2 = f(x2)
>> if f1 * f2 > 0
>> error("f(a) and f(b) must have different signs.")
>> elseif f1 == 0
>> return x1
>> elseif f2 == 0
>> return f2
>> end
>>
>> niter = 2
>> while niter < maxiter && abs(x1 - x2) > tol
>> xm = (x1 + x2)/2.0
>> fm = f(xm)
>> if fm == 0; return xm; end
>>
>> x3 = xm + (xm - x1) * sign(f1 - f2) * fm / sqrt(fm^2 - f1 * f2)
>> f3 = f(x3)
>> niter += 2
>> if f3 == 0; return x3; end
>>
>> if fm * f3 < 0
>> x1 = xm; f1 = fm
>> x2 = x3; f2 = f3
>> elseif f1 * f3 < 0
>> x2 = x3; f2 = f3
>> elseif f2 * f3 < 0
>> x1 = x3; f1 = f3
>> else
>> error("Inform the maintainer: you should never get here.")
>> end
>> end
>>
>> return (x1 + x2) / 2.0
>> end
>> #-- --------------------------------------------------------------------
>>
>>
