On Aug 12, 4:48 pm, Anartz Unamuno <[email protected]> wrote:
> Below is the outcome I get when I use myEqConstants.n(100):
>
> third variant:
> [-0.192261149873352 + .0e-21⋅ⅈ, -29.4436689512432 + .0e-19⋅ⅈ,
> 27.1235615040264 - .0e-19⋅ⅈ, 0, 0.192368597090154 - .0e-21
> ⋅ⅈ]
>
> the imaginary parts are zeros, but they are still in the solution. Should
> they not disappear?
When you have a solution with a quartic, there are parts of the
solution that have I*sqrt(3) in them. If you are careful (I suppose)
you should be able to get those to cancel out and end up with the real
answer only. When you start evaluating, though, when that term is
encountered, it adds in to the rest of the calculation as an imaginary
component (with some uncertainty) and you end up with that "almost
zero" artifact when you are done. Try the "chop=True" argument which
works most of the time:
>>> for s in mySols:
... print N(s,n=4,chop=True)
...
27.12
-1.923e-1
1.924e-1
-29.44
0
As for the symbolic approach, that seems good (especially if you just
want to change a few parameters) but I would recommend that even
though sympy *can* handle symbols, you do a little work to make it not
have to handle so many...at least until sympy is modified to know when
to ignore clusters of constants. Let's use the simpler manipulations
to look at your eq:
>>> myeq=(-4*a**2*r**2*e**2*(0.002*c*x*(o + v) + 2*a*b*r*e + 2*a*s*r*e + 2*a*w*r
*e + a*e*r**2)**2-4*a**2*r**2*e**2*(0.002*c*x*(o + v) - (2*a*b*r*e +
2*a*s*r*e +
2*a*w*r*e + a*e*r**2))**2+16*a**4*d**2*r**4*e**4 +
16*a**4*r**4*f**2*e**4 )
OK, let's change it to a poly so we can see the order clearly.
Although we could use expand, there's a nice method to pull off
coefficients:
>>> Poly(myeq,r)
Poly(-8*a**4*e**4*r**6 + (-32*b*a**4*e**4 - 32*s*a**4*e**4 -
32*w*a**4*e**4)*r**
5 + (-32*a**4*e**4*s**2 - 32*a**4*e**4*w**2 + 16*a**4*e**4*f**2 -
32*a**4*b**2*e
**4 + 16*a**4*d**2*e**4 - 64*b*s*a**4*e**4 - 64*b*w*a**4*e**4 -
64*s*w*a**4*e**4
)*r**4 + (-3.2e-5*a**2*c**2*e**2*o**2*x**2 -
3.2e-5*a**2*c**2*e**2*v**2*x**2 - 6
.4e-5*o*v*a**2*c**2*e**2*x**2)*r**2, r)
Now let's look at the coefficients in "monic" form (reducing the
leading term's coeff to 1):
>>> for xi in list(mypoly.as_monic().iter_all_coeffs()):
... print xi
...
1
4*b + 4*s + 4*w
8*b*s + 8*b*w + 8*s*w - 2*d**2 - 2*f**2 + 4*b**2 + 4*s**2 + 4*w**2
0
(3.2e-5*c**2*o**2*x**2 + 3.2e-5*c**2*v**2*x**2 + 6.4e-5*o*v*c**2*x**2)/
(8*a**2*e
**2)
0
0
Well that's interesting, it's got a 6th order term (there are 7
coeffs) but the constant and linear are 0 so you can factor out r**2
to get r**2*(r**4 + A*r**3 + B*r**2 + C); the linear term is also 0.
We'll capture the constants now for use later:
>>> one,AA,BB,z,CC,z,z = mypoly.as_monic().iter_all_coeffs()
So let's solve this (in an eyeblink, compared to the minutes it took
with the full form) and use cse on the results:
>>> solns = solve(x**4+A*x**3+B*x**2+C,x,simplified=False)
>>> for root in solns:
... rep,eq = cse(root)
... print rep
... print eq
... print
...
[(x0, A**2), (x1, 3*x0/256), (x2, B/16), (x3, x1 - x2), (x4, 32*x1),
(x5, B - x4
), (x6, x4/3), (x7, 8*x2), (x8, x6 - x7), (x9, x8**2), (x10, x5**2),
(x11, 3**(1
/2)), (x12, -64*x3*x6), (x13, 8*C), (x14, x12 + x13), (x15, -x13/4),
(x16, -x12/
4), (x17, x10 + x15 + x16), (x18, x14*x17), (x19, x10**(3/2)), (x20,
72*x6*x9),
(x21, x19 + x20), (x22, -x21*x5), (x23, x18 + x22), (x24, -8*x15),
(x25, -8*x16)
, (x26, x24 + x25), (x27, x23*x26), (x28, 3*x20), (x29, 8*x7), (x30,
12*x6), (x3
1, x29 - x30), (x32, x10*x31), (x33, x28 + x32), (x34, x28*x33/27),
(x35, x27 +
x34), (x36, x35**(1/2)), (x37, x11*x36/18), (x38, -x25/24), (x39, -
x24/24), (x40
, x10/54), (x41, x38 + x39 + x40), (x42, x29/2), (x43, x30/2), (x44,
x42 - x43),
(x45, x41*x44), (x46, x28/54), (x47, x37 + x45 + x46), (x48, x47**
(-1/3)), (x49
, x48**(-2)), (x50, 9*x49), (x51, -18*x39), (x52, -18*x38), (x53, x49**
(1/2)), (
x54, 3*x43), (x55, -3*x42), (x56, x54 + x55), (x57, x53*x56), (x58,
x10 + x50 +
x51 + x52 + x57), (x59, x48*x58), (x60, x59**(1/2))]
[-A/4 - x60/6 + (x48*(-x60*(x10 + x50 + x51 + x52 + x53*(-2*x54 -
2*x55)) + 54*A
*x53*x8))**(1/2)/(6*x60**(1/2))]
etc...
Since you have only 1 equation, you access it from eq as eq[0]. Those
reps are useful in their reversed order because you can feed them
directly to .subs()--again, in reversed order--and you will obtain
your original equation back.
So if you are going to be making some changes to parameters you could
do something like the following:
#do this once so it's always ready
>>> rep = list(reversed(rep))
>>> myrep=[(A,AA),(B,BB),(C,CC)] #AA etc...were captured previously
#define your constants
>>> constants={a:5.75,e:8e-6,c:0.25,x:20.,o:0.02,v:0.005,b:0.12,s:0.23,w:0.23,d:0.,f:20.,}
#update A, B, C
>>> abc = [(tmp[0], tmp[1].subs(constants)) for tmp in myrep]
#calculate your solution
>>> bigE = eq.subs(rep + abc)
>>> bigE
-0.58 - ((638203.292983406 + 4804.0368*(-18835904.2246641 +
3211610.35444874*I*3
**(1/2))**(1/3) + 9*(-18835904.2246641 + 3211610.35444874*I*3**(1/2))**
(2/3))/(-
18835904.2246641 + 3211610.35444874*I*3**(1/2))**(1/3))**(1/2)/6 + I*
((((638203.
292983406 + 4804.0368*(-18835904.2246641 + 3211610.35444874*I*3**(1/2))
**(1/3) +
9*(-18835904.2246641 + 3211610.35444874*I*3**(1/2))**(2/3))/
(-18835904.2246641
+ 3211610.35444874*I*3**(1/2))**(1/3))**(1/2)*(638203.292983406 -
9608.0736*(-18
835904.2246641 + 3211610.35444874*I*3**(1/2))**(1/3) + 9*
(-18835904.2246641 + 32
11610.35444874*I*3**(1/2))**(2/3)) - 50112.0*(-18835904.2246641 +
3211610.354448
74*I*3**(1/2))**(1/3))/(((638203.292983406 + 4804.0368*
(-18835904.2246641 + 3211
610.35444874*I*3**(1/2))**(1/3) + 9*(-18835904.2246641 +
3211610.35444874*I*3**(
1/2))**(2/3))/(-18835904.2246641 + 3211610.35444874*I*3**(1/2))**(1/3))
**(1/2)*(
-18835904.2246641 + 3211610.35444874*I*3**(1/2))**(1/3)))**(1/2)/6
>>> N(_)
-29.4350345287009 + 0.00248878631934675*I
And that was for the last of the 4 roots. {hmmm...that imaginary part
seems too big...I wonder if there's something lurking in there still.}
BUT...anyway, cse is your friend if you are doing repeated
calculation, and doing some human intelligent simplification (using
CAS to help you) is a "very good thing": the solution process goes
down from minutes (or hangs) to seconds.
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
-~----------~----~----~----~------~----~------~--~---
