[Date Index] [Thread Index] [Author Index]
Re: How accurate is the solution for high degree algebraic equation?
On 24 Oct 2012, at 09:32, Alexandra <watanabe.junzo at gmail.com> wrote: > I wanted to know all the solutions of f = (-z - 1)^d - (-z^d - 1)==0, where d=54. > I did the following: > > d = 54; f = (-z - 1)^d - (-z^d - 1); > sol = NSolve[f == 0,z]; > a = z /. sol; > > So a is a set of solutions. > > If I compute > f /. z -> a[] // N > It returns a number very close to zero. This is natural. > > But if I compute > f /. (z -> a[]) // N > > Then > Mathematica returns > 12.0047 + 14.7528 I > > I cannot say a[] is a solution of f=0. > > Many other elements in the solution set a does not seem to satisfy the equation. > Only the last few terms in a are satisfactory enough as solutions. > > Is the degree too high? > > You are using machine precision to do the computation. This will become unreliable for high degree equations. However, in a case like this there is no need to rely on machine precision at all (this is Mathematica, after all, not Excel or some hand held calculator.) First, Mathematica can find *exact* solutions: d=54;f=(-z-1)^d-(-z^d-1); sol=Solve[f==0,z]; a=z/.sol; Now, you can compute any of them to an arbitrary degree of precision: s=N[a[],100] = -0.0151922469877919406332569754104769863293567828645926442905752792132370242954156332994391135662242200+0.1736481776669303488517166267693147960003758733221472022345073025847854427518548781450433122117761692 I You can also verify your solution: f/.N[sol[],100]//Chop 0 This is so fast these days that in this particular case there is no need for another method. But in more difficult cases it may be faster to use NSolve with higher working precision: d = 54; f = (-z - 1)^d - (-z^d - 1); sol = NSolve[f == 0, z, WorkingPrecision -> 100]; a = z /. sol; Note, however, that NSolve orders the roots in a different way that Solve so a[] will also be a solution (with precision 100) but it will not be the same solution (although the set of solutions is the same, of course): a[] = -0.0420104876845111255626252330431971078102751675023564614914520080003610201692407202063501699086031192+0.2868032327110902531032801731669967109003063896092812580697966433037586988411218301430417590833776560 I Again verifying f/.sol[]//Chop 0 Without the Chop you will get very tiny residues in the above answers. Andrzej Kozlowski