Message Boards Message Boards

14 Replies
12 Total Likes
View groups...
Share this post:

how to get exact solutions (error free solution)

Posted 9 years ago

I have a small problem.

sk = List[](*Here an empty list 'sk' is generated*)

PTT = x^900 + 4 x^100 - 
  5(*Here I have defined a polynomial "PTT" in x having degree 900*)

sk = x /. NSolve[PTT == 0, x];
(*     Here "PTT" is solved for 900 roots.
And all the soltuins of polynomial "PTT" are stored in the list sk   *)

For[i = 1, i < 11, i = i + 1,
 x = sk[[i]];
 Print[{x, PTT}]

(*This loop Substitute the first 10 solutions from the array "sk" in the \
polynomial "PTT"*)

{-1.00218-0.0203695 I,-3.78364*10^-13+2.10498*10^-13 I}

{-1.00218+0.0203695 I,-3.78364*10^-13-2.10498*10^-13 I}

{-1.00218-0.0277838 I,-2.82441*10^-13+1.42109*10^-13 I}

{-1.00218+0.0277838 I,-2.82441*10^-13-1.42109*10^-13 I}

{-1.00194-0.0351982 I,-3.73035*10^-13-7.32747*10^-14 I}

{-1.00194+0.0351982 I,-3.73035*10^-13+7.32747*10^-14 I}

{-1.0019-0.0129806 I,-2.75779*10^-13+2.85105*10^-13 I}

{-1.0019+0.0129806 I,-2.75779*10^-13-2.85105*10^-13 I}

{-1.00148-0.0425983 I,2.35367*10^-13+1.20792*10^-13 I}

{-1.00148+0.0425983 I,2.35367*10^-13-1.20792*10^-13 I}

When solutions of the polynomial are substutituted in equation, they must completly satisfy the polynomial.

But As Observed in the output, it is not so.

How can this error be removed.
I need exact solutions, which can satisfy the orignal equation

Mathematica code of this file can be found in attachments.

POSTED BY: abhishek sharma
14 Replies

If x1, x2 are solution of a polynomial f(x),

then f(x1) and f(x2) must vanish.. i.e. f(x1)=0 and f(x2)=0. But in above example it is not so...

POSTED BY: abhishek sharma

Probably you should check some of the standard references on numerical analysis. Also, I've seen nothing thus far to indicate that the NSolve results are inconsistent with what one would expect (by "one", I mean "one with at least a passing familiarity with numerical approximation and truncation/roundoff error").

For a very quick example of what I refer to, let's find a root of a simple quadratic polynomial.

quadpoly = x^2 - 7;
rt = FindRoot[quadpoly == 0, {x, 3}]

(* Out[209]= {x -> 2.64575131106} *)

Now let's check that the residual is zero, that is to say, plugging the root into the polynomial actually makes it vanish.

In[210]:= quadpoly /. rt

(* Out[210]= 8.881784197*10^-16 *)

That's not an exact zero. Does that mean the root is somehow wrong?

There is another serious issue with your example. It is one of scaling. You have machine numbers and also numbers with exponents too large in magnitude to be represented as machine numbers. And furthermore the variation in the scales is far larger than the ULP of a machine number (look that up if it is not a familiar term). This is a situation where using either exact or high precision input is fairly important. Without doing that, basically you will be dealing with massive truncation error, to the point where you won't be able to tell if the results obtained are at all valid.

It would be a good idea to consider carefully what it is you are trying to solve, and what will be required in order to do so. Then consider what it is that would be desirable properties of a solution (absolute error, relative error, some feature of residuals, ...). It makes little sense to have a problem that is not carefully set up, obtain solutions that seem to give good (as in small) residuals, and then state that they are no good; you have to begin with a set up for which sensible results might be derived. And you have to understand what "sensible" might mean, in the context of numerical approximation.

POSTED BY: Daniel Lichtblau

Thanks sir,

Sir, u rightly ponied the variation in the scales in-between different coefficients of the polynomial. I am working on it.

kindly answer one wore thing. I have following one example.

aw = 23; bw = 34; cw = 12; nc = (aw/bw) + 10/cw

How can I set Precision to evaluate "nc"

POSTED BY: abhishek sharma

Your code exactly evaluates "nc" already: 77/51 is an exact rational number. It has no exact, finite decimal representation: 77/51 is the best you can have. If you want to see lots of digits, try something like N[nc,1000].

POSTED BY: John Doty

Consider an exact polynomial expression:

poly = -6 + 17 x - 32 x^2 + 26 x^3 - 20 x^4 + 14 x^5 - 8 x^6 + 3 x^7

Note that it contains no approximate numbers like 1.4: all of the coefficients are exact rational numbers (integers in this case). In Mathematica, any number containing a decimal point is an approximate number, not an exact number. You cannot expect exact results from approximate input.

Make an equation, solve it:

exact = Solve[poly == 0]
{{x -> 1/3 (1 - I Sqrt[2])}, {x -> 1/3 (1 + I Sqrt[2])}, {x -> 
   Root[-6 + 5 #1 - 4 #1^2 + 3 #1^3 - 2 #1^4 + #1^5 &, 1]}, {x -> 
   Root[-6 + 5 #1 - 4 #1^2 + 3 #1^3 - 2 #1^4 + #1^5 &, 2]}, {x -> 
   Root[-6 + 5 #1 - 4 #1^2 + 3 #1^3 - 2 #1^4 + #1^5 &, 3]}, {x -> 
   Root[-6 + 5 #1 - 4 #1^2 + 3 #1^3 - 2 #1^4 + #1^5 &, 4]}, {x -> 
   Root[-6 + 5 #1 - 4 #1^2 + 3 #1^3 - 2 #1^4 + #1^5 &, 5]}}

I get seven roots, all exact. Two are expressed as radicals, five have no such representation and are expressed using Root. These are the exact solutions, as we can verify:

Simplify[poly /. exact]
{0, 0, 0, 0, 0, 0, 0}

However, if I want numerical results:

approx = NSolve[poly == 0]
{{x -> -0.551685 - 1.25335 I}, {x -> -0.551685 + 1.25335 I}, {x -> 
   0.333333 - 0.471405 I}, {x -> 0.333333 + 0.471405 I}, {x -> 
   0.805786 - 1.2229 I}, {x -> 0.805786 + 1.2229 I}, {x -> 1.4918}}

The numerical solutions are necessarily approximate: exact numerical solutions of this equation are not possible. Most polynomial equations do not have exact numerical solutions: the solutions are irrational and thus have no finite numerical representation. Substituting the approximate solutions back into the original polynomial yields results that are approximately, but not exactly, zero:

poly /. approx
{1.84741*10^-13 + 1.91847*10^-13 I, 1.84741*10^-13 - 1.91847*10^-13 I,
  6.66134*10^-16 + 1.9984*10^-15 I, 
 6.66134*10^-16 - 1.9984*10^-15 I, -1.42109*10^-14 + 
  7.10543*10^-15 I, -1.42109*10^-14 - 
  7.10543*10^-15 I, -2.13163*10^-14}

And that is just as you should expect.

POSTED BY: John Doty

Thanks sir for your precious time.

POSTED BY: abhishek sharma

This is not making sense. You have an approximate set of values and plug them in to compute residuals. How can those residuals possibly be "exactly" zero?

This is a matter of numerical mathematics, not Mathematica.

POSTED BY: Daniel Lichtblau

Sincere thanks for your valuable time.

Here in attachments I am uploading my original problem. I hope Instructions in file attached are clear to understand.

POSTED BY: abhishek sharma

If you want exact results, you must use exact input. You are using approximate numbers like 0.51 in your input. Any number expressed with a decimal point is approximate in Mathematica. 0.51 (approximate) is not the same thing as 51/100 (exact, rational). If you don't actually want exact results, but rather results with better precision than you're getting, you can try higher WorkingPrecision as Daniel suggests.

POSTED BY: John Doty

What do you mean by "working well"? For high-order polynomials with exact coefficients, Solve generally yields results in terms of Root objects with exact coefficients. This is the best that Solve can do, but it's actually better than you may realize. Such objects are precisely determined algebraic numbers. Mathematica can operate on them as exact numbers, and can determine their decimal value to any practical precision.

POSTED BY: John Doty

Try NSolve with higher than default WorkingPrecision. And be prepared to wait.

POSTED BY: Daniel Lichtblau

If your problem is complex, Solve will take a long time. That's unavoidable.

POSTED BY: Frank Kampas

yes it worked for this polynomial...

But in my real problem the polynomial in more complex... Neither "Solve" nor "NSolve" are working well for that...

kindly suggest some another method...

If u wish I can upload my orignal equation

POSTED BY: abhishek sharma

NSolve is a numerical algorithm. For exact results, use Solve.

POSTED BY: Frank Kampas
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract