I don't know what all the problems are that systems run into, but as you see:
it happens. So I tend you lean on the CAS as an aid for the tedious part.
First, I get rid of the Reals and replace them with Rationals. I also throw
away the denominators which, in this case, do not contain any symbols. (I also
noticed the common term which I called z; I don't think the process is any
different if that is not identified.)
h[7] >>> [w.subs(dict([(r,nsimplify(r, rational=True)) for r in
w.atoms(Real)]))
.as_numer_denom()[0] for w in [z-(a0*(1 - x/2)*x),lam+2*y-z-0.005*x/2*x,
z-1*y-0
.743436700916726*y, x+y-conc]]
[2*z - a0*x*(2 - x), -400*z + 400*lam + 800*y - x**2, -174343670091673*y +
10000
0000000000*z, x + y - conc]
h[8] >>> neq = _
Pick the cherries: solve for the linear quantities:
h[9] >>> solve(neq[0], z)
[a0*x - a0*x**2/2]
h[10] >>> solve(neq[1].subs(z, _[0]), y)
[-lam/2 + a0*x/2 + x**2/800 - a0*x**2/4]
h[11] >>> neq[2].subs(z, h[9][0]).subs(y, _[0])
174343670091673*lam/2 + 25656329908327*a0*x/2 - 174343670091673*x**2/800 -
25656
329908327*a0*x**2/4
h[12] >>> solve(_, a0)
[(-69737468036669200*lam + 174343670091673*x**2)/(10262531963330800*x -
51312659
81665400*x**2)]
h[13] >>> neq[3].subs(y,h[10][0]).subs(a0, _[0])
x - conc - lam/2 + x*(-69737468036669200*lam +
174343670091673*x**2)/(2*(1026253
1963330800*x - 5131265981665400*x**2)) - x**2*(-69737468036669200*lam +
17434367
0091673*x**2)/(4*(10262531963330800*x - 5131265981665400*x**2)) + x**2/800
Let's separate it into numer and denom again:
h[14] >>> n,d=_.as_numer_denom()
h[14] >>> n
-1600*conc*(20525063926661600*x -
10262531963330800*x**2)*(41050127853323200*x -
20525063926661600*x**2) - 1600*x**2*(-69737468036669200*lam +
174343670091673*x
**2)*(20525063926661600*x - 10262531963330800*x**2) -
800*lam*(20525063926661600
*x - 10262531963330800*x**2)*(41050127853323200*x - 20525063926661600*x**2)
+ 2*
x**2*(20525063926661600*x - 10262531963330800*x**2)*(41050127853323200*x -
20525
063926661600*x**2) + 1600*x*(-69737468036669200*lam +
174343670091673*x**2)*(410
50127853323200*x - 20525063926661600*x**2) + 1600*x*(20525063926661600*x -
10262
531963330800*x**2)*(41050127853323200*x - 20525063926661600*x**2)
h[15] >>> n.expand()
-5254416365225369600000000000000000000*lam*x**2 -
134809039741934495380262692659
2000000*conc*x**2 + 1348090397419344953802626926592000000*x**3 +
134809039741934
4953802626926592000000*conc*x**3 +
5254416365225369600000000000000000000*lam*x**
3 - 1334954356506281529802626926592000000*x**4 -
1313604091306342400000000000000
000000*lam*x**4 - 337022599354836238450656731648000000*conc*x**4 +
3238865584417
72814450656731648000000*x**5 + 3284010228265856000000000000000000*x**6
It's not really x**6 difficult because it can be factored:
h[16] >>> factor(_)
13136040913063424000000*x**2*(2 - x)**2*(-100000000000000*lam -
25656329908327*c
onc + 25656329908327*x + 250000000000*x**2)
h[17] >>> solve(_, x)
[2, -25656329908327/500000000000 + (658247264364914528223938929 +
25656329908327
000000000000*conc + 100000000000000000000000000*lam)**(1/2)/500000000000,
0, -25
656329908327/500000000000 - (658247264364914528223938929 +
256563299083270000000
00000*conc + 100000000000000000000000000*lam)**(1/2)/500000000000]
Checking the denominator we see that two of those solutions will make the
denominator go to zero: the 2 and the 0.
h[18] >>> d
1600*(20525063926661600*x - 10262531963330800*x**2)*(41050127853323200*x -
20525
063926661600*x**2)
h[19] >>> factor(d)
337022599354836238450656731648000000*x**2*(2 - x)**2
We can substitute values into the others to test see what the solution is
numerically:
h[20] >>> [w.subs(lam,.1).subs(conc,.1) for w in h[17]]
[2.00000000000000, 0.487452050242467, 0, -103.112771683550]
If only a positive value for x makes sense then we have one solution after
tossing the 2, 0 and the negative one.
HTH,
/c
--
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.
