Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

Re: FindRoots?

38 views
Skip to first unread message

Daniel Lichtblau

unread,
Aug 31, 2010, 4:14:22 AM8/31/10
to
David Park wrote:
> Install Ted Ersek's package RootSearch from:
>
> http://library.wolfram.com/infocenter/MathSource/4482/
>
> 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. Ted's
> package is quite robust and returns all the roots in order.
>
>
> David Park
> djm...@comcast.net
> http://home.comcast.net/~djmpark/
>
>
>
> From: Sam Takoy [mailto:sam....@yahoo.com]
>
> Hi,
>
> Is there a command for numerically finding all roots of a function in a
> given interval. If not, what's the best way of accomplishing this task?
>
> Many thanks in advance,
>
> Sam

If I instead define RootSearch to use Reduce, as below, it handles quite
well all the examples from the Examples section of
RootSearchExamples.nb. I show them below. I modified one to use
Piecewise instead of a condition, and for the last I use Reduce directly
to get an expanded form of the solution set.

I think Mathematica is quite good at this sort of thing, actually.

Daniel Lichtblau
Wolfram Research

-------------------------------

In[33]:= RootSearch[
eqn_, {x_, lo_, hi_}] := {ToRules[Reduce[{eqn, lo <= x <= hi}]]}

In[34]:= f1[x_] := Sin[4 x] - (x + 1)/8
RootSearch[f1[x] == 0, {x, -13, 12}]

