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

about SOLVE

58 views
Skip to first unread message

Mohamad

unread,
Feb 25, 2005, 10:11:36 AM2/25/05
to
Hi every one,
I'm trying to obtain the solution of the equation:
z^20-exp(0.6*z+0.25*z^14+0.15*z^38-1=0;
Normally this equation has 20 roots inside the circle unit, but when
trying to use solve like the following:
syms z;
sol=solve(z^20-exp(0.6*z+0.25*z^14+0.15*z^38-1));
I obtain only one solution!! (the obvious root 1)
Can somebody help me?
Thank you in advance

dev...@math.udel.edu

unread,
Feb 25, 2005, 2:54:35 PM2/25/05
to
Mohamad wrote:
> I'm trying to obtain the solution of the equation:
> z^20-exp(0.6*z+0.25*z^14+0.15*z^38-1=0;
> Normally this equation has 20 roots inside the circle unit

What do you mean by "normally"?

> but when trying to use solve like the following:
> syms z;
> sol=solve(z^20-exp(0.6*z+0.25*z^14+0.15*z^38-1));

>From the syntax, I'll assume that you are using Matlab, which is
processing the solve through Maple. The following Maple code will find
the 20 roots inside the unit disk in about 2 minutes. I don't know the
exact syntax that you would use to process this Maple code in Matlab,
but I am pretty sure that it can be done.

.> Digits:= 14:
.> f:= z-> z^20-exp(6/10*z+1/4*z^14+15/100*z^38-1):
.> randNum:= rand(-2^30..2^30)/2^30:
.> randC:= ()-> evalf(randNum())+I*evalf(randNum()):
.> RootsFound:= {}: InDisk:= {}: failures:= 0: Failures:= 0:
.> do
.> r:= fsolve(
.> f(z)
.> ,{z= randC()}
.> ,complex
.> ,fulldigits
.> ,{z= -1-I..1+I}
.> ,avoid= RootsFound
.> );
.> if r=NULL
.> or r::specfunc(anything, fsolve)
.> or r in RootsFound
.> then
.> failures:= failures+1; Failures:= Failures+1;
.> if failures<2*nops(RootsFound) then next else
.> print("Giving up.");
.> break
.> fi
.> fi;
.> failures:= 0;
.> RootsFound:= RootsFound union r;
.> if abs(rhs(r[])) < 1.0 then
.> InDisk:= InDisk union r;
.> if nops(InDisk) >= 20 then break fi
.> fi
.> od:

.> failures, Failures, nops(InDisk), nops(RootsFound);

0, 22, 20, 25

.> RootsFound:= map(rhs, [RootsFound[]]);
.> InDisk:= map(rhs, [InDisk[]]);

RootsFound := [-0.88290480457775 + 0.27493030452949 I,
0.93023697804313 + 0.30061772249567 I,
0.96056333566004 + 0.81305565253490 I,
0.77216590847663 - 0.58400116661888 I,
0.63797763569112 + 0.90678424520243 I,
-0.56790125927015 - 0.74802543767007 I,
-0.024217404552146 + 0.94473405164932 I,
-0.88290480457775 - 0.27493030452949 I,
-0.75717314939774 - 0.53486583854938 I, -0.92967788113916,
1.0000000000000, 0.27936144072755 + 0.92143429733504 I,
-0.31876744421805 - 0.88531533906831 I,
0.27936144072755 - 0.92143429733504 I,
0.77216590847663 + 0.58400116661888 I,
-0.75717314939774 + 0.53486583854938 I,
-0.88290480457776 - 0.27493030452949 I,
-0.56790125927015 + 0.74802543767007 I,
0.58639974860150 + 0.89183223210159 I,
0.54834302036656 + 0.80381795254952 I,
0.54834302036656 - 0.80381795254952 I,
-0.024217404552146 - 0.94473405164932 I,
0.93023697804313 - 0.30061772249567 I,
-0.31876744421805 + 0.88531533906831 I,
0.58623534369954 + 0.97271088175519 I]

InDisk := [-0.88290480457775 + 0.27493030452949 I,
0.93023697804313 + 0.30061772249567 I,
0.77216590847663 - 0.58400116661888 I,
-0.56790125927015 - 0.74802543767007 I,
-0.024217404552146 + 0.94473405164932 I,
-0.88290480457775 - 0.27493030452949 I,
-0.75717314939774 - 0.53486583854938 I, -0.92967788113916,
0.27936144072755 + 0.92143429733504 I,
-0.31876744421805 - 0.88531533906831 I,
0.27936144072755 - 0.92143429733504 I,
0.77216590847663 + 0.58400116661888 I,
-0.75717314939774 + 0.53486583854938 I,
-0.88290480457776 - 0.27493030452949 I,
-0.56790125927015 + 0.74802543767007 I,
0.54834302036656 + 0.80381795254952 I,
0.54834302036656 - 0.80381795254952 I,
-0.024217404552146 - 0.94473405164932 I,
0.93023697804313 - 0.30061772249567 I,
-0.31876744421805 + 0.88531533906831 I]

Daniel Lichtblau

unread,
Feb 25, 2005, 3:00:52 PM2/25/05
to


Here are a couple of methods, coded in Mathematica. Both rely on the
function being analytic in a neighborhood of the region of interest
(the unit circle, in this case). If there were poles one could make
suitable adaptions with a bit of work.

The first approximates the function by a power series centered at the
origin, finds roots, then uses them as start points for a local root
finder.

In[1]:= expr = z^20-Exp[3/5*z+1/4*z^14+3/20*z^38-1];

In[2]:= ss = Normal[Series[expr, {z,0,50}]];

In[3]:= solns = Select[z /. {ToRules[NRoots[ss==0, z]]}, Abs[#]<1.05&];

In[4]:= InputForm[z/.Map[FindRoot[expr==0, {z,#},
WorkingPrecision->20]&, solns]]

Out[4]//InputForm=
{-0.929677881139164716562180936666819`20,
-0.882904804577755691681970227447832`20.1304 -
0.274930304529491056694246991095135`19.6237*I,
-0.882904804577755691681970227447832`20.1304 +
0.274930304529491056694246991095135`19.6237*I,
-0.757173149397737376799789065781946`20.0626 -
0.534865838549377151710863994078075`19.9117*I,
-0.757173149397737376799789065781946`20.0626 +
0.534865838549377151710863994078075`19.9117*I,
-0.567901259270137599223088874022957`19.932 -
0.74802543767005600536040769074953`20.0517*I,
-0.567901259270137599223088874022957`19.932 +
0.74802543767005600536040769074953`20.0517*I,
-0.318767444218031585425155502426662`19.6804 -
0.885315339068292687998676155660066`20.124*I,
-0.318767444218031585425155502426662`19.6804 +
0.885315339068292687998676155660066`20.124*I,
-0.024217404552093898129442902483429`18.5592 -
0.944734051649391938493972457293357`20.1504*I,
-0.024217404552093898129442902483429`18.5592 +
0.944734051649391938493972457293357`20.1504*I,
0.27936144072754617825343079997138`19.6131 -
0.921434297335035155373925397166638`20.1314*I,
0.27936144072754617825343079997138`19.6131 +
0.921434297335035155373925397166638`20.1314*I,
0.548343020366562508426129314846016`19.9014 -
0.803817952549523775797730086958279`20.0675*I,
0.548343020366562508426129314846016`19.9014 +
0.803817952549523775797730086958279`20.0675*I,
0.77216590847663063114324120810951`20.0523 -
0.584001166618883966672266436506201`19.931*I,
0.77216590847663063114324120810951`20.0523 +
0.584001166618883966672266436506201`19.931*I,
0.930236978043133062019727191986406`20.1289 -
0.300617722495673253943254682316175`19.6384*I,
0.930236978043133062019727191986406`20.1289 +
0.300617722495673253943254682316175`19.6384*I,
1.000000000000000000664805169598839`20}

This is probably the better approach in general. If the series
approximation required too many terms to be feasible as a polynomial
approximation in the entire range of interest you could break it apart
and use smaller series centered in different places, then merge all
solutions obtained in this manner.

A second approach is based on the Cauchy integral formula (CIF). We can
count roots in sectors (I will use rectangles), and isolate by
bisecting the sectors in each direction until they contain at most one
root. We can then evaluate the root inside, again by CIF. We use this
value as a starting point for a high precision local root finder.

This method uses a few different "epsilons". The first one (called eps
below, and set to 1/1000), gives how far beyond the rectangle
boundaries we will proceed in the quadrature. This is useful to avoid
the case where a root lies on a segment over which we wish to integrate
(for example, this will happen with the two real roots if one boundary
of our rectangle is right on the real axis).

Another, set to 1/100, is an approximation to how far outside the unit
circle we will go when we "partition" it into rectangles. The third is
a tolerance for deciding when two roots are to be regarded as the same.

eps = 10^(-3);

rectangleIntegrate[func_, a_, b_, c_, d_, z_, f2_:1] := Module[
{t1=0, t2=0, t3=0, t4=0, integrand=f2*D[func,z]/func},
t1 = NIntegrate[Evaluate[integrand], {z,a+I*c,b+I*c}];
t2 = NIntegrate[Evaluate[integrand], {z,b+I*c,b+I*d}];
t3 = NIntegrate[Evaluate[integrand], {z,b+I*d,a+I*d}];
t4 = NIntegrate[Evaluate[integrand], {z,a+I*d,a+I*c}];
(t1+t2+t3+t4)/(2*Pi*I)]

isolateRootsInRectangle[func_, a_, b_, c_, d_, z_, eps_] := Module[
{numroots, avgx=(a+b)/2, avgy=(c+d)/2, r1, r2, r3, r4, len=0,
allrts},
numroots =
Round[rectangleIntegrate[func,a-eps,b+eps,c-eps,d+eps,z]];
If[numroots==0, Return[{}]];
If[numroots==1, Return[{{a,b,c,d}}]];
r1 = isolateRootsInRectangle[func,a,avgx,c,avgy,z,eps];
r2 = isolateRootsInRectangle[func,a,avgx,c+avgy,d,z,eps];
r3 = isolateRootsInRectangle[func,a+avgx,b,c,avgy,z,eps];
r4 = isolateRootsInRectangle[func,a+avgx,b,c+avgy,d,z,eps];
Join[{r1,r2,r3,r4}]
]

refineRoots[func_, intervals_, z_, eps_]:=
Map[rectangleIntegrate[
func,#[[1]]-eps,#[[2]]+eps,#[[3]]-eps,#[[4]]+eps,z,z]&, intervals]

rootsInUnitCircle[func_, z_, tol_:10^(-5)] := Module[
{rts, r1, r2, r3, r4, r5, r6, r7, r8, eps=1/1000, x=1,
excess=1/100},
rts = isolateRootsInRectangle[func,-1/Sqrt[2],1/Sqrt[2],
-1/Sqrt[2],1/Sqrt[2],z,eps];
While[ArcCos[x]<Pi/4,
r1 = isolateRootsInRectangle[func,1/Sqrt[2],x,
Sin[ArcCos[x]],Sin[ArcCos[x-excess]],z,eps];
r2 = isolateRootsInRectangle[func,1/Sqrt[2],x,
-Sin[ArcCos[x-excess]],-Sin[ArcCos[x]],z,eps];
r3 = isolateRootsInRectangle[func,-x,-1/Sqrt[2],
Sin[ArcCos[x]],Sin[ArcCos[x-excess]],z,eps];
r4 = isolateRootsInRectangle[func,-x,-1/Sqrt[2],
-Sin[ArcCos[x-excess]],-Sin[ArcCos[x]],z,eps];
r5 = isolateRootsInRectangle[func,
Sin[ArcCos[x]],Sin[ArcCos[x-excess]],1/Sqrt[2],x,z,eps];
r6 = isolateRootsInRectangle[func,
-Sin[ArcCos[x-excess]],-Sin[ArcCos[x]],1/Sqrt[2],x,z,eps];
r7 = isolateRootsInRectangle[func,
Sin[ArcCos[x]],Sin[ArcCos[x-excess]],-x,-1/Sqrt[2],z,eps];
r8 = isolateRootsInRectangle[func,
-Sin[ArcCos[x-excess]],-Sin[ArcCos[x]],-x,-1/Sqrt[2],z,eps];
rts = Join[rts,r1,r2,r3,r4,r5,r6,r7,r8];
x -= excess;
];
rts = refineRoots[func, rts, z, eps];
rts = Map[
FindRoot[Evaluate[func]==0, {z,#}, WorkingPrecision->20]&, rts];
rts = z /. rts;
rts = Union[rts, SameTest->(Abs[#1-#2]<tol&)];
Select[rts, Abs[#]<=1&]
]

For the example in question,

InputForm[rootsInUnitCircle[expr, z]]

will give, after numerous messages about numeric integration
difficulties, the following result.

Out[26]//InputForm=
{-0.929677881139161812840909328858515`20.1505 + -1.9`-11.542*^-32*I,
-0.882904804577754980481792420591482`20.1304 +
0.274930304529490867913307754617456`19.6237*I,
-0.882904804577754980481792420591269`20.1304 -
0.274930304529490867913307754504865`19.6237*I,
-0.757173149397738341054699838305194`20.0626 -
0.534865838549377040197541416789759`19.9117*I,
-0.757173149397738341054699838305194`20.0626 +
0.534865838549377040197541416789759`19.9117*I,
-0.567901259270149215952731567846963`19.932 -
0.748025437670068092410036316911108`20.0517*I,
-0.567901259270149215952731567846963`19.932 +
0.748025437670068092410036316911108`20.0517*I,
-0.31876744421804628882956346413738`19.6804 -
0.885315339068314097443391391163738`20.124*I,
-0.31876744421804628882956346413738`19.6804 +
0.885315339068314097443391391163738`20.124*I,
-0.024217404552145731181306270583393`18.5592 +
0.944734051649321729733624338861818`20.1504*I,
-0.024217404552145731181306270413747`18.5592 -
0.944734051649321729733624338858926`20.1504*I,
0.279361440727543971441802663968991`19.6131 -
0.921434297335036031806331752726211`20.1314*I,
0.279361440727543971441802663968991`19.6131 +
0.921434297335036031806331752726211`20.1314*I,
0.548343020366562508436597927402393`19.9014 -
0.803817952549523775796628805828892`20.0675*I,
0.548343020366562508436597927402393`19.9014 +
0.803817952549523775796628805828892`20.0675*I,
0.772165908476630631142931317757618`20.0523 -
0.584001166618883966672076437695646`19.931*I,
0.772165908476630631142931317937326`20.0523 +
0.584001166618883966672076438156072`19.931*I,
0.930236978043133061796819282913374`20.1289 -
0.300617722495673254002463738296876`19.6384*I,
0.930236978043133061796819282913374`20.1289 +
0.300617722495673254002463738296876`19.6384*I,
1.`20.1505 + 3.91461932421`0*^-21*I}

One can use quadrature settings that will remove the messages in this
case at the expense of much slower evaluation (and the same final
result). For example, the settings below suffice.

SetOptions[NIntegrate, SingularityDepth->10, MinRecursion->4,
MaxRecursion->16, WorkingPrecision->25]

This method has its appeal but it is likely to fall prey to numeric
problems when roots cluster or otherwise cause trouble for the
quadrature. It has the advantage that one need not be concerned with
how large a series approximation to take, or whether and where to
change the centering of the series.


Daniel Lichtblau
Wolfram Research

A N Niel

unread,
Feb 26, 2005, 8:57:48 AM2/26/05
to
> > syms z;
> > sol=solve(z^20-exp(0.6*z+0.25*z^14+0.15*z^38-1));
> > I obtain only one solution!! (the obvious root 1)

"solve" tries to find symbolic solutions. If you want numerical
solutions, try something else.

Did you try RootFinding[Analytic] ?

dev...@math.udel.edu

unread,
Feb 26, 2005, 3:26:02 PM2/26/05
to
A N Niel wrote:
> "solve" tries to find symbolic solutions. If you want numerical
> solutions, try something else.
>
> Did you try RootFinding[Analytic] ?

I tried RootFinding[Analytic], and after several hours I still have no
answers. I also tried an approach similar to Daniel Litchblau's second
approach -- the one that counts the roots in a region by integrating
around the boundary -- but Maple seems unable to handle any of the
numeric integrations. The approach of using fsolve with random initial
points or using a series approximation to get the initial points works
in about a minute.

Alec Mihailovs

unread,
Feb 26, 2005, 5:17:40 PM2/26/05
to
> The approach of using fsolve with random initial
> points or using a series approximation to get the initial points works
> in about a minute.

That is a good approach in general. For this particular example, one can
just use roots of 1 to find 20 roots in the disc,

> tt:=time():Digits:=14:
> s:=[seq(evalf(exp(I*k*Pi/10)),k=0..19)]:
> sol:=[seq(fsolve(z->z^20-exp(6/10*z+1/4*z^14+15/100*z^38-1),z,complex)
> ,z in s)]:
> sort(sol,(x,y)->Re(x)<Re(y));

[-0.92967788113916, -0.88290480457775 - 0.27493030452949 I,

-0.88290480457775 + 0.27493030452949 I,

-0.75717314939774 - 0.53486583854938 I,

-0.75717314939774 + 0.53486583854938 I,

-0.56790125927015 - 0.74802543767007 I,

-0.56790125927015 + 0.74802543767007 I,

-0.31876744421805 - 0.88531533906831 I,

-0.31876744421805 + 0.88531533906831 I,

-0.024217404552146 - 0.94473405164932 I,

-0.024217404552146 + 0.94473405164932 I,

0.27936144072755 - 0.92143429733504 I,

0.27936144072755 + 0.92143429733504 I,

0.54834302036656 - 0.80381795254952 I,

0.54834302036656 + 0.80381795254952 I,

0.77216590847663 - 0.58400116661888 I,

0.77216590847663 + 0.58400116661888 I,

0.93023697804313 - 0.30061772249567 I,

0.93023697804313 + 0.30061772249567 I, 1.]

> time()-tt;

0.266

> nops(sol);

20

Alec


Carl Devore

unread,
Feb 26, 2005, 8:47:03 PM2/26/05
to
On Fri, 25 Feb 2005, Daniel Lichtblau wrote:
> This is probably the better approach in general. If the series
> approximation required too many terms to be feasible as a polynomial
> approximation in the entire range of interest you could break it apart
> and use smaller series centered in different places, then merge all
> solutions obtained in this manner.
>
> A second approach is based on the Cauchy integral formula (CIF). We can
> count roots in sectors (I will use rectangles), and isolate by
> bisecting the sectors in each direction until they contain at most one
> root. We can then evaluate the root inside, again by CIF. We use this
> value as a starting point for a high precision local root finder.

The benefit of the second method is that it counts the number of roots.
In general, you can't be sure that the first method is finished because
you don't know the number of roots. In this particular problem, it is
easy to apply Rouche's Theorem to prove that there are 20 roots in and on
the unit circle, but that won't work in general.

A hybrid of the two methods -- much faster than the second and more robust
than the first -- is to count the roots in large areas with the integral
method, and then uses the series method to get that many initial values.
This will work well if there are not multiple roots are nearly equal
roots, which is the usual case. And even if there are multiple or nearly
equal roots, this hybrid method will find most of the roots quickly, and
you only need to use the more expensive methods for a few that have
already been fairly closely located.

Daniel Lichtblau

unread,
Feb 28, 2005, 11:40:45 AM2/28/05
to


In addition to the various methods shown thus far, you might do this
via homotopy. This will work because you know in advance that the
expression is a perturbation of z^20-1/E.

We set up the homotopy in the variable t so that at t=0 we have the
simple polynomial and at t=1 we have the full expression. For this to
work well in Mathematica I found it expedient to separate into explicit
real and imaginary parts. This might not be necessary using a different
nondefault method or a different program (by setting SolveDelayed->True
below I am invoking a DAE solver behind the scenes).

expr = z^20 - Exp[3/5*z+1/4*z^14+3/20*z^38-1];

htpy[t_] := z^20 - (1-t)*Exp[-1] - t*Exp[3/5*z+1/4*z^14+3/20*z^38-1]

{re,im} = ComplexExpand[{Re[htpy[t]],Im[htpy[t]]}, z] /.
{Re[z]->x[t],Im[z]->y[t]}

We find the roots to the simple polynomial and separate into real and
imaginary parts. These will be our starting values.

rts = Apply[List, Map[Last, NRoots[z^20==1/E, z]]];
rts = Map[{Re[#],Im[#]}&, rts];

The defining condition for our homotopy is that at each value of time
t, we have a solution to expr above. Hence the derivative vanishes.

ode[{x0_,y0_}] := {D[re,t]==0, D[im,t]==0, x[0]==x0, y[0]==y0}

Now for each start value we tope our way over to the solution at t=1.

desolns = {x[t],y[t]} /.
Map[NDSolve[ode[#], {x[t],y[t]}, {t,0,1}, SolveDelayed->True]&, rts];

We evaluate the solutions at t=1 and convert to explicit complex
values.

In[146]:= InputForm[solns = Map[#[[1]]+I*#[[2]]&, Map[First, desolns /.
t->1]]]

Out[146]//InputForm=
{-0.9296777638034137 + 0.*I, -0.882904773967465 -
0.27493015468435966*I,
-0.882904773967465 + 0.27493015468435966*I,
-0.7571731919487987 - 0.534865809706014*I, -0.7571731919487987 +
0.534865809706014*I, -0.5679012692271274 - 0.748025473331044*I,
-0.5679012692271274 + 0.748025473331044*I, -0.31876740911815493 -
0.8853153319249604*I, -0.31876740911815493 + 0.8853153319249604*I,
-0.024217395962507344 - 0.9447340488047146*I,
-0.024217395962507344 + 0.9447340488047146*I,
0.279361427623412 - 0.9214342550330425*I,
0.279361427623412 + 0.9214342550330425*I, 0.5483429835284926 -
0.8038179339079018*I, 0.5483429835284926 + 0.8038179339079018*I,
0.7721659184486211 - 0.5840011515292556*I,
0.7721659184486211 + 0.5840011515292556*I,
0.930236972846258 - 0.3006177067291801*I,
0.930236972846258 + 0.3006177067291801*I, 1.000000544225191 + 0.*I}

These are just machine arithmetic approximations and moreover NDSolve
uses automatic PrecisionGoal equal to half the the WorkingPrecision. In
this case we used machine arithmetic so the solutions are onlyy
expected to hold to about 8 digits or so. This is of course more than
adequate to use as initial values for a local root finder to
arbitrarily high precision.


Daniel Lichtblau
Wolfram Research

0 new messages