Re: Re: RootSearch Performance
- To: mathgroup at smc.vnet.net
- Subject: [mg72421] Re: [mg72382] Re: RootSearch Performance
- From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
- Date: Wed, 27 Dec 2006 05:43:10 -0500 (EST)
- References: <emggf1$2jk$1@smc.vnet.net> <200612231114.GAA20929@smc.vnet.net>
On 23 Dec 2006, at 20:14, akozlowski at gmail.com wrote:
>
>
> On Dec 22, 8:44 pm, "Ted Ersek" <ted.er... at tqci.net> wrote:
>> About a week ago we had the thread:
>> "FindRoot anomoly (example from Mathematica"
>> --------------------------------------------------------------
>> In that thread galordloo wrote:
>>
>> Ersek's RootSearch function finds only seven roots to the equation
>> between x == 1 and x == 100:
>>
>> sol == RootSearch[3*Cos[x] ==== Log[x], {x, 1, 100}]
>>
>> <snip> (* Rules for seven roots, all between 1 and 20 were
>> returned. *)
>>
>> ---------------------------------
>> Then Carl Woll (of Wolfram Research) replied:
>>
>> An alternative method is possible using IntervalBisection:
>> < snip >
>>
>> The nice thing about IntervalBisection is that we were able to use an
>> initial range of {0,Infinity} instead of {1,100}. The other nice
>> thing
>> about IntervalBisection is that we are guaranteed that all roots
>> lie in
>> the interval given in the result, something that is not true with
>> RootSearch.
>>
>> The price you pay is that the only transcendental functions
>> allowed in
>> the input are trigonometric/exponential functions and their inverses,
>> i.e., only functions which support Interval arguments.
>>
>> =====================================================================
>> ====== =
>> I wrote the RootSearch package, and it's designed to never look
>> for a root outside the range specified in the input. While
>> RootSearch has many lines of code, the ideas I use to limit the
>> range being searched are straight forward. I doubt you can find
>> an example where it returns a root outside the specified range.
>> However, if you can find an example that shows otherwise,
>> I would like to know about it.
>>
>> Ted Ersek
>> Reply to: Ted.Er... at navy.mil
>
> I think the problem here is a simple linguistic ambiguity. What Carl
> meant is that IntervalBisection returns a list of sub-intervals of the
> original interval that contain all the roots that lie in the original
> interval. In other words: no roots are ever missed, unlike the case
> with RootSearch that may sometimes miss some roots. This is quite
> different form saying that all rots that are found lie in the
> originally specified interval, which of course is true for RootSearch.
>
> Andrzej Kozlowski
>
I would like to add that all the above in no way implies that
IntervalBisection or any other root interval based root finding
method is better or preferable to RootSearch. They just owrk
diffeently and fail in different ways on certain type of problems.
The "guarantee" that IntervalBisection gives of "not missing any
roots" is not always as impressive as it may sound. Here is a rather
striking example. Consider the polynomial:
In[1]:=
f[x_] = ExpandAll[Product[x - (1 + 1/n), {n, 100, 110}]*
Product[x - (1 - 1/n), {n, 100, 110}]];
This is a polynomial with very densely clustered roots; the kind of
thing that will generally be very difficult to deal with for most
numerical solvers.
Let's first find all the roots using Solve and NSolve. This is of
course not fair to the other solvers, since we are solving a
polynomial equation, and Solve and NSolve are designed specifically
for such equations. Anyway, Solve, as one might expect has no problems:
In[2]:=
ls = x /. N[Solve[f[x] == 0, x]]
Out[2]=
{0.99, 0.9900990099009901, 0.9901960784313726, 0.9902912621359223,
0.9903846153846154, 0.9904761904761905, 0.9905660377358491,
0.9906542056074766, 0.9907407407407407, 0.9908256880733946,
0.990909090909091, 1.009090909090909, 1.0091743119266054,
1.0092592592592593, 1.0093457943925233, 1.009433962264151,
1.0095238095238095, 1.0096153846153846, 1.0097087378640777,
1.0098039215686274, 1.00990099009901, 1.01}
In[3]:=
Length[ls]
Out[3]=
22
NSolve finds the problem much harder and needs high working precision
to get the exact set of roots:
In[4]:=
ls1=Chop[x/.NSolve[f[x]==0,x,WorkingPrecision->100]];
In[5]:=
Length[ls1]
Out[5]=
22
RootSearch also needs high precision. Even using a large number of
initial samples, it finds only two out of the 22 roots in the
interval {-1,2}.
<< Ersek`RootSearch`
In[7]:=
RootSearch[f[x] == 0, {x, -1, 2}, PrecisionGoal -> 10,
InitialPrecision -> 100, InitialSamples -> 1000]
Out[7]=
{{x -> 0.9907407407407407}, {x -> 1.009090909090909}}
As for IntervalBisection, it is even less impressive.
<< NumericalMath`IntervalRoots`
Let's ask IntervalBisection to find all the roots in the interval [0,2]:
In[9]:=
IntervalBisection[f[x], x, Interval[{-1, 2}], 10^(-3),
MaxRecursion -> 50, WorkingPrecision -> 100]
Out[9]=
Interval
[{0.104492187500000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000000452060442841134`100., 2}]
As we see, all the roots are indeed contained in the interval
returned, but the output is hardly a great improvement on the input!
Andrzej Kozlowski
- References:
- Re: RootSearch Performance
- From: akozlowski@gmail.com
- Re: RootSearch Performance