MathGroup Archive 2010

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: FindRoots?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg112213] Re: FindRoots?
  • From: "David Park" <djmpark at comcast.net>
  • Date: Sun, 5 Sep 2010 05:26:34 -0400 (EDT)

First, here is a packaged version of Daniel's method, which might be more
convenient to use:

ClearAll[RealFunctionRoots];
RealFunctionRoots::usage = 
  "RealFunctionRoots[lhs \[Equal] rhs, {var, low, high}, precision, \
wrapper] will return roots of real equations in the specified real \
domain.\nRealFunction[expr, {var, low, high}] calculates the roots \
for expr \[Equal] 0.\nprecision specifies the precision of real \
variables to be used in the display, The default is MachinePrecision. \
It is necessary to set the precision of real approximate constants \
using number marks.\nwrapper is a function, default None, that wraps \
the Reduce statement. It can be used to implement extra specification \
for the Reduce routine.\nOptions for Reduce may be passed. The \
message Reduce::ratnz is turned off for the duration of the \
calculation.";

Options[RealFunctionRoots] = Options[Reduce];
SyntaxInformation[
   RealFunctionRoots] = {"ArgumentsPattern" -> {_, {_, _, _}, _., _., 
     OptionsPattern[]}};
RealFunctionRoots[equation_, {x_, low_, high_}, 
  precision : (_Integer?Positive | MachinePrecision) : 
    MachinePrecision, wrapperfunction_: None, 
  opts : OptionsPattern[]] :=
 Module[{work, workequation},
  Off[Reduce::ratnz];
  workequation = If[Head[equation] === Equal, equation, equation == 0];
  work =
   If[wrapperfunction === None,
    Reduce[{workequation, low <= x <= high}, x, Reals, opts],
    wrapperfunction[
     Reduce[{workequation, low <= x <= high}, x, Reals, opts]]];
  work = N[{ToRules[work]}, precision] // Chop;
  On[Reduce::ratnz];
  work
  ]

It only handles real functions, but this is the common case. (So it doesn't
handle the Zeta zeros case, but that is a little quirky because we have a
Mathematica routine to get them directly.)

This will handle many common cases very well, but not everything and not
always efficiently. Gianluca's example appears to be a case that exposes a
bug in Reduce. Here is another case where Reduce is quite inefficient.
First, with RootSearch:

Needs["Ersek`RootSearch`"] 

f9[x_] := BesselJ[5, x] - Sin[x] 
Plot[f9[x], {x, -1, 10 \[Pi]}] 
Timing[roots = RootSearch[f9[x] == 0, {x, -1, 10 \[Pi]}]] 

{0.078, {{x -> 0.}, {x -> 3.09271}, {x -> 6.6617}, {x -> 
    9.60579}, {x -> 12.6278}, {x -> 15.7099}, {x -> 18.8209}, {x -> 
    21.946}, {x -> 25.0783}, {x -> 28.2145}, {x -> 31.3531}}}

Then using Reduce:

Timing[roots = RealFunctionRoots[f9[x], {x, -1, 10 \[Pi]}]] 

Reduce::incs: Warning: Reduce was unable to prove that the solution set
found is complete.

{100.824, {{x -> 3.09271}, {x -> 6.6617}, {x -> 9.60579}, {x -> 
    12.6278}, {x -> 15.7099}, {x -> 18.8209}, {x -> 21.946}, {x -> 
    25.0783}, {x -> 28.2145}, {x -> 31.3531}}}

Not only is the solution set incomplete, but Reduce was almost 1300 times
slower.

1) The new Reduce/Root functionality was new with Version 7. Before that
there was not an adequate method for extracting multiple roots of a real
function on a real domain. Ted's RootSearch routine was based on the
literature, was professionally done, is robust and served a real need. And
still does.

2) To say that RootSearch is "amateurish", not "bleeding edge" and unworthy
of being in Mathematica is not helpful. Is it let them eat cornbread and
keep silent, until the elites are ready to hand them cake? 

3) The Reduce/Root method, for all its marvels, is not all that it's cracked
up to be. Already cases have been found where it does not perform at all, is
slow, or performs erratically. Dollars to donuts that more cases will pop
up. Furthermore, it is not packaged in an accessible way nor well documented
(linking from FindRoot say). Few users would think of using Reduce, or know
at first what to do with Root objects. This is a case where people may have
a simple straightforward problem to solve but are trapped into an afternoon
of trying to figure out what is going on, or guessing at various options and
permutations without any clue as to what will get them to the answer, or
even if there is a method. RootSearch is far more intuitive. There are no
references or documentation of the "bleeding edge" technology of Reduce and
even if there were, are users who want to find roots of a real function
really going to master it? Not to speak of any innovations or quirks that
may have been introduced by WRI programmers.

4) I consider Interpolating function to be mathematical objects. They can be
fairly accurate and are certainly suitable objects for root finding. There
are many applied problems that fall in this class. To suggest that users
shouldn't use the Mathematica interpolating functions is not helpful to
users who might actually waste a lot of their time following such advice.

So give us both advanced equation solving algorithms and a good numerical
root finder. They both have a place in Mathematica.


David Park
djmpark at comcast.net
http://home.comcast.net/~djmpark/  




From: Gianluca Gorni [mailto:gianluca.gorni at uniud.it] 


Very bizarre! Try the following inputs in sequence,
in a fresh kernel:

eq == 2 x + Log[-((-1 + 2 x)/(-1 + 2 x^2))] ==== 0;
Reduce[eq && -2 < x < -1/Sqrt[2], x, Reals]
Reduce[eq && x < 0, x, Reals]
Reduce[eq && -2 < x < -1/Sqrt[2], x, Reals]

I get the first Reduce output correct, but the next two outputs
are False.
If I re-evaluate the inputs in the same kernel,
I get False three times.

Strange results also from the following, in a fresh kernel:

eq == 2 x + Log[-((-1 + 2 x)/(-1 + 2 x^2))] ==== 0;
Reduce[eq && -2 < x < -1/Sqrt[2], x, Reals]
Reduce[eq, x, Reals]
Reduce[eq && -2 < x < -1/Sqrt[2], x, Reals]

$Version
7.0 for Mac OS X x86 (64-bit) (February 19, 2009)

Best regards,
Gianluca Gorni



On 03/set/2010, at 14.38, Bob Hanlon wrote:

>
> $Version
>
> 7.0 for Mac OS X x86 (64-bit) (February 19, 2009)
>
> Reduce[2 x + Log[-((-1 + 2 x)/(-1 + 2 x^2))] ==== 0 &&
>     -2 < x < -1/Sqrt[2], x, Reals] //
>  N // ToRules
>
> {x->-0.861936}
>
> sol == x /. First@NDSolve[
>     {x'[t] ==== x[t] + 2, x[0] ==== -1},
>      x, {t, 0, 2}];
>
> Minimize[{Abs[sol[t]], 0 < t < 2}, t][[2]]
>
> {t->0.693147}
>
>
> Bob Hanlon
>
> ---- Gianluca Gorni <gianluca.gorni at uniud.it> wrote:
>
> ==========================
>
> In my opinion Reduce can replace RootSearch in some
> cases but not in others.
>
> First of all, Reduce has bugs. Here is an analytic function
> that clearly has a real root:
>
> Plot[2 x + Log[-((-1 + 2 x)/(-1 + 2 x^2))], {x, -2, -1/Sqrt[2]}]
>
> Still, Reduce does not see it (as of version 7.0.1):
>
> Reduce[2 x + Log[-((-1 + 2 x)/(-1 + 2 x^2))] ======== 0 && -2 <
>   x < -1/Sqrt[2], x, Reals]
> False
>
> (I reported this example to wolfram last year).
>
> Next, Reduce has problems with inexact input, and with
> InterpolatingFunction, so that it won't work with
> the output of NDSolve:
>
> sol ==== x /.
>  First@NDSolve[{x'[t] ======== x[t] + 2, x[0] ======== -1}, x, {t, 0, 2}]=
;
> Reduce[sol[t] ======== 0 && 0 < t < 2, t]
>
> or, say, with functions obtained by interpolating between Locators.
> RootSearch works fine in these cases.
>
> I have made some interactive panels where I can change the Locators
> with the mouse and I get in real time the roots of the interpolating
> function as big Points in the plot: I can do this with RootSearch,
> but not with Reduce.
> Unfortunately, I can't give these panels to other users, because
> I can't assume that they have RootSearch installed.
>
> I endorse the wish that the functionality of RootSearch were available
> in the kernel.
>
> Best regards,
> Gianluca Gorni
>
>
>
> On 01/set/2010, at 
>>
>> On 30 Aug 2010, at 12:19, David Park wrote:
>>
>>> I don't know why Wolfram doesn't acquire the rights to this package and
>>> incorporate it as part of regular Mathematica. Finding all the roots of=
 a
>>> real function on a 1-dimensional domain is the most common case of root
>>> finding, and it is the one thing regular Mathematica is very poor at. T=
ed's
>>> package is quite robust and returns all the roots in order.
>>
>>
>> Mathematica's Reduce can do this (in the case of complex analytic functi=
ons, which covers most important practical cases) since version 7. I have p=
osted examples of this several times, but somehow nobody seems to notice an=
d we keep reading strange remarks like the above. Well, I don't want it to =
sound like I am trying to diminish Ted's achievement, but quite frankly, fr=
om the mathematical point of view, the methods used in the RootSearch packa=
ge are rather primitive, relatively to the current state of knowledge in th=
is field, of course. Not only does Reduce solve most transcendental equatio=
ns much faster but it also does not require setting the values of any optio=
ns by hand, and, what I think is most important, the results thus obtained =
are provably correct, which is not the case with Ted's package. In fact, mu=
ch more general methods of solving transcendental equations exist (in gener=
al the functions need not be analytic and can involve systems of n-real equ=
ations in n-vari
> a!
> b
>> l!
>> es with non-vanishing Jacobian). An example of such a method can be foun=
d in one of my demonstrations on the demonstrations site (the algorithm can=
 easily be converted to a practical method for real life probl ems.
>>
>> In view of this, how could there be any justification for the suggestion=
 that Wolfram incorporates what is actually an impressive but amateurish pa=
ckage into what is supposed to be a sophisticated mathematical program?
>> Sometimes certain capabilities are missing from Mathematica simply becau=
se implementing "bleeding edge" algorithms in a way that is satisfactory fo=
ra program that has the kind of aspirations that Mathematica obviously does=
 takes time.
>>
>> These sort of "helpful suggesting", particularly when coming from non-ex=
p==
> erts, have no chance whatever of being accepted and even can produce the =
im==
> pression of an ulterior motive being involved.
>>
>



  • Prev by Date: Mathlink problem in C
  • Next by Date: Re: problem with RandomInteger
  • Previous by thread: Re: FindRoots?
  • Next by thread: Re: FindRoots?