Archive for June, 2009

Code cleanup and updating of examples

Tuesday, June 30th, 2009

A few days ago, I discovered the .get() method of the dictionary data type. It allows you to request the value corresponding to a key which a dictionary may not have, in which case you can specify what should be returned instead of rasing a KeyError. In PyDy, I had numerous cases where I was using if/else statements to check if a dictionary had a certain value, but by using the .get() method, I was able to eliminate about 50 or 60 lines of code and make the code a lot more readable and easier to maintain.

The progress on the trigonometric functions has been slow. Ondrej and I worked on some code together that handles most of the cases in the table by Fu et al., but not all. One of the things that has been a little challenging is to identify pi shifts in the argument to the trig functions and map the shift into the first quadrant (0, pi/2) and return the appropriate result, e.g.: sin(x+17*pi/18) –> sin(x + pi/18), etc. The approach has been to match x and r in the following generic argument:

x + r*pi

As I see it, there are three cases to consider: 1) x==0, r!=0, 2) r==0, x!=0, and 3) x!=0, and r!=0. I am splitting the logic up in this fashion and am trying to deal with each case in the appropriate manner. In most cases, r will be of the Rational class, and therefore the modulo operator is needed to bring it to the interval (0, 2*pi), so I implemented the necessary code and sent a patch to the Sympy mailing.

Although the trig functions aren’t quite as they need to be yet, I was still able to manually get the simplifications that are presented by Fu et al. to work using the builtin Sympy commands, which hopefully will mean that simplifying the implementing the algorithm will be feasible with the existing Sympy tools and not too many extra functions will need to be created.

Still working on trigonometric.py

Monday, June 22nd, 2009

I am still working on modifying trigonometric.py to make all of the trig methods behave as they would in Mathematica/Matlab/Maple. This involves getting the .eval method of each of them correct, and then from there I can work on implmenting a new trigsimp.

NSF Success!

Wednesday, June 17th, 2009

Today I found out that our NSF proposal to the control systems division was accepted! This means I will be funded for two years to study and work bicycles and control! This is a dream come true for me and I am very excited about it.

Soon, I hope to use PyDy to derive the model of the bicycle for our studies. Once trigsimp() is made better, it will be fully capable of this :)

sin, cos, tan, cot, sec, csc

Monday, June 15th, 2009

In my trigsimp branch of Sympy, I’ve modified all the .eval methods of sin(), cos(), tan(), and cot() so that they give the same behavior of Mathematica. Next on the list is to implement sec() and csc() in a similar fashion. Currently sec() and csc() do not exist in Sympy, so I will have to familiarize myself with all the other code that is in the trig functions that do exist. It seems to me that much of the code is redundant and may not be tested for, so I will attempt to do some housecleaning in that section while I am at it. The rule based trig simplification should then work fine, I have already started some basic tests of applying the rules that Fu et al. show in their examples, and Sympy handles them just fine.

Ondrej also helped me get trigsimp branch up and running on github:
git://github.com/hazelnusse/sympy.git

My trigsimp branch is there and I have created an issue for it here:
http://code.google.com/p/sympy/issues/detail?id=1475

Dyad and Inertia classes implemented

Sunday, June 7th, 2009

In order to be able to specify the inertia of a rigid body, I implemented the Dyad and Inertia classes. The Dyad class works with the UnitVector class and supports left and right dot multiplication by any other UnitVector or Vector expression (assuming the relative orientations have been specified). The Inertia class is subclassed from Dyad, and allows the user to specify the reference frame, along with the 6 inertia scalars that fully define the mass distribution of a rigid body.

An example of use for the Dyad class would be:

>>> from sympy import *
>>> from pydy import *
>>> t = Symbol(‘t’)
>>> q1 = Function(‘q1′)(t)
>>> A = ReferenceFrame(‘A’)
>>> B = A.rotate(‘B’, 1, q1)
>>> var(‘a b c’)
(a, b, c)
>>> a_Dyad = Dyad(a*A[1]*B[2] + b*A[3]*B[1] + c*A[2]*B[3])
>>> print a_Dyad
Dyad(a·a₁·b₂ + b·a₃·b₁ + c·a₂·b₃)

