Numerical evaluation, how to speed up?

[Vinzenz]

Hi,

I am looking for a way to efficiently calculate a matrix (ca. 1-2 MB
pickled version) with quite large trig. expressions numerically, i.e.
for numerical values of it's symbols (ca. 20).

First I tried using subs, now I am using lambdify, to get a numpy
function. This function evaluates fast. However, lambdify takes a
really long time for this task and takes a lot of memory.

At the moment, I can only think of using lambdify element-wise. Also,
simplification would be an option, but I didn't find a good and fast
approach / combinations of functions to simplify the matrix.

Are there other options?

Kind Regards,
Vinzenz

[Chris Smith]

perhaps cse could help so you dn't have to keep recomputing
thing...can you pastebin the matrix so we can see it?

[Michael Hoffman]

Or at least an particularly nasty example of a matrix entry. pasting a 1 to 2 MB matrix is probably not advisable.

Also, if you can give the desired input and output that would be helpful.

Mike

On Tue, Oct 18, 2011 at 8:25 AM, Chris Smith <smi...@gmail.com> wrote:

perhaps cse could help so you dn't have to keep recomputing
thing...can you pastebin the matrix so we can see it?

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To post to this group, send email to sy...@googlegroups.com.
To unsubscribe from this group, send email to sympy+un...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/sympy?hl=en.

--
Michael Hoffman
ROSE Compiler Group
Lawrence Livermore National Laboratories
Department of Energy

[Aaron Meurer]

Apparently there's a bug in cse with matrices, though it's easy to
work around. It returns a list instead of a Matrix (the workaround is
just to pass it as Matrix(2, 2, list)).

In [183]: cse(Matrix([[sin(x), cos(x)], [-cos(x), sin(x)]]))
Out[183]: ([(x₀, cos(x)), (x₁, sin(x))], [x₁, x₀, -x₀, x₁])

But it's important that you cse the entire matrix at once. Otherwise,
you won't really be removing all the common subexpressions.

Aaron Meurer

[Michael Hoffman]

What's the complexity of cse

[Michael Hoffman]

Ugh, accidentally sent too fast.

What I meant to say was that if the cse is NP of any sort, and we're dealing with millions of entries I don't think cse is tenable.

Mike

[Chris Smith]

Not sure, but it uses the skip feature of the traversal so it doesn't
re-enter expressions that have been seen before.

[Chris Smith]

I think it would help a lot to see a section of what you've got and
what you want, even a 4x4

[Vinzenz]

On 18 Okt., 16:47, Michael Hoffman <archangel.associ...@gmail.com>
wrote:

Thanks for your answers.
Here is a small expression of one element of the matrix:
http://pastebin.com/AYbUu66t
How would you simplify it?

Here is a simplified (to about 1-2% length in few seconds) version of
the expression obtained with Mathematica:
http://pastebin.com/VMiGK0dJ

I already found an old thread, stating that the algorithms to combine
sin/cos arguments is not implemented yet.
Is this still true?

Running lambdify with parallel python on a 6-core Xeon I get the
result much much faster than on a i5. I don't know yet, if this is
just because of the machines or also a software matter. Does atlas or
similar have an effect on lamdify?

However, I need the expression for further processing in python as
well as in Matlab in a simplified version. I already tried the cse
method, but got strange results, i.e. some subexpressions were not
used at all and the result was not the same. Besides that, the
computation became slowlier (there were many small subexpressions). I
will investigate it again, maybe my fault somewhere. I am also
checking if I can simplify the expressions numerically, i.e. replace
summands which have little effect with constants.

Vinzenz

[Aaron Meurer]

Unfortunately, it still is. See
http://code.google.com/p/sympy/issues/detail?id=1676.

You can sometimes still do this by rewriting the expression in terms
of complex exponentials. Here, a is your expression:

In [94]: print expand(expand(a.rewrite(exp)).rewrite(cos))
-2*ddq1*sin(2*q2)*cos(q3)*cos(q4)/5 + 2*ddq1*sin(q4)*cos(2*q2)/5 -
2*ddq1*sin(q4)/5 + 2*ddq2*sin(q3)*cos(q2)*cos(q4)/5 -
2*ddq3*sin(q2)*cos(q3)*cos(q4)/5 + 2*ddq4*sin(q2)*sin(q3)*sin(q4)/5 -
4*dq1*dq2*sin(2*q2)*sin(q4)/5 - 4*dq1*dq2*cos(2*q2)*cos(q3)*cos(q4)/5
+ 2*dq1*dq3*sin(2*q2)*sin(q3)*cos(q4)/5 +
2*dq1*dq4*sin(2*q2)*sin(q4)*cos(q3)/5 + 2*dq1*dq4*cos(2*q2)*cos(q4)/5
- 2*dq1*dq4*cos(q4)/5 - 2*dq2**2*sin(q2)*sin(q3)*cos(q4)/5 +
2*dq3**2*sin(q2)*sin(q3)*cos(q4)/5 +
4*dq3*dq4*sin(q2)*sin(q4)*cos(q3)/5 +
2*dq4**2*sin(q2)*sin(q3)*cos(q4)/5

This might be slow, because expand() is slow, but it works. If it is
too slow, you might make it a little faster by turning off the hints
you don't need in expand.

You could then get a result similar to the Mathematica one by using
collect(). I don't know if there is a good general way to do
collection to reduce the number of times an expression appears. Maybe
someone else could comment on that.

By the way, it appears there is a bug in factor():

In [99]: factor(a.rewrite(exp))
Out[99]: nan

Aaron Meurer

>
> Running lambdify with parallel python on a 6-core Xeon I get the
> result much much faster than on a i5. I don't know yet, if this is
> just because of the machines or also a software matter. Does atlas or
> similar have an effect on lamdify?
>
> However, I need the expression for further processing in python as
> well as in Matlab in a simplified version. I already tried the cse
> method, but got strange results, i.e. some subexpressions were not
> used at all and the result was not the same. Besides that, the
> computation became slowlier (there were many small subexpressions). I
> will investigate it again, maybe my fault somewhere. I am also
> checking if I can simplify the expressions numerically, i.e. replace
> summands which have little effect with constants.

If the result was not the same, then there is a bug in cse.

>
> Vinzenz

[Aaron Meurer]

I created http://code.google.com/p/sympy/issues/detail?id=2818 for this.

Aaron Meurer

[Vinzenz]

On 2 Nov., 01:04, Aaron Meurer <asmeu...@gmail.com> wrote:

> Unfortunately, it still is.  Seehttp://code.google.com/p/sympy/issues/detail?id=1676.

I am impressed, thanks a lot.