Out[35]= {{x ->
Root[{1 - 8 Sin[4 #1] + #1 &, -8.3482893782997216496}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, -8.1289315341182707392}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, -6.8629511940895252882}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, -6.4714706256130258036}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, -5.3539165396979340111}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, -4.8374622685220237568}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, -3.8363831890811367726}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, -3.21161777565464349626}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, -2.3149160596763252331}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, -1.58922633877627429127}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, -0.79190196071777069054}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, 0.032351189053103634282}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, 0.73087710353047892233}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, 1.65538160232955244046}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, 2.2515547207396359640}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, 3.2828240626757282491}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, 3.7673870629720211202}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, 4.9206885635879098258}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, 5.2724915094753945872}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, 6.5961054327905341135}]}, {x ->
Root[{1 - 8 Sin[4 #1] + #1 &, 6.7398184613474433969}]}}

In[36]:= f2[x_] := If[x < 3/2, Sqrt[3/2 - x] - Exp[-4 x], 2* x Exp[-x]]
RootSearch[f2[x] == 0, {x, -25, 25}]

Out[37]= {{x ->
Root[{-Sqrt[2] +
E^(4 #1) Sqrt[3 - 2 #1] &, -0.055200564731044670903}]}, {x ->
Root[{-Sqrt[2] + E^(4 #1) Sqrt[3 - 2 #1] &,
1.49999385548561363888}]}}

In[38]:= RootSearch[Exp[x - \[Pi]] == 1 - \[Pi] + x, {x, -3, 4}]

Out[38]= {{x -> \[Pi]}}

In[39]:= f4[x_] := (1.5 - Erf[x] - Erf[2 - x]) Exp[-Abs[x - 1]]
RootSearch[f4[x] == 0, {x, -400, 200}]

During evaluation of In[39]:= Reduce::ratnz: Reduce was unable to solve
the system with inexact coefficients. The answer was obtained by solving
a corresponding exact system and numericizing the result. >>

Out[40]= {{x -> 0.517891}, {x -> 1.48211}}

In[41]:= f5[x_] :=
3.2 - 3*Cos[x] - Exp[-(x - 8 \[Pi])^2] - Exp[-(x - 14 \[Pi])^2]
s5 = RootSearch[f5[x] == 0, {x, 0, 30 \[Pi]}]

During evaluation of In[41]:= Reduce::ratnz: Reduce was unable to solve
the system with inexact coefficients. The answer was obtained by solving
a corresponding exact system and numericizing the result. >>

Out[42]= {{x -> 24.5432}, {x -> 25.7223}, {x -> 43.3928}, {x ->
44.5718}}

In[43]:= f6[x_] :=
Piecewise[{{ Sqrt[Abs[x] - 1.2], x <= 1.4}}, Indeterminate]
RootSearch[f6[x] == 0, {x, -8, 8}]

During evaluation of In[43]:= Reduce::ratnz: Reduce was unable to solve
the system with inexact coefficients. The answer was obtained by solving
a corresponding exact system and numericizing the result. >>

Out[44]= {{x -> -1.2}, {x -> 1.2}}

In[45]:= RootSearch[Zeta[1/2 + I y] == 0, {y, 0, 237}]

Out[45]= {ToRules[
C[1] \[Element] Integers && 1 <= C[1] <= 100 &&
y == -(1/2) I (-1 + 2 ZetaZero[C[1]])]}

In[48]:= explicitsol =
Resolve[Exists[C[1],
Reduce[{Zeta[1/2 + I y] == 0, 0 <= y <= 237}, y]]]

Out[48]= y == Im[ZetaZero[1]] || y == Im[ZetaZero[2]] ||
y == Im[ZetaZero[3]] || y == Im[ZetaZero[4]] ||
y == Im[ZetaZero[5]] || y == Im[ZetaZero[6]] ||
y == Im[ZetaZero[7]] || y == Im[ZetaZero[8]] ||
y == Im[ZetaZero[9]] || y == Im[ZetaZero[10]] ||
y == Im[ZetaZero[11]] || y == Im[ZetaZero[12]] ||
y == Im[ZetaZero[13]] || y == Im[ZetaZero[14]] ||
y == Im[ZetaZero[15]] || y == Im[ZetaZero[16]] ||
y == Im[ZetaZero[17]] || y == Im[ZetaZero[18]] ||
y == Im[ZetaZero[19]] || y == Im[ZetaZero[20]] ||
y == Im[ZetaZero[21]] || y == Im[ZetaZero[22]] ||
y == Im[ZetaZero[23]] || y == Im[ZetaZero[24]] ||
y == Im[ZetaZero[25]] || y == Im[ZetaZero[26]] ||
y == Im[ZetaZero[27]] || y == Im[ZetaZero[28]] ||
y == Im[ZetaZero[29]] || y == Im[ZetaZero[30]] ||
y == Im[ZetaZero[31]] || y == Im[ZetaZero[32]] ||
y == Im[ZetaZero[33]] || y == Im[ZetaZero[34]] ||
y == Im[ZetaZero[35]] || y == Im[ZetaZero[36]] ||
y == Im[ZetaZero[37]] || y == Im[ZetaZero[38]] ||
y == Im[ZetaZero[39]] || y == Im[ZetaZero[40]] ||
y == Im[ZetaZero[41]] || y == Im[ZetaZero[42]] ||
y == Im[ZetaZero[43]] || y == Im[ZetaZero[44]] ||
y == Im[ZetaZero[45]] || y == Im[ZetaZero[46]] ||
y == Im[ZetaZero[47]] || y == Im[ZetaZero[48]] ||
y == Im[ZetaZero[49]] || y == Im[ZetaZero[50]] ||
y == Im[ZetaZero[51]] || y == Im[ZetaZero[52]] ||
y == Im[ZetaZero[53]] || y == Im[ZetaZero[54]] ||
y == Im[ZetaZero[55]] || y == Im[ZetaZero[56]] ||
y == Im[ZetaZero[57]] || y == Im[ZetaZero[58]] ||
y == Im[ZetaZero[59]] || y == Im[ZetaZero[60]] ||
y == Im[ZetaZero[61]] || y == Im[ZetaZero[62]] ||
y == Im[ZetaZero[63]] || y == Im[ZetaZero[64]] ||
y == Im[ZetaZero[65]] || y == Im[ZetaZero[66]] ||
y == Im[ZetaZero[67]] || y == Im[ZetaZero[68]] ||
y == Im[ZetaZero[69]] || y == Im[ZetaZero[70]] ||
y == Im[ZetaZero[71]] || y == Im[ZetaZero[72]] ||
y == Im[ZetaZero[73]] || y == Im[ZetaZero[74]] ||
y == Im[ZetaZero[75]] || y == Im[ZetaZero[76]] ||
y == Im[ZetaZero[77]] || y == Im[ZetaZero[78]] ||
y == Im[ZetaZero[79]] || y == Im[ZetaZero[80]] ||
y == Im[ZetaZero[81]] || y == Im[ZetaZero[82]] ||
y == Im[ZetaZero[83]] || y == Im[ZetaZero[84]] ||
y == Im[ZetaZero[85]] || y == Im[ZetaZero[86]] ||
y == Im[ZetaZero[87]] || y == Im[ZetaZero[88]] ||
y == Im[ZetaZero[89]] || y == Im[ZetaZero[90]] ||
y == Im[ZetaZero[91]] || y == Im[ZetaZero[92]] ||
y == Im[ZetaZero[93]] || y == Im[ZetaZero[94]] ||
y == Im[ZetaZero[95]] || y == Im[ZetaZero[96]] ||
y == Im[ZetaZero[97]] || y == Im[ZetaZero[98]] ||
y == Im[ZetaZero[99]] || y == Im[ZetaZero[100]]

In[49]:= N[explicitsol]

Out[49]= y == 14.1347 || y == 21.022 || y == 25.0109 || y == 30.4249 ||
y == 32.9351 || y == 37.5862 || y == 40.9187 || y == 43.3271 ||
y == 48.0052 || y == 49.7738 || y == 52.9703 || y == 56.4462 ||
y == 59.347 || y == 60.8318 || y == 65.1125 || y == 67.0798 ||
y == 69.5464 || y == 72.0672 || y == 75.7047 || y == 77.1448 ||
y == 79.3374 || y == 82.9104 || y == 84.7355 || y == 87.4253 ||
y == 88.8091 || y == 92.4919 || y == 94.6513 || y == 95.8706 ||
y == 98.8312 || y == 101.318 || y == 103.726 || y == 105.447 ||
y == 107.169 || y == 111.03 || y == 111.875 || y == 114.32 ||
y == 116.227 || y == 118.791 || y == 121.37 || y == 122.947 ||
y == 124.257 || y == 127.517 || y == 129.579 || y == 131.088 ||
y == 133.498 || y == 134.757 || y == 138.116 || y == 139.736 ||
y == 141.124 || y == 143.112 || y == 146.001 || y == 147.423 ||
y == 150.054 || y == 150.925 || y == 153.025 || y == 156.113 ||
y == 157.598 || y == 158.85 || y == 161.189 || y == 163.031 ||
y == 165.537 || y == 167.184 || y == 169.095 || y == 169.912 ||
y == 173.412 || y == 174.754 || y == 176.441 || y == 178.377 ||
y == 179.916 || y == 182.207 || y == 184.874 || y == 185.599 ||
y == 187.229 || y == 189.416 || y == 192.027 || y == 193.08 ||
y == 195.265 || y == 196.876 || y == 198.015 || y == 201.265 ||
y == 202.494 || y == 204.19 || y == 205.395 || y == 207.906 ||
y == 209.577 || y == 211.691 || y == 213.348 || y == 214.547 ||
y == 216.17 || y == 219.068 || y == 220.715 || y == 221.431 ||
y == 224.007 || y == 224.983 || y == 227.421 || y == 229.337 ||
y == 231.25 || y == 231.987 || y == 233.693 || y == 236.524

David Park

unread,
Aug 31, 2010, 4:17:23 AM8/31/10
to
True, Sin[1/x] will eventually break any root finding algorithm because we
will always run out of time or space to store the answers if we push too
close to zero.

However, this does provide a good example to show how well Ted's algorithm
works. In this case the exact roots are 1/(n Pi) where n is a positive
interger that runs between limits that span the region we are searching.

So, let's try some cases numerically. To search for roots between 0.01 and
0.1 n runs from 4 to 21 and there are 28 roots. Here are the results from
RootSearch.

Needs["Ersek`RootSearch`"]

RootSearch[Sin[1/x] == 0, {x, 0.01, 0.1}]
Length[%]

{{x -> 0.0102681}, {x -> 0.0106103}, {x -> 0.0109762}, {x ->
0.0113682}, {x -> 0.0117893}, {x -> 0.0122427}, {x ->
0.0127324}, {x -> 0.0132629}, {x -> 0.0138396}, {x ->
0.0144686}, {x -> 0.0151576}, {x -> 0.0159155}, {x ->
0.0167532}, {x -> 0.0176839}, {x -> 0.0187241}, {x ->
0.0198944}, {x -> 0.0212207}, {x -> 0.0227364}, {x ->
0.0244854}, {x -> 0.0265258}, {x -> 0.0289373}, {x ->
0.031831}, {x -> 0.0353678}, {x -> 0.0397887}, {x ->
0.0454728}, {x -> 0.0530516}, {x -> 0.063662}, {x -> 0.0795775}}

28

For the domain from 0.001 to 0.01 n runs from 32 to 318 and there are 287
roots.

RootSearch[Sin[1/x] == 0, {x, 0.001, 0.01}, InitialPrecision :> 40,
InitialSamples -> 5000, MaxBrentSteps -> 1000,
MaxSecantSteps -> 1000, PrecisionGoal :> 10];
Length[%]

287

For the domain from 0.0001 to 0.001 n runs from 319 to 3183 and there are
2865 roots. RootSearch took about 20 minutes on this. (I'm not certain if I
picked optimum values for the parameters to get the best efficiency.)

results3 =
RootSearch[Sin[1/x] == 0, {x, 0.0001, 0.001},
InitialPrecision :> 40, InitialSamples -> 30000,
MaxBrentSteps -> 1000, MaxSecantSteps -> 1000, PrecisionGoal :> 10];
Length[%]

2865

I would consider this pretty good performance.

From: Mark Adler [mailto:mad...@alumni.caltech.edu]


On 2010-08-28 23:51:42 -0700, Sam Takoy said:
> Is there a command for numerically finding all roots of a function in a
> given interval. If not, what's the best way of accomplishing this task?

NSolve and NRoots does that, but only for polynomials. I think you'd
have to roll your own to search for the roots in an interval of any
function. If you consider functions like Sin[1/x] from 0 to 0.1, you
can see that it might be quite difficult to do in general.

Mark

David Park

unread,
Sep 1, 2010, 6:27:17 AM9/1/10
to
Now that Daniel has shown the way here is a slightly reformulated form of
his routine that provides a single convenient command that should handle
most cases. The precision of the real variables in the calculation can be
specified. I don't know if the Reduce options would ever come into play
here, but I allowed for them. I also allowed control of messages from Reduce
with the default being Quiet because most of the time we wouldn't want to
see the warnings. Finally I allowed a wrapper function for Reduce for cases
such as the zeros of zeta example.

ClearAll[RealFunctionRoots];
RealFunctionRoots::usage =
"RealFunctionRoots[lhs \[Equal] rhs, {var, low, high}, precision, \
messages, wrapper] will return roots of the equation in the specified \
real domain.\nprecision specifies the precision of real variables to \
be used in the calculation and display, The default is \
MachinePrecision.\nmessages is a function, default Quiet, that \
controls messages from the routine.\nwrapper is a function, default \
Identity, that wraps the Reduce statement. It can be used to \
implement extra specification for the Reduce routine.\nOptions for \
Reduce may be passed.";

Options[RealFunctionRoots] = Options[Reduce];
SyntaxInformation[
RealFunctionRoots] = {"ArgumentsPattern" -> {_, {_, _, _}, _., _., \
_., OptionsPattern[]}};
RealFunctionRoots[equation_, {x_, low_, high_},
precision : (_Integer?Positive | MachinePrecision) :
MachinePrecision, messagewrapper : (_Function | _Symbol) : Quiet,
wrapperfunction_: Identity, opts : OptionsPattern[]] :=
Module[{work, equationp, lowp, highp},
{equationp, lowp, highp} = {equation, low, high} /.
r_Real :> SetPrecision[r, precision];
work = messagewrapper@
wrapperfunction@Reduce[{equationp, lowp <= x <= highp}, opts];
work = {ToRules[work]};
N[work, precision]]

Here are some of the examples.

f1[x_] := Sin[4 x] - (x + 1)/8

Plot[f1[x], {x, -13, 12}]
RealFunctionRoots[f1[x] == 0, {x, -13, 12}, 30]
f1[x] /. %

f2[x_] := If[x < 3/2, Sqrt[3/2 - x] - Exp[-4 x], 2 x Exp[-x]]
Plot[f2[x], {x, -25, 25},
PlotRange -> {{-0.5, 2.0}, {-2, 2}}]
RealFunctionRoots[f2[x] == 0, {x, -25, 25}, 30]
f2[x] /. %

f3[x_] := Exp[x - \[Pi]] - (1 - \[Pi] + x)
Plot[f3[x], {x, -3, 4},
PlotRange -> Automatic]
RealFunctionRoots[f3[x] == 0, {x, -3, 4}, 40]
f3[x] /. %

f4[x_] := (1.5 - Erf[x] - Erf[2 - x]) Exp[-Abs[x - 1]]

Plot[f4[x], {x, -400, 200},
PlotRange -> {{0, 2}, Automatic}]
RealFunctionRoots[f4[x] == 0, {x, -400, 200}, 40]
f4[x] /. %

f5[x_] :=
3.2 - 3 Cos[x] - Exp[-(x - 8 \[Pi])^2] - Exp[-(x - 14 \[Pi])^2]
Plot[f5[x], {x, 0, 30 \[Pi]}, PlotRange -> {{0, 30 \[Pi]}, Automatic}]
RealFunctionRoots[f5[x] == 0, {x, 0, 30 \[Pi]}, 30]
f5[x] /. %

f6[x_] := Piecewise[{{Sqrt[Abs[x] - 1.2], x <= 1.4}}, Indeterminate]

Plot[f6[x], {x, -8, 8}, PlotRange -> {Full, {-1, 1}},
MaxRecursion -> 4,
PlotPoints -> 25]
RealFunctionRoots[f6[x] == 0, {x, -8, 8}, 30]
f6[x] /. %

We might try the zeta function and find it involves a constant and doesn't
return the desired roots.

RealFunctionRoots[Zeta[1/2 + I y] == 0, {y, 0, 257}]

Then, if we were smart enough to think of Exist and Resolves, we could use
the wrapper version:

RealFunctionRoots[Zeta[1/2 + I y] == 0, {y, 0, 257}, 30, Quiet,
Resolve[Exists[C[1], #]] &]
Zeta[1/2 + I y] /. % // Chop[#, 10.^-27] &

Here is a type of equation that one gets in single reaction chemical
equilibrium equations, where y is the reaction extent and 0 and 1 are the
extent boundaries for positive reactant concentrations. Using FindRoot on
such functions is quite iffy because the solution is often very close to a
boundary point and there are often non-useful roots on the other side.

RealFunctionRoots[(4 y^2)/((2 - 2 y)^2 (2 - y)) == 100000, {y, 0,
1}, 30]
(4 y^2)/((2 - 2 y)^2 (2 - y)) == 100000 /. %

Just playing around of course. So don't worry Andrzej, us kiddies will stay
in the sandbox and out of the way of the real professionals.


From: Daniel Lichtblau [mailto:da...@wolfram.com]


David Park wrote:
> Install Ted Ersek's package RootSearch from:
>
> http://library.wolfram.com/infocenter/MathSource/4482/
>
> 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.
Ted's
> package is quite robust and returns all the roots in order.
>
>

> From: Sam Takoy [mailto:sam....@yahoo.com]
>
> Hi,
>

> Is there a command for numerically finding all roots of a function in a
> given interval. If not, what's the best way of accomplishing this task?
>

Andrzej Kozlowski

unread,
Sep 1, 2010, 6:28:32 AM9/1/10
to


Perhaps. But note that:

sols =
N[Reduce[Sin[1/x] == 0 && 10^-4 < x < 10^-3, x], 10]; // Timing

{2.88873,Null}

Length[List @@ sols]

2865

So what do you need RootSearch for?


Andrzej Kozlowski

David Park

unread,
Sep 1, 2010, 6:28:10 AM9/1/10
to
Three comments:

1) This appears to arise from new capabilities added to Root in Version 7,
capabilities that may have slipped under the radar screen of many users,
certainly under mine.

2) Nice work!

3) WRI can't provide a specific command for everything and often users will
have to add their own definitions, but 1-dimensional root searching of real
functions seems one place a conveniently packaged routine would be useful.

Andrzej Kozlowski

unread,
Sep 1, 2010, 6:28:43 AM9/1/10
to

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. Ted's
> package is quite robust and returns all the roots in order.


Mathematica's Reduce can do this (in the case of complex analytic functions, which covers most important practical cases) since version 7. I have posted examples of this several times, but somehow nobody seems to notice and 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, from the mathematical point of view, the methods used in the RootSearch package are rather primitive, relatively to the current state of knowledge in this field, of course. Not only does Reduce solve most transcendental equations much faster but it also does not require setting the values of any options 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, much more general methods of solving transcendental equations exist (in general the functions need not be analytic and can involve systems of n-real equations in n-variables with non-vanishing Jacobian). An example of such a method can be found 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 package into what is supposed to be a sophisticated mathematical program?
Sometimes certain capabilities are missing from Mathematica simply because implementing "bleeding edge" algorithms in a way that is satisfactory for a program that has the kind of aspirations that Mathematica obviously does takes time.

These sort of "helpful suggesting", particularly when coming from non-experts, have no chance whatever of being accepted and even can produce the impression of an ulterior motive being involved.

Gianluca Gorni

unread,
Sep 3, 2010, 6:10:14 AM9/3/10
to

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 12.28, Andrzej Kozlowski wrote:

>
> 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. Ted's
>> package is quite robust and returns all the roots in order.
>
>
> Mathematica's Reduce can do this (in the case of complex analytic functions, which covers most important practical cases) since version 7. I have posted examples of this several times, but somehow nobody seems to notice and 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, from the mathematical point of view, the methods used in the RootSearch package are rather primitive, relatively to the current state of knowledge in this field, of course. Not only does Reduce solve most transcendental equations much faster but it also does not require setting the values of any options 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, much more general methods of solving transcendental equations exist (in general the functions need not be analytic and can involve systems of n-real equations in n-variab

> l!


> es with non-vanishing Jacobian). An example of such a method can be found 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 package into what is supposed to be a sophisticated mathematical program?

> Sometimes certain capabilities are missing from Mathematica simply because implementing "bleeding edge" algorithms in a way that is satisfactory fora program that has the kind of aspirations that Mathematica obviously does takes time.
>
> These sort of "helpful suggesting", particularly when coming from non-exp=
erts, have no chance whatever of being accepted and even can produce the im=

Andrzej Kozlowski

unread,
Sep 4, 2010, 3:58:19 AM9/4/10
to

On 3 Sep 2010, at 12:10, Gianluca Gorni 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).

You seem to be using ====== instead of ==== which is a very basic error. In fact

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

x ==== Root[{Log[-((2*#1 - 1)/(2*#1^2 - 1))] + 2*#1 & ,
-0.86193624643066461859672257230652325787900736031\

58816695213`20.308604836334766}]

In[4]:== $Version

Out[4]== "7.0 for Mac OS X x86 (64-bit) (February 19, 2009)"

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

What are these ======== doing in your code? I don't think they could have been in it when you run it?

Reduce used exact methods so you have to rationalize the output or use equivalent approaches and it works fine in such cases. You are right that you can't use it with Interpolating functions since of course they are not analytic.

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


I concede that you make a reasonable case for this sort of capability. However, computations of this kind with Interpolating functions are generally quite unreliable, since the conditions required by FindRoot to work are often not satisfied. I don't think this sort of hit or miss approach is appropriate for a mathematical solver. Perhaps this sort of capabilities should be available in a package designed specifically for dealing with interpolating functions (which I don't consider as "mathematical objects").

Best regards

Andrzej

Matthias Bode

unread,
Sep 4, 2010, 3:58:51 AM9/4/10
to
Hola Gianluca:

In[1]:== Reduce[2*x + Log[-((-1 + 2*x)/(-1 + 2*x^2))] ==== 0 &&


-2 < x < -Sqrt[2]^(-1), x, Reals]

Out[1]== x ==== Root[{Log[-((-1 + 2*#1)/(-1 + 2*#1^2))] + 2*#1 & ,
-0.86193624643066461859672257230652325788`20.308604836334762}]

In[2]:== 2*x + Log[-((-1 + 2*x)/(-1 + 2*x^2))] /.
x -> -0.86193624643066461859672257230652325788`20.308604836334762

Out[2]== 0``19.380562970987107

== 0. * 10^ -20

Whereas (see your original equation below),

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

Syntax::sntxf:"2x+Log[-((-1+2x)/(-1+2x^2))]======" cannot be followed by "=
==0&&-2<x<-1/Sqrt[2]".

Syntax::tsntxi:"==0&&-2<x<-1/Sqrt[2]" is incomplete; more input is needed.

Syntax::sntxi:Incomplete expression; more input is needed.

Best regards,

MATTHIAS BODE
COCHABAMBA/BOLIVIA


> Date: Fri, 3 Sep 2010 06:10:23 -0400
> From: gianluc...@uniud.it
> Subject: Re: FindRoots?
> To: math...@smc.vnet.net


>
>
> 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 12.28, Andrzej Kozlowski wrote:
>
> >
> > On 30 Aug 2010, at 12:19, David Park wrote:
> >

> >> I don't know why Wolfram doesn't acquire the rights to this package an=
d
> >> incorporate it as part of regular Mathematica. Finding all the roots o=
f a
> >> real function on a 1-dimensional domain is the most common case of roo=
t
> >> finding, and it is the one thing regular Mathematica is very poor at. =


Ted's
> >> package is quite robust and returns all the roots in order.
> >
> >

> > Mathematica's Reduce can do this (in the case of complex analytic funct=
ions, which covers most important practical cases) since version 7. I have =
posted examples of this several times, but somehow nobody seems to notice a=
nd 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, f=
rom the mathematical point of view, the methods used in the RootSearch pack=
age are rather primitive, relatively to the current state of knowledge in t=
his field, of course. Not only does Reduce solve most transcendental equati=
ons much faster but it also does not require setting the values of any opti=
ons 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, m=
uch more general methods of solving transcendental equations exist (in gene=
ral the functions need not be analytic and can involve systems of n-real eq=
uations in n-varia!
> b
> > l!
> > es with non-vanishing Jacobian). An example of such a method can be fou=
nd in one of my demonstrations on the demonstrations site (the algorithm ca=


n easily be converted to a practical method for real life probl ems.
> >

> > In view of this, how could there be any justification for the suggestio=
n that Wolfram incorporates what is actually an impressive but amateurish p=


ackage into what is supposed to be a sophisticated mathematical program?

> > Sometimes certain capabilities are missing from Mathematica simply beca=
use implementing "bleeding edge" algorithms in a way that is satisfactory f=
ora program that has the kind of aspirations that Mathematica obviously doe=
s takes time.
> >
> > These sort of "helpful suggesting", particularly when coming from non-e=
xp==
> erts, have no chance whatever of being accepted and even can produce the =
im==


> pression of an ulterior motive being involved.
> >
>

--_637439e7-ac01-4858-84ae-ee894c0a0aa9_
Content-Type: text/html; charset="iso-8859-1"
Content-Transfer-Encoding: quoted-printable
X-Sun-Content-Length: 6413

<html>
<head>
<style><!--
.hmmessage P
{
margin:0px;
padding:0px
}
body.hmmessage
{
font-size: 10pt;
font-family:Tahoma
}
--></style>
</head>
<body class='hmmessage'>
Hola Gianluca:<br><br>In[1]:= Reduce[2*x + Log[-((-1 + 2*x)/(-1 + 2*x^2))=
] == 0 &amp;&amp; <br>&nbsp;&nbsp;&nbsp;&nbsp; -2 &lt; x =
&lt; -Sqrt[2]^(-1), x, Reals]<br><br>Out[1]= x == Root[{Log[-((=
-1 + 2*#1)/(-1 + 2*#1^2))] + 2*#1 &amp; , <br>&nbsp;&nbsp;&nbsp;&=
nbsp;&nbsp;&nbsp; -0.86193624643066461859672257230652325788`20.308604=
836334762}]<br><br>In[2]:= 2*x + Log[-((-1 + 2*x)/(-1 + 2*x^2))] /. <br>&=
nbsp;&nbsp; x -&gt; -0.86193624643066461859672257230652325788`20.3086=
04836334762<br><br>Out[2]= 0``19.380562970987107<br><br>= 0. * 10^ -20<=
br><br>Whereas (see your original equation below),<br><br>Reduce[2 x + Lo=
g[-((-1 + 2 x)/(-1 + 2 x^2))] === = <br>&nbsp;&nbsp; 0 &amp;&=
amp; -2 &lt; x &lt; -1/Sqrt[2], x, Reals]<br><br>Syntax::sntxf:"2=
x+Log[-((-1+2x)/(-1+2x^2))]===" cannot be followed by "=0&amp;&am=
p;-2&lt;x&lt;-1/Sqrt[2]".<br><br>Syntax::tsntxi:"=0&amp;&amp;-2=
&lt;x&lt;-1/Sqrt[2]" is incomplete; more input is needed.<br><br>Synt=
ax::sntxi:Incomplete expression; more input is needed.<br><br>Best regard=
s,<br><br><div><strong><em></em></strong><style>
</style><div><font face="Arial" size="2">MATTHIAS BODE<br>COCHABAMBA/BO=
LIVIA</font><br></div></div><br><br>&gt; Date: Fri, 3 Sep 2010 06:10:23=
-0400<br>&gt; From: gianluc...@uniud.it<br>&gt; Subject: [mg112177=
] Re: FindRoots?<br>&gt; To: math...@smc.vnet.net<br>&gt; <br>&gt; =
<br>&gt; In my opinion Reduce can replace RootSearch in some<br>&gt; ca=
ses but not in others.<br>&gt; <br>&gt; First of all, Reduce has bugs=
. Here is an analytic function<br>&gt; that clearly has a real root:<br>&=
gt; <br>&gt; Plot[2 x + Log[-((-1 + 2 x)/(-1 + 2 x^2))], {x, -2, =
-1/Sqrt[2]}]<br>&gt; <br>&gt; Still, Reduce does not see it (as of ve=
rsion 7.0.1):<br>&gt; <br>&gt; Reduce[2 x + Log[-((-1 + 2 x)/(-1 + 2 x^=
2))] ==== 0 &amp;&amp; -2 &lt;<br>&gt; x &lt; -1/Sqrt[=
2], x, Reals]<br>&gt; False<br>&gt; <br>&gt; (I reported this exa=
mple to wolfram last year).<br>&gt; <br>&gt; Next, Reduce has problem=
s with inexact input, and with<br>&gt; InterpolatingFunction, so that=
it won't work with<br>&gt; the output of NDSolve:<br>&gt; <br>&gt; s=
ol == x /.<br>&gt; First@NDSolve[{x'[t] ==== x[t] + 2, x[=
0] ==== -1}, x, {t, 0, 2}];<br>&gt; Reduce[sol[t] ==
=== 0 &amp;&amp; 0 &lt; t &lt; 2, t]<br>&gt; <br>&gt; o=
r, say, with functions obtained by interpolating between Locators.<br>&=
gt; RootSearch works fine in these cases.<br>&gt; <br>&gt; I have mad=
e some interactive panels where I can change the Locators<br>&gt; with th=
e mouse and I get in real time the roots of the interpolating<br>&gt; fun=
ction as big Points in the plot: I can do this with RootSearch,<br>&gt;=
but not with Reduce.<br>&gt; Unfortunately, I can't give these panels =
to other users, because<br>&gt; I can't assume that they have RootSearc=
h installed.<br>&gt; <br>&gt; I endorse the wish that the functionality=
of RootSearch were available<br>&gt; in the kernel.<br>&gt; <br>&gt;=
Best regards,<br>&gt; Gianluca Gorni<br>&gt; <br>&gt; <br>&gt; <=
br>&gt; On 01/set/2010, at 12.28, Andrzej Kozlowski wrote:<br>&gt; =
<br>&gt; &gt;<br>&gt; &gt; On 30 Aug 2010, at 12:19, David Park=
wrote:<br>&gt; &gt;<br>&gt; &gt;&gt; I don't know why Wolfram do=
esn't acquire the rights to this package and<br>&gt; &gt;&gt; incorpo=
rate it as part of regular Mathematica. Finding all the roots of a<br>&gt=
; &gt;&gt; real function on a 1-dimensional domain is the most common=
case of root<br>&gt; &gt;&gt; finding, and it is the one thing reg=
ular Mathematica is very poor at. Ted's<br>&gt; &gt;&gt; package is q=
uite robust and returns all the roots in order.<br>&gt; &gt;<br>&gt; =
&gt;<br>&gt; &gt; Mathematica's Reduce can do this (in the case of co=
mplex analytic functions, which covers most important practical cases) si=
nce version 7. I have posted examples of this several times, but somehow =
nobody seems to notice and we keep reading strange remarks like the above. =
Well, I don't want it to sound like I am trying to diminish Ted's achieve=
ment, but quite frankly, from the mathematical point of view, the met=
hods used in the RootSearch package are rather primitive, relatively to t=
he current state of knowledge in this field, of course. Not only does Red=
uce solve most transcendental equations much faster but it also does not re=
quire setting the values of any options by hand, and, what I think is m=
ost important, the results thus obtained are provably correct, which is=
not the case with Ted's package. In fact, much more general methods of s=
olving transcendental equations exist (in general the functions need not be=
analytic and can involve systems of n-real equations in n-varia!<br>&gt;=
b<br>&gt; &gt; l!<br>&gt; &gt; es with non-vanishing Jacobian). An=
example of such a method can be found in one of my demonstrations on the d=
emonstrations site (the algorithm can easily be converted to a practical me=
thod for real life probl ems.<br>&gt; &gt;<br>&gt; &gt; In view of =
this, how could there be any justification for the suggestion that Wolfra=
m incorporates what is actually an impressive but amateurish package into w=
hat is supposed to be a sophisticated mathematical program?<br>&gt; &gt=
; Sometimes certain capabilities are missing from Mathematica simply beca=
use implementing "bleeding edge" algorithms in a way that is satisfactory f=
ora program that has the kind of aspirations that Mathematica obviously doe=
s takes time.<br>&gt; &gt;<br>&gt; &gt; These sort of "helpful sugg=
esting", particularly when coming from non-exp=<br>&gt; erts, have =
no chance whatever of being accepted and even can produce the im=<br>&gt=
; pression of an ulterior motive being involved.<br>&gt; &gt;<br>&gt=
; <br> </body>
</html>=

--_637439e7-ac01-4858-84ae-ee894c0a0aa9_--

Syd Geraghty

unread,
Sep 4, 2010, 4:00:06 AM9/4/10
to
Gianluca,

Wolfram has evidently fixed the 1st problem you mentioned.

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

x==Root[{Log[-((-1+2 #1)/(-1+2 #1^2))]+2 #1&,-0.86193624643066461860}]

Yours truly ... Syd Geraghty

Syd Geraghty B.Sc., M.Sc.
sydge...@me.com
San Jose, CA

Mathematica 7.0.1.0 for Mac OS X x86 (64 - bit) (12 September 2009)
Licenses: L2983-5890, L3028-2592
MacOS X V 10.6.1 Snow Leopard
MacBook Pro 2.33 Ghz Intel Core 2 Duo 2GB RAM

On Sep 3, 2010, at 3:10 AM, Gianluca Gorni 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 <

ADL

unread,
Sep 4, 2010, 4:00:37 AM9/4/10
to
In fact, I am a little surprised by the suggestions related to using
Reduce for the purposes of soving numeric equations.

First, it has been repeated in this Group several times (and also by
Andrej) that doing so "can cause all sorts of weird problems".
Second, the impossibility to deal with interpolations with Reduce is
quite clear.
Third, Reduce has been reported a number of times to have strange
bugs, possibly due to its intrinsic complexity.
For example, please, try the following consecutive inputs on a fresh
kernel (I have 7.0.1)

(* x is confined to be negative and the result is correct *)
In[1]:= Reduce[2 x + Log[-((-1 + 2 x)/(-1 + 2 x^2))] == 0 && -2 < x <


-1/Sqrt[2], x, Reals]

Out[1]=
x == Root[{Log[-((-1 + 2*#1)/(-1 + 2*#1^2))] + 2*#1 & ,
-0.86193624643066461859672257230652325788`20.308604836334762}]


(* x can be negative and positive. The last root is lost *)
In[2]:= Reduce[{2 x + Log[-((-1 + 2 x)/(-1 + 2 x^2))] == 0, -2 < x <
2}, x, Reals]

Out[2]=
x == Root[{Log[-((-1 + 2*#1)/(-1 + 2*#1^2))] + 2*#1 & ,
0.56055055440035189796824605607850099728`20.301071037894882}]


(* This is the first input. The root is lost and will be never
recovered again *)
In[3]:= Reduce[2 x + Log[-((-1 + 2 x)/(-1 + 2 x^2))] == 0 && -2 < x <


-1/Sqrt[2], x, Reals]

Out[3]=
False


So, my understanding is that the concept of Rationalizing results and
applying N at the end to solve a numeric equation, which even Reduce
does internally, appears contrary to the expressed opinion of several
experts about Mathematica.

The introduction of an "NReduce" function (just to make a hypothesis)
could perhaps be more coherent.

ADL

Mark Adler

unread,
Sep 4, 2010, 4:01:42 AM9/4/10
to
On 2010-09-03 03:10:14 -0700, Gianluca Gorni said:
> 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 don't know what's with the quadruple equals, but if I put in a double
equals I get:

Reduce[2 x + Log[-((-1 + 2 x)/(-1 + 2 x^2))] == 0 && -2 < x <


-1/Sqrt[2], x, Reals]

x == Root[{Log[-((-1 + 2 #1)/(-1 + 2 #1^2))] + 2 #1 &,
-0.86193624643066461860}]

Which seems like a quite reasonable answer to an exact equation without
a closed form solution. The -0.8619... is the location of the root.

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

> Next, Reduce has problems with inexact input,

Indeed, Reduce will often complain when given numerical problems and say:

"Reduce was unable to solve the system with inexact coefficients or the
system obtained by direct rationalization of inexact numbers present in
the system. Since many of the methods used by Reduce require exact
input, providing Reduce with an exact version of the system may help."

Clearly Reduce was not designed to be a numerical equation solver, at
least so far. So the recommendations to use it as such are, I think, a
little off base.

Mark

Bob Hanlon

unread,
Sep 4, 2010, 4:01:53 AM9/4/10
to

$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 <gianluc...@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 12.28, Andrzej Kozlowski wrote:

>
> 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. Te=


d's
>> package is quite robust and returns all the roots in order.
>
>

> Mathematica's Reduce can do this (in the case of complex analytic functio=
ns, which covers most important practical cases) since version 7. I have po=
sted examples of this several times, but somehow nobody seems to notice and=
we keep reading strange remarks like the above. Well, I don't want it to s=
ound like I am trying to diminish Ted's achievement, but quite frankly, fro=
m the mathematical point of view, the methods used in the RootSearch packag=
e are rather primitive, relatively to the current state of knowledge in thi=
s field, of course. Not only does Reduce solve most transcendental equation=
s much faster but it also does not require setting the values of any option=
s by hand, and, what I think is most important, the results thus obtained a=
re provably correct, which is not the case with Ted's package. In fact, muc=
h more general methods of solving transcendental equations exist (in genera=
l the functions need not be analytic and can involve systems of n-real equa=
tions in n-varia!
!
b
> l!
> es with non-vanishing Jacobian). An example of such a method can be found=
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 pac=


kage into what is supposed to be a sophisticated mathematical program?

> Sometimes certain capabilities are missing from Mathematica simply becaus=
e implementing "bleeding edge" algorithms in a way that is satisfactory for=
a program that has the kind of aspirations that Mathematica obviously does =
takes time.
>
> These sort of "helpful suggesting", particularly when coming from non-exp=
==
erts, have no chance whatever of being accepted and even can produce the im=

Gianluca Gorni

unread,
Sep 4, 2010, 4:03:19 AM9/4/10
to

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

> 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 12.28, Andrzej Kozlowski wrote:
>
>>
>> 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==

Matthias Bode

unread,
Sep 5, 2010, 5:27:36 AM9/5/10
to
Hola:

In each re-posting of the original messages, cf. below, the Equal signs appear to increase in number according to two to the power n, n positive integer.

Here == I wrote one equal sign. Just ONE.

Best regards,

MATTHIAS BODE
COCHABAMBA/BOLIVIA

> Date: Sat, 4 Sep 2010 03:58:24 -0400
> From: akozl...@gmail.com


> Subject: Re: FindRoots?
> To: math...@smc.vnet.net
>
>

> On 3 Sep 2010, at 12:10, Gianluca Gorni 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).
>

> You seem to be using ============ instead of ======== which is a very bas=
ic error. In fact
>
> Reduce[2*x + Log[-((-1 + 2*x)/(-1 + 2*x^2))] ======== 0 &&


> -2 < x < -Sqrt[2]^(-1), x, Reals]
>

> x ======== Root[{Log[-((2*#1 - 1)/(2*#1^2 - 1))] + 2*#1 & ,
> -0.86193624643066461859672257230652325787900736031\
>
> 58816695213`20.308604836334766}]
>
> In[4]:==== $Version
>
> Out[4]==== "7.0 for Mac OS X x86 (64-bit) (February 19, 2009)"


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

> What are these ================ doing in your code? I don't think they co=


uld have been in it when you run it?
>

> Reduce used exact methods so you have to rationalize the output or use eq=
uivalent approaches and it works fine in such cases. You are right that you=
can't use it with Interpolating functions since of course they are not ana=
lytic.


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

> I concede that you make a reasonable case for this sort of capability. Ho=
wever, computations of this kind with Interpolating functions are generally=
quite unreliable, since the conditions required by FindRoot to work are of=
ten not satisfied. I don't think this sort of hit or miss approach is appro=
priate for a mathematical solver. Perhaps this sort of capabilities should =
be available in a package designed specifically for dealing with interpolat=

David Park

unread,
Sep 5, 2010, 5:26:28 AM9/5/10
to
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.


From: Gianluca Gorni [mailto:gianluc...@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 <gianluc...@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
>>

Ingolf Dahl

unread,
Sep 5, 2010, 5:29:01 AM9/5/10
to
Andrzej,
I must ask some question related to your answer. Better to ask and appear
unwise than not to ask and remain unwise.
The questions are interlaced below.

> -----Original Message-----
> From: Andrzej Kozlowski [mailto:akozl...@gmail.com]
> Sent: den 4 september 2010 09:58
> To: math...@smc.vnet.net
> Subject: Re: FindRoots?
>
> (snipped)


>
> Reduce used exact methods so you have to rationalize the output or use

equivalent
> approaches and it works fine in such cases. You are right that you can't
use it with
> Interpolating functions since of course they are not analytic.

I want to put a question mark on "of course". Say, an interpolation
polynomial of fifth degree - is that not "analytic"? Or a RBF (radial basis
function) interpolation, with a Gaussian radial basis function. And most
other interpolation methods (also those returned by NDSolve) are
differentiable any number of times everywhere except in isolated points or
along some lines. Then it should be a bookkeeping problem to analyze each
interval or region, one at a time? It should not be too difficult to
exchange the interpolation method, if that is preferred.

>
Gianluca Gorni wrote:
> > 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.

In some cases it should as an alternative be possible to interpolate the
inverse function to find roots. That requires in general scattered-point
interpolation, "available in a package designed specifically for dealing
with interpolating functions ", maybe found at
http://www.familydahl.se/mathematica . (It is possible that we have already
have solved the hard part of the problem when we have found that there is an
inverse function for some region.)

Andrzej Kozlowski wrote:
> I concede that you make a reasonable case for this sort of capability.
However,
> computations of this kind with Interpolating functions are generally quite
unreliable, since
> the conditions required by FindRoot to work are often not satisfied. I


don't think this sort

> of hit or miss approach is appropriate for a mathematical solver. Perhaps
this sort of
> capabilities should be available in a package designed specifically for
dealing with
> interpolating functions (which I don't consider as "mathematical
objects").

What are interpolating functions then? Please elaborate in detail, because I
am curious!

I think that different point sets in abstract vector spaces are important in
many places for our description of the world and of the universe of
thoughts. Interpolation functions make it possible to apply various
mathematical methods to these point sets. And I think that abstract point
set tools might become important in the further development of machine
intelligence.

Best regards

Ingolf Dahl


Andrzej Kozlowski

unread,
Sep 5, 2010, 7:08:47 AM9/5/10
to

On 5 Sep 2010, at 10:49, Ingolf Dahl wrote:

> Andrzej,
> I must ask some question related to your answer. Better to ask and appear
> unwise than not to ask and remain unwise.
> The questions are interlaced below.
>
>> -----Original Message-----
>> From: Andrzej Kozlowski [mailto:akozl...@gmail.com]
>> Sent: den 4 september 2010 09:58
>> To: math...@smc.vnet.net
>> Subject: Re: FindRoots?
>>
>> (snipped)
>>
>> Reduce used exact methods so you have to rationalize the output or use
> equivalent
>> approaches and it works fine in such cases. You are right that you can't
> use it with
>> Interpolating functions since of course they are not analytic.
>
> I want to put a question mark on "of course". Say, an interpolation
> polynomial of fifth degree - is that not "analytic"? Or a RBF (radial basis
> function) interpolation, with a Gaussian radial basis function. And most
> other interpolation methods (also those returned by NDSolve) are
> differentiable any number of times everywhere except in isolated points or
> along some lines. Then it should be a bookkeeping problem to analyze each
> interval or region, one at a time? It should not be too difficult to
> exchange the interpolation method, if that is preferred.

Of course interpolating polynomials are analytic but if you use Reduce with an interpolating polynomial you will certainly not encounter the problems Gianluca referred to. As for the other types of interpolating functions: remember that for Reduce to work it needs to be a complex analytic function, i.e. differentiable when considered as a function of a complex variable (i.e. the Cauchy-Riemann equations have to be satisfied), and of course that is not true for a non-polynomial interpolating function. Remember that the purpose of Reduce is to return mathematically "provable" results (and, in those cases when it can only produce a partial result it will display a message to the effect that the methods available to it are not sufficient to obtain the complete solution and either display a partial solution or none).


>
>
> Andrzej Kozlowski wrote:
>> I concede that you make a reasonable case for this sort of capability.
> However,
>> computations of this kind with Interpolating functions are generally quite
> unreliable, since
>> the conditions required by FindRoot to work are often not satisfied. I
> don't think this sort
>> of hit or miss approach is appropriate for a mathematical solver. Perhaps
> this sort of
>> capabilities should be available in a package designed specifically for
> dealing with
>> interpolating functions (which I don't consider as "mathematical
> objects").
>
> What are interpolating functions then? Please elaborate in detail, because I
> am curious!
>
> I think that different point sets in abstract vector spaces are important in
> many places for our description of the world and of the universe of
> thoughts. Interpolation functions make it possible to apply various
> mathematical methods to these point sets. And I think that abstract point
> set tools might become important in the further development of machine
> intelligence.

Well, I wrote that quickly and did not give it much thought. In fact, in a certain sense everything that Mathematica does can be given a rigorous mathematical description, and hence it is all "mathematical". However, some of the basic objects that Mathematica deals with not the kind of objects one deal with in "usual" mathematica yet in some sense they can be viewed as some sort of "simulations" of them. The main examples of such things are approximate numbers in Mathematica. One can certainly produce a rigorous mathematical theory of Mathematica's real or complex numbers but it will be in many ways different from what we mean by real and complex numbers in mathematics. As has been often mentioned on this forum, Mathematica's real numbers are really "fuzz balls", similar to intervals but not quite. The rules for working with them are perfectly mathematical but usually they are not known to the users (one can read about them in full detail in some technical publications but not in the Mathematica documentation). So what users do is regard them as some sort of approximations to certain standard mathematical concepts. In practice this is almost always enough. The problem occurs of course when you try to apply mathematical theorems, proven for the "real things" to these approximations. In general you can't expect to get reliable results.

The problem with (non-polynomial) interpolating functions is that they are constructed out of and applied to Mathematica's non-exact numbers. It seems clear that you can't expect to obtain mathematically valid results by applying to these entities theorems proved for the "standard" mathematica objects which they approximate of "simulate". For really reliable results you would need methods that fully took into account the actual "fuzzy" nature of these things. That is why functions like Reduce or Solve which rely on "standard mathematics" are not usually appropriate.

This is, of course, not an argument against having in Mathematica the sort of capabilities that Gianluca wrote about. Rather, I don't think the existing RootSearch actually offers them, since I think it relies on the built in FindRoot function. Although FindRoot is indeed intended for solving numerical equations I do not think it can be reliably used with InterpolatingFunction objects that arise as solutions of differential equations. I say this on the basis of my own experience. There is often a way to avoid this by using the "EventLocator" value of the option Method in NDSolve.


Andrzej Kozlowski


Andrzej Kozlowski

unread,
Sep 6, 2010, 4:13:21 AM9/6/10
to
It has been a long time since I looked at the source code of the RootSearchpackage (Ted sent it to me several years ago). I forgot the details of it and when writing this simply assumed it used FindRoot or an essentially equivalent method. In fact it implements Brent's method, which is of course one of the values of the Method option in FindRoot. I don't know exactly why Ted did wanted to implement this independently, but in any case it is not relevant to my comments. Everything "mathematical" that I wrote applies just as much to RootSearch as to FindRoot.

The point of my original post was not that Mathematica should not have a functionality to produce complete solutions of numerical equations but that leasing packages from users is not the right approach for Wolfram. The arguments for this view are unrelated to mathematics and I do not want to enter into them here (I know of only one exception that WRI has made to this rule: the Combinatorica package, co-authored by Steven Skiena and although there were very good reasons for doing so once I doubt that this is something WRI would wish to repeat).

From the "mathematical" viewpoint the problem with the RootSearch package is that it is unable to prove that all roots have been found. To know that the user has to rely on other ad hoc methods, like looking at the graph. These will of course fail whenever one is dealing with clusters of roots so close together than they cannot be distinguished by the eye.

The answers given by Reduce, on the other hand, are mathematically guaranteed to be correct (barring bugs, of course). It's pointless to compare an example where Reduce is trying to obtain a mathematical proof (which may take a very long time in some cases) with a method which will return "roots" which may not be roots at all.

Of course Reduce, as currently written, is limited only to complex analytic functions. But there exist methods which are much more general, and will work in all the cases where SearchRoot works, but, unlike SearchRoot, will return mathematically provable answers. One such method is due to Semenev and is the subject of two of my demonstrations on the Demosntrations site. It deals not just with a singe equation but with a system of n real equations in n variables. The functions are only required to be of class C2 - that is have continuous second derivatives. My implementation uses Interval and FindRoot, and the results are "guaranteed" (without the need for any options) in all cases where these two functions work. I don't think this includes the case of Interpolating functions because I don't think Interval arithmetic works in this case. But I don't believe that RootSearch is at all reliable in such cases - in fact I can provide real life important examples where it fails no matter what options are used.


Andrzej

On 5 Sep 2010, at 13:41, David Park wrote:

> I do not find any reference to or instance of FindRoot in the RootSearch
> code.

> From: Andrzej Kozlowski [mailto:akozl...@gmail.com]
>

David Park

unread,
Sep 6, 2010, 4:13:10 AM9/6/10
to

Andrzej Kozlowski

unread,
Sep 6, 2010, 4:42:28 AM9/6/10
to

On 5 Sep 2010, at 20:10, Gianluca Gorni wrote:

>
> On 03/set/2010, at 18.57, Andrzej Kozlowski wrote:
>
>>
>> On 3 Sep 2010, at 12:10, Gianluca Gorni wrote:
>>
>>>

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

>>> False
>>
>> You seem to be using === instead of == which is a very basic error. In fact
>>
>>
>> What are these ==== doing in your code? I don't think they could have been in it when you run it?
>
> Hi Andrzei!
>
> I had made copy-and-paste from Mathematica to Mail,
> and my original message had the correct double "==".
> The incorrect "====" must have crept in during
> transmission, I wonder where. How is this message
> showing up on your side?
>
> Best regards,
> Gianluca Gorni
>

I think the problem is caused by the MathGroup mail server. It seems
that when you send a post to the MathGroup with an equality sign the
same post is sent to other MathGroup users with the equality sign
doubled. In other words = is replaced by ==. This is very
inconvenient because it means that a lot of code when pasted into
Mathematica will turn into nonsense. The problem is clearly on the
server side and I don't think that we can do anything about it. I don't
know how long it has existed as I have not been very active in the
MathGroup recently and I only noticed it because of your post.

[I have checked the raw emails that came into my mailbox and I
saw ======== in the mail instead of ===. My server is not
doing anything to the posts as far as I can see. Note that
the posts come with "=" as hex 3D and other hex codes instead
of ASCII. This may be the problem with some server before
it gets to me - Moderator]


Andrzej

Mark McClure

unread,
Sep 7, 2010, 2:02:13 AM9/7/10
to
On Mon, Sep 6, 2010 at 4:13 AM, Andrzej Kozlowski <akozl...@gmail.com> wrote:

> The point of my original post was not that Mathematica should not have
> a functionality to produce complete solutions of numerical equations but
> that leasing packages from users is not the right approach for Wolfram.

> ...


> I know of only one exception that WRI has made to this rule: the
> Combinatorica package, co-authored by Steven Skiena

Actually, Mathematica contains quite a lot of code authored by users.
The palettes introduced in V7 were created by Eric Schultz of Walla
Walla Community College. Much of the geodetic capabilities were
authored by Thomas Mayer of The University of Connecticut. Stan Wagon
of Macalester tells me that he's contributed functionality to the
kernel. Sometimes, Wolfram even hires such authors; the author of the
IMTEK Mathematica supplement is such and example. If you click on
"Credits" in the About Mathematica box, you'll find a huge amount of
third party software used by Mathematica.

As far as this specific discussion goes, I must say that I find Reduce
to be rather mysterious. The documentation actually states that
"Reduce[expr,vars] reduces the statement...". It's clearly very
powerful but it's primarily algebraic. Thus it doesn't seem
appropriate in many situations - if a function is defined via a
numerical technique, for example.

Mark McClure

Murray Eisenberg

unread,
Sep 7, 2010, 2:02:56 AM9/7/10
to
Perhaps that's the case with add-on packages, literally. But my
understanding is that WRI has licensed some technology from others. In
fact, I recall an explicit statement at a Mathematica conference that
they saw no reason to "reinvent the wheel" and were therefore going to
use (I think this is what the item is:) a 3rd party floating-point
library from then on instead of continuing to "roll their own".


On 9/6/2010 4:13 AM, Andrzej Kozlowski wrote:

> The point of my original post was not...


> but that leasing packages from users is not the right approach for Wolfram.

> The arguments for this view are unrelated to mathematics and I do not want

> to enter into them here (I know of only one exception that WRI has made to
> this rule: the Combinatorica package...)

--
Murray Eisenberg mur...@math.umass.edu
Mathematics & Statistics Dept.
Lederle Graduate Research Tower phone 413 549-1020 (H)
University of Massachusetts 413 545-2859 (W)
710 North Pleasant Street fax 413 545-1801
Amherst, MA 01003-9305

Andrzej Kozlowski

unread,
Sep 7, 2010, 6:08:18 AM9/7/10
to

On 7 Sep 2010, at 06:03, Mark McClure wrote:

> On Mon, Sep 6, 2010 at 4:13 AM, Andrzej Kozlowski <akozl...@gmail.com> wrote:
>
>> The point of my original post was not that Mathematica should not have

>> a functionality to produce complete solutions of numerical equations but


>> that leasing packages from users is not the right approach for Wolfram.

>> ...


>> I know of only one exception that WRI has made to this rule: the

>> Combinatorica package, co-authored by Steven Skiena
>
> Actually, Mathematica contains quite a lot of code authored by users.
> The palettes introduced in V7 were created by Eric Schultz of Walla
> Walla Community College. Much of the geodetic capabilities were
> authored by Thomas Mayer of The University of Connecticut. Stan Wagon
> of Macalester tells me that he's contributed functionality to the
> kernel. Sometimes, Wolfram even hires such authors; the author of the
> IMTEK Mathematica supplement is such and example. If you click on
> "Credits" in the About Mathematica box, you'll find a huge amount of
> third party software used by Mathematica.

Well, as far as I can tell, none of these concern not "packages" that are already available for download from the Internet (free or not). I think such cases are very different and normally the authors should take full responsibility for them themselves. There may be of course exceptions (arguably Combinatorica was one of them) but they should be v. rare.


>
> As far as this specific discussion goes, I must say that I find Reduce
> to be rather mysterious. The documentation actually states that
> "Reduce[expr,vars] reduces the statement...". It's clearly very
> powerful but it's primarily algebraic. Thus it doesn't seem
> appropriate in many situations - if a function is defined via a
> numerical technique, for example.
>

It really very much depends what you mean by "algebraic". If you think that all branches of mathematics where you actually prove things are "algebraic" than I agree. The methods Reduce uses to solve non-polynomial equations are baed on theorems in complex analysis, such as Rouche's theorem. Also, if you evaluate this

Reduce[AiryAi[E^x] ==== E^-x + 1 && Abs[x] < 2, x]

you will see a message referring to something that is not normally considered "algebraic".

Andrzej

Ingolf Dahl

unread,
Sep 7, 2010, 6:07:14 AM9/7/10
to
Andrzej states that a non-polynomial interpolating function cannot be
analytic, and therefore cannot be used together with Reduce. I had asked
about RBF interpolation with a Gaussian radial basis function. As I
understand it, such an interpolation is a sum of Gaussians, centered around
different points, and this sum can be written in such a way that it is
analytic (The exponent of the Gaussian should then contain x squared, and
not x times the conjugate of x.) Such interpolation functions cannot in
general be handled by Reduce, but this is not related to being analytic or
not, they are simply too complicated.

Reduce can handle the following function with simple coefficients and three
terms, if given some time:

Reduce[0.3 E^(-(1/2) (-2.23 + x)^2) - 0.7 E^(-(1/2) (-1.45 + x)^2) +
1.1 E^(-(1/2) (0.` + x)^2) == 0, {x}, Reals]

but does not seem to get ready with this, not within my patience:

Reduce[-16.367079553136787` E^(-(1/2) (-2.33` + x)^2) +
20.460124007820784` E^(-(1/2) (-2.23` + x)^2) -
5.524386343053091` E^(-(1/2) (-1.45` + x)^2) +
2.4125732736976717` E^(-(1/2) (0.` + x)^2) == 0, {x}, Reals]

, both expressions equally analytic. The roots of this second equation are
not in any way problematic, and are robust to the "fuzziness" of the real
numbers in Mathematica. But to me it seems as if Reduce does inadequate
things in this case.
Maybe I am generalizing from too few examples. Please correct me if that is
the case!
But anyway, I consider it as snobbery to say that Reduce is superior to
RootSearch. They do not play in the same division, and perform different
tasks.

Best regards

Ingolf Dahl

> -----Original Message-----
> From: Andrzej Kozlowski [mailto:akozl...@gmail.com]

> Sent: den 5 september 2010 13:09
> To: math...@smc.vnet.net
> Subject: Re: FindRoots?
>
>

> On 5 Sep 2010, at 10:49, Ingolf Dahl wrote:
>
> > Andrzej,
> > I must ask some question related to your answer. Better to ask and
appear
> > unwise than not to ask and remain unwise.
> > The questions are interlaced below.
> >
> >> -----Original Message-----

> >> From: Andrzej Kozlowski [mailto:akozl...@gmail.com]

> no!


> t in the Mathematica documentation). So what users do is regard them as
some sort of
> approximations to certain standard mathematical concepts. In practice this
is almost always
> enough. The problem occurs of course when you try to apply mathematical
theorems, proven
> for the "real things" to these approximations. In general you can't expect
to get reliable
> results.
>
> The problem with (non-polynomial) interpolating functions is that they are
constructed out
> of and applied to Mathematica's non-exact numbers. It seems clear that you
can't expect to
> obtain mathematically valid results by applying to these entities theorems
proved for the
> "standard" mathematica objects which they approximate of "simulate". For
really reliable
> results you would need methods that fully took into account the actual
"fuzzy" nature of
> these things. That is why functions like Reduce or Solve which rely on
"standard
> mathematics" are not usually appropriate.
>

Andrzej Kozlowski

unread,
Sep 8, 2010, 12:57:47 AM9/8/10
to

On 7 Sep 2010, at 10:02, Ingolf Dahl wrote:

> But anyway, I consider it as snobbery to say that Reduce is superior to
> RootSearch. They do not play in the same division, and perform different
> tasks.

Why should it matter to me what you call it? Is that an argument?

If you read what I have written on this subject you may noticed that there
are ways to do essentially the same things that RootSearch attempts to do and obtain the complete set of roots (provably complete) (I am not referring to using Reduce). That, I do call "superior" and really I don't care if anyone thinks it's "snobbery" or not.

Andrzej Kozlowski

Andrzej Kozlowski

unread,
Sep 8, 2010, 12:59:15 AM9/8/10
to


In case my meaning is still unclear, I will once again repeat something that I have written already on this forum several times on various occasions.

If Wolfram really wanted a method that solves numerical transcendental equations that could implement one of several existing sophisticated modern algorithms, including the one due to Semenev (2007) that is the subject of my demonstrations:

http://demonstrations.wolfram.com/SolvingSystemsOfTranscendentalEquations/

and

http://demonstrations.wolfram.com/SemenovsAlgorithmForSolvingSystemsOfNonlinearEquations/

Unlike the SearchRoot package this algorithm finds the complete set of solutions of any system of n equations in n unknowns satisfying certain mild conditions. I have no problem at all calling it "superior" to what is implemented in the SearchRoot package. Moreover, it can solve systems of equations (yes, systems) which Mathematica is at present is unable even to attempt (one such example is given in the first demonstration above) and neither of course can SearchRoot. The set of solutions is provably complete and the user does not need to set any options.

My implementation of Semenev's algorithm is probably inefficient as it was only meant to demonstrate the method and not to be used in practice. Semenev's paper gives accurate complexity estimates for the method which show that it is quite practical if the number of equations (and unknowns) is not too large.

It would be trivial for Wolfram to implement this algorithm and it would not require leasing anything from anyone, with all the problems that this usual entails.

Andrzej Kozlowski

Mark McClure

unread,
Sep 8, 2010, 1:00:19 AM9/8/10
to
On Tue, Sep 7, 2010 at 3:26 AM, Andrzej Kozlowski <akozl...@gmail.com> wrote:

> Well, as far as I can tell, none of these concern not "packages" that are
> already available for download from the Internet (free or not). I think such
> cases are very different and normally the authors should take full
> responsibility for them themselves. There may be of course exceptions
> (arguably Combinatorica was one of them) but they should be v. rare.

I'm not sure what your point is, since I can't really tell what your
first sentence means. I can say that all the specific examples I
referred to began life as independent packages (available either on
the internet or as book supplements) and have subsequently become part
of the Mathematica proper.

I'm not particularly arguing that the RootSearch package should be
licensed by Wolfram Research. I'm simply stating the fact that there
is plenty of precedent for such a move.


> It really very much depends what you mean by "algebraic".

I mean that, when given exact input, it returns exact results, rather
than numerical approximations. Even if a numerical integral is used
to determine the number of zeros (as in your example), the result is
exact. (Perhaps, I should have said exact. :) )

Mark

Andrzej Kozlowski

unread,
Sep 9, 2010, 4:22:28 AM9/9/10
to

On 8 Sep 2010, at 19:08, Richard Fateman wrote:

> Just a brief question...
>
> Does this method work only for zero-dimensional solutions (i.e. points)?
> In which case it is far less interesting than Reduce, it seems to me.
>
>
> RJF
>

Yes, indeed. It works only for zero dimensional solutions and, of
course, it is much less interesting than Reduce. But still, it does
something that Reduce does not do and for which there seems to be a
demand.

Andrzej Kozlowski

Ingolf Dahl

unread,
Sep 9, 2010, 5:31:04 AM9/9/10
to
To Andrzej Kozlowski and MathGroup,
Yes, I was too harsh in using the word "snobbery". But I feel that it is not
good to let the best (the Semenov method?) be an enemy of the good
(RootSearch), when the Semenov method not yet is fully implemented. You have
argued for the capabilities of Reduce to find zeroes, so I wonder of course
if a routine based on the Semenov method would fail in a similar way. So
much better if it does not. Maybe we should ask MathGroup, if there is
anyone willing to implement your demonstrations as a more general package?
The code does not look too difficult, even if I not immediate gasp all
details (I usually don't). And as usual, the hell might be in the details.

However, in order to implement the Semenov method and to obtain provable
results, you state that there are "certain mild conditions" that must be
fulfilled. These conditions seem to be that the program needs to obtain
explicit and valid upper bounds for the first and second derivatives. The
task to find these bounds cannot be automated in a provable way in the
general case, so one problem has effectively been replaced by another. One
could put the responsibility to solve this second problem on the user, but I
cannot see such a solution as superior. Nevertheless, there seems to be
niches where the Semenov method might be very useful, especially in the
multidimensional case, so it should be worth the effort to create a full
implementation. I doubt that such an implementation will be option-free -
would it not be nice to have an "Automatic" way of estimating the bounds of
the derivatives? But then some roots might be missed.

Best regards

Ingolf Dahl

> -----Original Message-----
> From: Andrzej Kozlowski [mailto:akozl...@gmail.com]
> Sent: den 8 september 2010 06:59
> To: math...@smc.vnet.net
> Subject: Re: FindRoots?
>
>

Andrzej Kozlowski

unread,
Sep 9, 2010, 5:31:26 AM9/9/10
to

On 9 Sep 2010, at 09:50, Ingolf Dahl wrote:

> To Andrzej Kozlowski and MathGroup,
>

> However, in order to implement the Semenov method and to obtain provable
> results, you state that there are "certain mild conditions" that must be
> fulfilled. These conditions seem to be that the program needs to obtain
> explicit and valid upper bounds for the first and second derivatives. The
> task to find these bounds cannot be automated in a provable way in the
> general case, so one problem has effectively been replaced by another. One
> could put the responsibility to solve this second problem on the user, but I
> cannot see such a solution as superior.

This tasks *is* automated in Mathematica. The function Interval does it for you.

Andrzej Kozlowski

0 new messages