Most of the time in classical mechanics, dyads are used to keep track of the inertia of rigid bodies, so I also implemented an inertia class:

>>> I_dyad = Inertia(A, (1, 2, 3, 4, 5, 6))
print I_dyad
Dyad(a₁**2 + 2·a₂**2 + 3·a₃**2 + 4·a₁·a₂ + 4·a₂·a₁ + 5·a₂·a₃ + 5·a₃·a₂ + 6·a₁·a₃ + 6·a₃·a₁)

I need to polish the display of the dyad and inertia classes, as they currently don’t display quite right (a₁**2 should be a₁*a₁), but the mathematics is there and behaves as it should:

>>> print dot(A[1], a_Dyad)
a·b₂
>>> print dot(a_Dyad, A[1])
b·b₁

In addition to these classes and class methods, I also implemented a couple of convenience functions: InertiaForce() and InertiaTorque(), which simply form -m*a and -dot(alpha, I) – cross(omega, dot(I, omega)), respectively.

These new functions will be incorporated into the frame work that will allow for the derivation of the equations of motion for a system of particles and rigid bodies.

My finals will be done by Tuesday, I’m looking forward to some serious work on PyDy once exams are over :)

Progress on trigsimp(), .eval() method of sin, cos, tan, and PyDy printing

Tuesday, June 2nd, 2009

A few things I’ve been working on in the last week have been the .eval() methods for the sin, cos, and tan classes inside of sympy.functions.elementary.trigonometric. Each of these class methods control how something like

>>> sin(x + pi/2)

is evaluated. Ideally, this would nicely evaluate to cos(x), by default. Similar results are desirable for cos and tan. This is the behavior that Mathematica has. So far, I have implemented this behavior into the sin and cos functions, and should have tan done in the next few days. A few snags along the way have been dealing with the case when some of the arguments are Function or Derivative types, rather than Symbol, but I think I have figured this out. This will be a nice improvement of the default Sympy behavior and will make implementing trigsimp() easier.

Once the behavior of Table 1 of the paper by Fu et al. (see last post) has been implemented into these .eval() methods of sin, cos, tan, my next task will be to implement the rule based algorithm for trigsimp(). The first step will be to be able to compare the length of two trigonometric expressions, so I will be writing a function that will act as a metric for comparing trig expressions. Once this is done, one can compare a trig expression before and after a trig identity has been applied to it, and if the new expression is smaller, according to the metric, then we can accept it, otherwise we reject it. The algorithm is very customizable because you can try any sequence of trig identities to see if one approach works better than another, although hopefully the sequences of trig identities in the paper will prove to be sufficient for most problems.

Other developments on the PyDy front include adding a subs method to the Vector class, for the purpose of replacing time derivatives of the generalized coordinates with the expressions involving the generalized speeds. This was a simple matter of looping through the expressions corresponding to each of the UnitVector keys in the dictionary that is used to represent each Vector quantity.

The PyDy printer has been customized and now it by default will display sin(q1) as s1, cos(q1) as c1, etc. This really makes things more readable. I had a lot of issues getting unicode to work properly, but now the multiplication ‘*’ symbol is being displayed with a nice small unicode dot, and the display looks a lot cleaner:

(-r1·q5′ – r1·q3′·s4 – r2·q5′·c4)·b₁ + (r1·q4′ + r2·q4′·c4)·b₂ – r2·q4′·s4·b₃

The above output is the velocity of the mass center of the rolling torus example that I am using to test many of the functions I implement in PyDy. What isn’t visible above is that the b₁, b₂, and b₃ symbols are in bold and bright white, because they are unit vectors, so they stand out make the output easy to read. Another idea I have been toying with is to make the time derivative of the generalized coordinates a different color, and similarly with the time derivatives of the generalized speeds. Additionally I have customized the printer to display derivatives with an apostrophe after them, which is a lot easier to read and takes up a lot less space than D(q1(t), t), which was the default Sympy behavior. What would be really cool was if there was a Unicode dot syntax, like \dot{} in LaTeX, then things would be really small. Creating a custom LaTeX printer is on my list of things to do this summer, so that all of the equations that are need to include in a paper can be generated in a LaTeX form that is consistent with the syntax used in the dynamics literature.