On Aug 14, 3:37 pm, New2Sympy <[email protected]> wrote:
> Hi again,
>
> I am getting confused now. See the script and the results below.
>
You are running into numerical issues (underflow, I belileve) because
of scaling issues. If you don't normalize your polynomial then the the
coefficients of the polynomial are on the order of 1e-15. Since your
roots are on the order of 1 that means that if you substitute 1 into
your polynomial you will get something very close to zero:
>>> mypoly=myEqConstants.as_poly(r)
>>> mypoly.subs(r,1)
1.85794794291198e-16
So if you are using a numerical solver and you say that if the
function gets closer than 1e-15 to the x-axis that that represents a
root for you. Well...any function multiplied by a small enough factor
will give you something close to zero! :-) Let sympy help you get this
into normalized form:
>>> mymonic=mypoly.as_monic();mymonic
Poly(1.0*r**6 + 2.32*r**5 - 126.6544*r**4 + 118.147448015123*r**2, r)
>>> mymonic.subs(r,1)
-5.18695198487708
OHH...so it's not so close to zero after all. Your root numerical
solver is going to do a lot better for you now:
first variant with monic:
Solutions using nsolve function in the range are:
-12.4399173136662
0.978377565046911
-0.960891724232247
10.1024314728516
> Second Variant (using Poly) kind of works, but there is a small fix
> required. "mypoly.coeffs" ommits the zeros, therefore if the
> coefficients for x^0, x^1 and x^3 are 0, they would not appear. Is
> this a bug or a feature?
Feature...it only gives the coeffs of terms that are there, not those
that aren't. That way if you are working with x**100 - 1 you don't get
a lot of zeros. It is incovenient to type these in by hand for mpmath,
however, so use sympy like this to do it for you:
>>> #get the dictionary of monomials and coefficients
>>> cod=mymonic.as_dict()
>>> cod
{(2,): 118.147448015123, (5,): 2.32000000000000, (6,):
1.00000000000000, (4,): -126.654400000000}
Now there must be a better way to do this, but not knowing it we can
construct a list of monomials from 6 to 0 asking d to give them to us.
But what do we do if a monomial is not there? Use the ability to have
the dictionary give you a default value:
>>> cod = [d.get((j,),0) for j in reversed(range(0,po.LM[0]+1))]
^ ^ ^
| | +- that's the
leading monomial, (6,)
| |
| _we need to get them in 6 to 0
order
|
get that monomial, but if it's not there, get 0
>>> cod
[1.00000000000000, 2.32000000000000, -126.654400000000, 0,
118.147448015123, 0, 0]
>>>
> Third Variant (using solve) is completely messed now, solutions don't
> get even close to the real solutions.
>
feeding the "basic" form of the MONIC poly (an expression rather than
a Poly) into solve with quartics=True takes a long time but finally
gives:
>>> ans = solve(mymonic.as_basic(), r, simplified=False, quartics=True)
>>>
> I plotted the function in a spreadsheet, and found that the solutions
> provided by the Second Variant (with the manual fix) are the correct
> ones.
>
I plotted it, too, but without the monic form I just got a graph that
looked like it didn't go through the zero axis at all. This
underscores the importance of working with a monic form of the
polynomial.
> Any ideas on how to fix all this to get the same results from all
> methods? how can I believe any of the results otherwise?
>
Every method is doing the best it can:
1) you gave it something that was so close to zero everywhere,
especially in the minima region around 0 that it looked like there
were several possibilities that were within the computational limit
and within your specified tolerance. Give it something that doesn't
look so close to zero and you will get results consistent with your
expectation: give it the monic form.
2) poly (I'm guessing) probably did some normalizing for you and found
the roots as you expected.
3) the current quartic routine which would be needed, I believe, to
solve your equation, is broken and is awaiting an update. It's doing
the best it can :-) Using the version that is awaiting review and
patching into the current sympy (by others who are also doing the best
they can with a hundred+ chnages to review and add into sympy) gives
the following results:
-12.4569900815341 - 0.0126318166896176*I
-0.719450735477597 + 0.239542505387883*I
0
0.72799186899517 - 0.246077177232919*I
10.1284489480165 + 0.0191664885346538*I
Those imaginary parts are too large for my liking and I'm looking into
why that is so, but I have a couple of issues that I am tracking down
right now that might be the solution.
An introductory numerical methods book would be a great help to you if
you are interested in numerical procedures...Burden and Faires
(Numerical Methods) is one I used, but a basic Chemical Engineering
first-year book for numerical methods would be great, too:
At http://tinyurl.com/kt486x there are some amazon books:
I taught from Numerical Methods for Engineers by Steven C. Chapra and
Raymond P. Canale and this book had lots of sample code and
discussions of a wide variety of basic numerical procedures. They
teach about the pitfalls like you are running into.
Burden and Faires is on the list, too.
An analysis of numerical methods is described at http://tinyurl.com/p53bad
and that might have some good discussions about why methods fail (and
I didn't make the url end in 'bad'!). If you understand what the
methods are trying to do and how they do it, you can help them do
their best, using your intelligence until the artificial intelligence
of the CAS catches up. But even that raises an issue...if you say you
want to solve an equation close to the computational representation of
the computer, should the CAS not allow you to do it? Maybe you really
*do* want to see how it fails. So that's why sometime there are
higher-level function calls (like solve) and lower-level calls (like
to mpmath) to do the work for you.
Best regards,
/c
Christopher Smith
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to [email protected]
For more options, visit this group at http://groups.google.com/group/sympy?hl=en
-~----------~----~----~----~------~----~------~--~---
