Re: FindRoots?
- To: mathgroup at smc.vnet.net
- Subject: [mg112129] Re: FindRoots?
- From: "David Park" <djmpark at comcast.net>
- Date: Wed, 1 Sep 2010 06:27:26 -0400 (EDT)
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.
David Park
djmpark at comcast.net
http://home.comcast.net/~djmpark/
From: Daniel Lichtblau [mailto:danl at 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.
>
>
> David Park
> djmpark at comcast.net
> http://home.comcast.net/~djmpark/
>
>
>
> From: Sam Takoy [mailto:sam.takoy at 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