Here is my current problem, still related to numerical evaluation, ...
somehow:
(see here for code: http://pastebin.com/QCRaKz5q )

What I have:
Equation of motion in form
tau = Y * pv = f(q, dq, ddq), e.g.

Y = Matrix([S('[4*ddq1/5 - 981*sin(q1)/100, 981*cos(q1)/100, ddq1,
4*ddq1*cos(q2)/5 + 117*ddq1/125 - 2*ddq2*cos(q2)/5 - 117*ddq2/125 -
4*dq1*dq2*sin(q2)/5 + 2*dq2**2*sin(q2)/5 - 981*sin(q1)*cos(q2)/100 +
981*sin(q2)*cos(q1)/100, -4*ddq1*sin(q2)/5 + 2*ddq2*sin(q2)/5 -
4*dq1*dq2*cos(q2)/5 + 2*dq2**2*cos(q2)/5 + 981*sin(q1)*sin(q2)/100 +
981*cos(q1)*cos(q2)/100, ddq1 - ddq2]'),
S('[ 0, 0,
0, -2*ddq1*cos(q2)/5 -
117*ddq1/125 + 117*ddq2/125 + 2*dq1**2*sin(q2)/5 + 981*sin(q1)*cos(q2)/
100 - 981*sin(q2)*cos(q1)/
100, 2*ddq1*sin(q2)/5 +
2*dq1**2*cos(q2)/5 - 981*sin(q1)*sin(q2)/100 - 981*cos(q1)*cos(q2)/
100, -ddq1 + ddq2]')
])

pv = Matrix(S('3.3113, 0.0513, -1.2524, 1.3106, 0.0065, -1.1819'))

tau = []
ddq = []
dq = []
for k in range(Y.rows):
tau += [Symbol('tau'+str(k+1))]
ddq += [Symbol('ddq'+str(k+1))]
dq += [Symbol('dq'+str(k+1))]

tau = Matrix(tau)


What I need:
Equation in different forms, such as

ddq = f( tau, q, dq )
with
ddq = [Symbol('ddq1'),Symbol('ddq2')]


First I tried to use sympy.solve, but it had problems with floats:
......
File "/usr/lib/python2.7/site-packages/sympy/polys/densearith.py",
line 1005, in dmp_pow
raise ValueError("can't raise polynomial to a negative power")
ValueError: can't raise polynomial to a negative power

I then used nsimplify to convert to rationals (also pv only).

indirectEq = (Y * pv) - tau

for k in range(indirectEq.rows):
indirectEq[k] = nsimplify(indirectEq[k], rational=True)
print indirectEq[k]

ddqExpr = solve(indirectEq, ddq, simplify=False, check=False)

Well, solve works then but I get something with really long longs,
e.g:
>>> print ddqExpr
{ddq1: -80000000000000000000000000*(-6500*sin(q2) + 1310600*cos(q2) +
1801827)*(13*sin(q2)/5000 - 6553*cos(q2)/12500 -
448216000000001/10000000000000000)*(327650000000000000000000000000*dq1**2*(-619888305664062500000000000000000000000000000000000000000000000000000000000000000000*sin(q2)
+....

As I need numerical evaluation I do:
# define variables for lambdify
my_vars = []
for link in range(Y.rows):
my_vars += ["q"+str(link+1),"dq"+str(link+1),"tau"+str(link+1)]

# create state equations
stateEq = zeros(2*Y.rows, 1)
for k in range(Y.rows):
print k
stateEq[k*2,0] = dq[k]
stateEq[k*2+1,0] = ddqExpr[ddq[k]]

# stateEq[k*2+1,0] =
(eqmot[ddq[k]].subs({sympy.conjugate(sympy.Symbol('q2')):sympy.Symbol('q2'),
#
sympy.conjugate(sympy.Symbol('q1')):sympy.Symbol('q1')}))
# stateEq[k*2+1,0] =
sympy.expand(sympy.expand((stateEq[k*2+1,0]).rewrite(sympy.exp)).rewrite(sympy.cos))



eqmotLambdified = lambdify(my_vars, stateEq)

import numpy as np

testvals = np.random.rand(3*Y.rows,)
eqmotLambdified(*testvals)

>>> eqmotLambdified(*testvals)
Traceback (most recent call last):
File "<console>", line 1, in <module>
File "<string>", line 1, in <lambda>
TypeError: unsupported operand type(s) for *: 'long' and
'numpy.float64'


Then I tried to use evalf(5) before lambdifying; the expression then
has huge floats and conjugates (am I missing some assumptions?!)
# create state equations
stateEq = zeros(2*Y.rows, 1)
for k in range(Y.rows):
print k
stateEq[k*2,0] = dq[k]
stateEq[k*2+1,0] = ddqExpr[ddq[k]].evalf(5)

>>> print stateEq[1].args[0]
-2.0454e+290*dq2**2*sin(conjugate(q2))*cos(q2)**3/
((2.704e-5*sin(q2)*sin(conjugate(q2)) - 0...............

I also tried to subs those conjugates(X) with x.
#.subs({sympy.conjugate(sympy.Symbol('q2')):sympy.Symbol('q2'),
sympy.conjugate(sympy.Symbol('q1')):sympy.Symbol('q1')})


Because of these problems and since solve already takes some time for
this small case, I use the
Second option:
It is possible to write:
tau = B(q) * ddq + rest(q,dq)
I determine B by coefs and invert it to get the expression for ddq.

This works, but stateEq conains very small values, which I like to
'chop'
0.0469945511680003*dq1*dq2*sin(q2)/(0.001363024*sin(2*q2) +
1.11022302462516e-16*cos(2*q1)*cos(2*q2) +
1.11022302462516e-16*cos(2*q1) - 2.77555756156289e-17*cos(4.......

In another post I found a chop function such as
def chop(expr, tol):
return sympy.Add(*[term for term in expr.as_Add if
term.as_coeff_terms()[0] >= tol])

but as_Add isn't available (tried .args instead). Any ideas how to get
rid of terms like 1e-17*x?

I really appreciate your comments.

Kind Regards,
Vinzenz

[Aaron Meurer]

Hi.

This is http://code.google.com/p/sympy/issues/detail?id=2814.

>
> I then used nsimplify to convert to rationals (also pv only).

In the git master, you can also just do solve(rational=True) to do
this automatically. Actually, apparently this is the default (though
I'm not sure if that's a good idea; Chris, why is this?).

>
> indirectEq = (Y * pv) - tau
>
> for k in range(indirectEq.rows):
>    indirectEq[k] = nsimplify(indirectEq[k], rational=True)
>    print indirectEq[k]
>
> ddqExpr = solve(indirectEq, ddq, simplify=False, check=False)
>
> Well, solve works then but I get something with really long longs,
> e.g:
>>>> print ddqExpr
> {ddq1: -80000000000000000000000000*(-6500*sin(q2) + 1310600*cos(q2) +
> 1801827)*(13*sin(q2)/5000 - 6553*cos(q2)/12500 -
> 448216000000001/10000000000000000)*(327650000000000000000000000000*dq1**2*(-619888305664062500000000000000000000000000000000000000000000000000000000000000000000*sin(q2)
> +....

This is because nsimplify chose the closest matching rationals, and
you had 0.0448216000000001, which it converted to
448216000000001/10000000000000000. You can chop off that last 1 by
evalf'ing your Matrix to the precision you want before calling
nsimplify (or solve(rational=True)), like

indirectEq = indirectEq.evalf(10)

If you want to assume your variables are real (or positive), then set
that when you initialize them, using Symbol('q2', positive=True). By
the way, you can also kill all conjugates by doing
expr.replace(conjugate, Id):

In [102]: (x + x.conjugate()).replace(conjugate, Id)
Out[102]: 2⋅x

This replaces all instances of the conjugate() function with the
identity function. It's a useful trick for getting rid of functions
you don't want when you forgot to set assumptions (I also use it
sometimes to kill heavisides when I forgot to set my symbols to be
positive).

You can just use expr.evalf(chop=True).

Aaron Meurer

[Vinzenz]

Hi,

> > This works, but stateEq conains very small values, which I like to
> > 'chop'
> > 0.0469945511680003*dq1*dq2*sin(q2)/(0.001363024*sin(2*q2) +
> > 1.11022302462516e-16*cos(2*q1)*cos(2*q2) +
> > 1.11022302462516e-16*cos(2*q1) - 2.77555756156289e-17*cos(4.......
>
> > In another post I found a chop function such as
> > def chop(expr, tol):
> >    return sympy.Add(*[term for term in expr.as_Add if
> > term.as_coeff_terms()[0] >= tol])
>
> You can just use expr.evalf(chop=True).
>
> Aaron Meurer
>

I tried evalf with chop a couple of times, but somehow it doesn't chop
in some cases.
E.g in this case I would expect a different result:

>>> myexpr = S('-1.954e-14*cos(2*q1)*cos(2*q2) + 5.3291e-15*cos(2*q1) - 14.365*cos(2*q2) + 274.5')
>>> myexpr
-1.954e-14⋅cos(2⋅q₁)⋅cos(2⋅q₂) + 5.3291e-15⋅cos(2⋅q₁) -
14.365⋅cos(2⋅q₂) + 274.5
>>> myexpr.evalf(5, chop=True)
-1.954e-14⋅cos(2⋅q₁)⋅cos(2⋅q₂) + 5.3291e-15⋅cos(2⋅q₁) -
14.365⋅cos(2⋅q₂) + 274.5
>>> myexpr.evalf(chop=True)
-1.954e-14⋅cos(2⋅q₁)⋅cos(2⋅q₂) + 5.3291e-15⋅cos(2⋅q₁) -
14.365⋅cos(2⋅q₂) + 274.5
>>>

How could I chop the small values in this case?

Regards,
Vinzenz

[Chris Smith]

Was there a problem with the chop function other than it didn't use
abs and was outdated?

>>> def chop(expr, tol):
... return Add(*[term for term in Add.make_args(expr) if
abs(term.as_coeff_Mul()[0]) >= tol])
...
>>> myexpr

-1.954e-14*cos(2*q1)*cos(2*q2) + 5.3291e-15*cos(2*q1) -
14.365*cos(2*q2) + 274.5

>>> chop(myexpr, 1e-12)
-14.365*cos(2*q2) + 274.5

/c

[Chris Smith]

This is long...hopefully you can follow and ask questions as theyarise.
First, the reason that Rationals come back is because you didn't
select checking. But it's a problem that the conjugates come back;
I'll track this down and fix it. I'll also fix it so that floats are restored
even if you don't check the answers.

That aside, here are your equations to solve:
```python>>> import textwrap>>> text = lambda x:
textwrap.TextWrapper(...     width=60,...
break_long_words=False...     ).fill(str(x))>>>>>> print
text(eqs)[-0.0052*ddq1*sin(q2) + 1.04848*ddq1*cos(q2) +1.4414616*ddq1
+ 0.0026*ddq2*sin(q2) - 0.52424*ddq2*cos(q2)- 0.0448216000000001*ddq2
- 1.04848*dq1*dq2*sin(q2) -0.0052*dq1*dq2*cos(q2) +
0.52424*dq2**2*sin(q2) +0.0026*dq2**2*cos(q2) - tau1 +
0.063765*sin(q1)*sin(q2) -12.856986*sin(q1)*cos(q2) -
32.483853*sin(q1) +12.856986*sin(q2)*cos(q1) +
0.063765*cos(q1)*cos(q2) +0.503253*cos(q1), 0.0026*ddq1*sin(q2) -
0.52424*ddq1*cos(q2)- 0.0448216000000001*ddq1 +
0.0448216000000001*ddq2 +0.52424*dq1**2*sin(q2) +
0.0026*dq1**2*cos(q2) - tau2 -0.063765*sin(q1)*sin(q2) +
12.856986*sin(q1)*cos(q2) -12.856986*sin(q2)*cos(q1) -
0.063765*cos(q1)*cos(q2)]```
Ans the variables that are contained therein:
```python>>> list(reduce(set.union, [eq.free_symbols for eq in
eqs]))[tau1, dq1, ddq1, tau2, dq2, q2, q1, ddq2]>>>
var(','.join([str(x) for x in _]))(tau1, dq1, ddq1, tau2, dq2, q2, q1,
ddq2)```
Use cse to simplify things:
```python>>> r,e=cse(eqs)>>> for i,ri in enumerate(r):...  print i,
ri...0 (x0, cos(q2))1 (x1, cos(q1))2 (x2, sin(q2))3 (x3, dq2**2)4 (x4,
sin(q1))5 (x5, 0.0448216000000001*ddq2)6 (x6, dq1**2)7 (x7, dq1*dq2)8
(x8, 0.0026*x0)9 (x9, 0.0026*x2)10 (x10, ddq1*x0)11 (x11,
0.52424*x2)12 (x12, 12.856986*x0*x4)13 (x13, 0.063765*x2*x4)14 (x14,
0.063765*x0*x1)15 (x15, 12.856986*x1*x2)16 (x16, x12 + x5)17 (x17, x13
+ x14 + x15)```
We want to keep ddq1 and ddq2,, so restore those values from x5 and x10:
```python>>> var('x5 x10')(x5, x10)>>> e=[ei.subs((r[5],r[10])) for ei in e]
>>> ans=solve(e,[ddq1,ddq2],check=False, simplify=False)
```

This will make coefficients into reals so we can look at the
real versions of the constants:
```python>>> def realcoeff(expr, deep=False, expo=False,
first=True):...     if first:...         c, p =
expr.as_content_primitive()...         return c.n()*realcoeff(p, deep,
expo, first=False)...     if expr.is_Rational:...         return
expr.n()...     elif not expr.args:...         return expr...     elif
expr.func is exp:...         return exp(realcoeff(expr.args[0], deep,
expo, first=False))...     elif expr.is_Pow:...         b, e =
expr.as_base_exp()...         if expo:...             e = e.n()...
    return realcoeff(b, deep, expo, first=False)**e...     elif deep
or expr.is_Mul:...         return expr.func(*[realcoeff(a, deep, expo,
first=False) for a in expr.args])...     elif expr.is_Add:...
return Add(*[Mul(float(c), realcoeff(m, deep, expo, first=False))...
             for c, m in [a.as_coeff_Mul() for a in
Add.make_args(expr)]])...     return expr...>>> print
text(realcoeff(ans[ddq1]))
1.25*(1000000.0*tau1 + 5200.0*x0*x7 - 503253.0*x1 -
1000000.0*x11*x3 + 1000000.0*x16 - 1000000.0*x17 +
1048480.0*x2*x7 - 1000000.0*x3*x8 +
32483853.0*x4)/(1310600.0*x0 - 6500.0*x2 + 1801827.0) -
1.25*(-8.0e+15*tau2*(1310600.0*x0 - 6500.0*x2 + 1801827.0) +
8.0e+15*x11*x6*(1310600.0*x0 - 6500.0*x2 + 1801827.0) +
8.0e+15*x16*(1310600.0*x0 - 6500.0*x2 + 1801827.0) -
8.0e+15*x17*(1310600.0*x0 - 6500.0*x2 + 1801827.0) +
8.0e+15*x6*x8*(1310600.0*x0 - 6500.0*x2 + 1801827.0) +
1.0*(-5.2424e+15*x0 + 1.0e+16*x9 -
448216000000001.0)*(1000000.0*tau1 + 5200.0*x0*x7 -
503253.0*x1 - 1000000.0*x11*x3 + 1000000.0*x16 -
1000000.0*x17 + 1048480.0*x2*x7 - 1000000.0*x3*x8 +
32483853.0*x4))/((-5.2424e+15*x0 + 1.0e+16*x9 -
448216000000001.0)*(1310600.0*x0 - 6500.0*x2 + 1801827.0))
>>> print text(realcoeff(ans[ddq2]))
0.0125*(-8.0e+15*tau2*(1310600.0*x0 - 6500.0*x2 + 1801827.0)
+ 8.0e+15*x11*x6*(1310600.0*x0 - 6500.0*x2 + 1801827.0) +
8.0e+15*x16*(1310600.0*x0 - 6500.0*x2 + 1801827.0) -
8.0e+15*x17*(1310600.0*x0 - 6500.0*x2 + 1801827.0) +
8.0e+15*x6*x8*(1310600.0*x0 - 6500.0*x2 + 1801827.0) +
1.0*(-5.2424e+15*x0 + 1.0e+16*x9 -
448216000000001.0)*(1000000.0*tau1 + 5200.0*x0*x7 -
503253.0*x1 - 1000000.0*x11*x3 + 1000000.0*x16 -
1000000.0*x17 + 1048480.0*x2*x7 - 1000000.0*x3*x8 +
32483853.0*x4))/((-6553.0*x0 + 12500.0*x9)*(-5.2424e+15*x0 +
1.0e+16*x9 - 448216000000001.0))
```

Now you could do substitutions based on the reps that were captured
with cseafter evaluating them at the desired numerical values.
Another alternative would be to use the condense function in pull
requesthttps://github.com/sympy/sympy/pull/598 :
```python>>> e1, r1=condense(eqs[0],ddq1)>>> solve(e1,
ddq1)[-C0/C1]>>> a1=_>>> e2,
r2=condense(eqs[1].subs(ddq1,a1[0]).subs(r1),ddq2)>>> solve(e2,
ddq2)[(C0 + C1)/(-C2 - C3)]>>> a2=_```
Now I'm going to substitute the constants back in to show that ddq1
and ddq2 havebeen cleared, but you would probably want to numerically
evaluate them or evenuse cse to identify the basic things that need
only be computed once.
```python>>> ddq1ans = a1[0].subs(r1).subs(ddq2, a2[0].subs(r2))>>>
ddq2ans = a2[0].subs(r2)>>> ddq2ans.has(ddq1,ddq2)False>>>
ddq1ans.has(ddq1,ddq2)False```

[Chris Smith]

Well that was horrid. It's a better at http://pastebin.com/ZnAqpbHD.
Do we have any markup control here?

[Chris Smith]

OK, I've modified solve to handle Floats better: if rational=False,
nothing happens to Floats (None, the default, allows the Floats to be
made Rational but they are converted back at the end; True converts
them and leaves them as Rationals in the answer):

Most of this is review but just so you can see that the time is much quicker:

>>> from time import time

>>> t=time();r,e = cse(eqs);eqsnew=[ei.subs([r[5],r[10]]) for ei in e];ans = solve(eqsnew, ddq1, ddq2, rational=False, check=0, simplify=0, manual=1);r2, e2 = cse(flatten(ans));print time()-t

0.725000143051

The manual option is exercised since polys doesn't like the floats but
the system is trivially solvable since it is linear in ddq1 and ddq2.

[Aaron Meurer]

If you use gist.github.com and choose Markdown as the language, then
you can use full markdown markup, which is exactly the same as in
GitHub pull request comments. See for example
https://gist.github.com/1358969.

Aaron Meurer

[Chris Smith]

On Sat, Nov 12, 2011 at 1:16 AM, Aaron Meurer <asme...@gmail.com> wrote:
> If you use gist.github.com and choose Markdown as the language, then
> you can use full markdown markup, which is exactly the same as in
> GitHub pull request comments.  See for example
> https://gist.github.com/1358969.

OK, thanks.

Here's about as compact as I can present the numerical calculation
with the new solve pull request mods:

variables and equations

v=var('tau1, dq1, ddq1, tau2, dq2, q2, q1, ddq2')
eqs = [-0.0052*ddq1*sin(q2) + 1.04848*ddq1*cos(q2) +1.4414616*ddq1

+ 0.0026*ddq2*sin(q2) - 0.52424*ddq2*cos(q2)- 0.0448216000000001*ddq2
- 1.04848*dq1*dq2*sin(q2) -0.0052*dq1*dq2*cos(q2) +
0.52424*dq2**2*sin(q2) +0.0026*dq2**2*cos(q2) - tau1 +
0.063765*sin(q1)*sin(q2) -12.856986*sin(q1)*cos(q2) -
32.483853*sin(q1) +12.856986*sin(q2)*cos(q1) +
0.063765*cos(q1)*cos(q2) +0.503253*cos(q1), 0.0026*ddq1*sin(q2) -
0.52424*ddq1*cos(q2)- 0.0448216000000001*ddq1 +
0.0448216000000001*ddq2 +0.52424*dq1**2*sin(q2) +
0.0026*dq1**2*cos(q2) - tau2 -0.063765*sin(q1)*sin(q2) +
12.856986*sin(q1)*cos(q2) -12.856986*sin(q2)*cos(q1) -
0.063765*cos(q1)*cos(q2)]

Put equations in cse form

r,e=cse(eqs)
r.reverse()

Restore the symbols of interest by putting those variables back

restore=[ddq1,ddq2]
rr=[ri for ri in r if ri[1].free_symbols & set(restore)]
e = [ei.subs(rr) for ei in e]

Solve this cse system

ans=solve(e,restore,check=0,simplify=0,rational=False,manual=1)

cse the results

r2,e2=cse(ans[0]) # there's only one solution
r2.reverse()

make the replacements

import random
reps=zip(v,[random.random() for vi in v])

replacements of equation cse

rn=[(ri[0],ri[1].subs(reps)) for ri in r]

replacements for answer cse

r2n=[(ri[0],ri[1].subs(rn)) for ri in r2]

final numerical answer

en=[ei.subs(r2n).subs(reps) for ei in e2]

>>> en
[2.09579142350940, 30.2909076894234]

>>> for i in reps:
... print i
...
(tau1, 0.18622320944903792)
(dq1, 0.77302534514532184)
(ddq1, 0.63839976360929995)
(tau2, 0.55417626041067558)
(dq2, 0.39107653059970293)
(q2, 0.15440282862677301)
(q1, 0.027723538970289119)
(ddq2, 0.38258056784276506)

[Chris Smith]

>>>> for i in reps:
> ...  print i
> ...
> (tau1, 0.18622320944903792)
> (dq1, 0.77302534514532184)
> (ddq1, 0.63839976360929995)
> (tau2, 0.55417626041067558)
> (dq2, 0.39107653059970293)
> (q2, 0.15440282862677301)
> (q1, 0.027723538970289119)
> (ddq2, 0.38258056784276506)
>

Of course the ddq1 and ddq2 values here were not used...I didn't
remove them from v when assigning random variables. But when checking
this numerically, I wonder if there is some mistake in something of
the preceding since a straight solution after substituting in values
doesn't agree with the above:

>>> v=var('tau1, dq1, ddq1, tau2, dq2, q2, q1, ddq2')
>>> eqs = [-0.0052*ddq1*sin(q2) + 1.04848*ddq1*cos(q2) +1.4414616*ddq1

... + 0.0026*ddq2*sin(q2) - 0.52424*ddq2*cos(q2)- 0.0448216000000001*ddq2
... - 1.04848*dq1*dq2*sin(q2) -0.0052*dq1*dq2*cos(q2) +
... 0.52424*dq2**2*sin(q2) +0.0026*dq2**2*cos(q2) - tau1 +
... 0.063765*sin(q1)*sin(q2) -12.856986*sin(q1)*cos(q2) -
... 32.483853*sin(q1) +12.856986*sin(q2)*cos(q1) +
... 0.063765*cos(q1)*cos(q2) +0.503253*cos(q1), 0.0026*ddq1*sin(q2) -
... 0.52424*ddq1*cos(q2)- 0.0448216000000001*ddq1 +
... 0.0448216000000001*ddq2 +0.52424*dq1**2*sin(q2) +
... 0.0026*dq1**2*cos(q2) - tau2 -0.063765*sin(q1)*sin(q2) +
... 12.856986*sin(q1)*cos(q2) -12.856986*sin(q2)*cos(q1) -
... 0.063765*cos(q1)*cos(q2)]
>>> reps=[
... (tau1, 0.18622320944903792),
... (dq1, 0.77302534514532184),
... (tau2, 0.55417626041067558),
... (dq2, 0.39107653059970293),
... (q2, 0.15440282862677301),
... (q1, 0.027723538970289119),
... ]
>>> neq=[ei.subs(reps) for ei in eqs]
>>> nsolve(neq,[ddq1,ddq2],[1,1])
matrix(
[['-5.7720493053172'],
['-23.5214066555216']])
>>> neq
[2.47666869757445*ddq1 - 0.562425148787224*ddq2 + 1.06642319741892, -0.562425148
787224*ddq1 + 0.0448216000000001*ddq2 - 2.19207860879909]
>>> solve(neq,[ddq1,ddq2])
{ddq1: -5.77204930531720, ddq2: -23.5214066555216}

hmmm...