Re: [Q] accuracy control in eq. solving?
- To: mathgroup at smc.vnet.net
- Subject: [mg23716] Re: [mg23678] [Q] accuracy control in eq. solving?
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Mon, 5 Jun 2000 01:09:12 -0400 (EDT)
- References: <200005291624.MAA17934@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
James wrote:
>
> Hello!
>
> I have a question regarding to the accuracy control of
> the equation solution.
>
> I'm currently solving polynomial equations(degrees are at least 5),
> and wondering if I can control the accuracy level
> of the solution.
> Specifically, when I recomputed the original equation
> by replacing the answers that mathematica gives,
> I got the error of roughly 10^-7 or 10^-8 for degree 9.
> Regarding this, I have the following questions.
>
> (1) Is there any way to control(increase) the accuracy?
> (2) What would be the best known math tool (in terms of accuracy)
> in solving high-degree polynimial equations?
> (3) I'm a newbie in mathematica,
> so I manually type in, when I check the answer(polynomial eqs.).
> Is there any systax for that in mathematica?
>
> Thank you.
>
It is hard to respond without a specific example because there may be
features peculiar to your equations. For example, are leading or
trailing coefficients large or small compared to the rest? Do you
typically encounter multiple roots? These issues can affect accuracy.
Here are some ideas for question 1. I'll demonstrate with random
equations of degree 9, where coefficients are 6 digits or so.
In[39]:= randomEqn[deg_, var_] := Table[Random[Integer, {-10^6,10^6}],
{deg+1}] . var^Range[0,deg]
In[40]:= InputForm[eq9 = randomEqn[9,x]]
Out[40]//InputForm=
582120 - 943992*x - 119444*x^2 + 873998*x^3 + 906859*x^4 + 692184*x^5 -
270059*x^6 - 28395*x^7 - 288276*x^8 + 80551*x^9
Here we give a machine precision solution (the default behavior of
NSolve/NRoots). This uses numeric root-finding in machine arithmetic.
In[41]:= solmach = NSolve[eq9, x];
We'll check the residuals.
In[42]:= InputForm[eq9 /. solmach]
Out[42]//InputForm=
{-7.421476766467094*^-10, -2.6193447411060333*^-10 + 0.*I,
-2.6193447411060333*^-10 + 0.*I, 3.259629011154175*^-9 -
3.725290298461914*^-9*I, 3.259629011154175*^-9 +
3.725290298461914*^-9*I,
2.5295321393059567*^-11 + 1.1459633242338896*^-10*I,
2.5295321393059567*^-11 - 1.1459633242338896*^-10*I, 0., 0.}
Not bad, around 10^(-10).
Here is a way to work with 50 digits of precision. We can assess
accuracy after the fact by again checking residuals.
In[43]:= solbig = NSolve[eq9, x, WorkingPrecision->50];
In[44]:= InputForm[eq9 /. solbig]
Out[44]//InputForm=
{2.002213`-8.0872*^-51, 3.206964`-8.6281*^-51 + 6.26678`-9.3372*^-52*I,
3.206964`-8.6281*^-51 + -6.26678`-9.3372*^-52*I,
5.32312`-8.6198*^-51 + -1.254737`-9.2474*^-51*I,
5.32312`-8.6198*^-51 + 1.254737`-9.2474*^-51*I,
1.5`-12.6168*^-56 + 2.6`-12.3835*^-56*I, 1.5`-12.6168*^-56 +
-2.6`-12.3835*^-56*I, 1.04`-9.1174*^-50, 1.95`-11.3376*^-50}
This tells us we have accuracy to around 40 or so digits (a bit more for
some solutions, a bit less for others).
I do not have an answer for question 2. Mathematica is using an
implementation of Jenkins-Traub polynomial root-finding. One alternative
might be to form a companion matrix and extract numeric eigenvalues. I
do not know whether this is likely to be more or less accurate. I
believe both methods have trouble at machine precision when faced with
multiple roots.
I'm not exactly certain what you are asking in question 3 but possibly
the replacements used above to check residuals is along the lines you
seek.
Daniel Lichtblau
Wolfram Research