Posts Tagged ‘Sympy’

Recursive methods for computing equations of motion

Tuesday, July 28th, 2009

In the last week, I added methods to the NewtonianReferenceFrame class of PyDy which start at the top of the reference frame tree and the point tree and do several things. The first method, recursive_subs(expr), is designed to take kinematic differential equations (which arise through a combination of linear nonholonomic constraints and/or holonomic constriants, and definitions of generalized speeds) and substitute them into all expressions for relative angular velocity between ReferenceFrames and relative velocity between Points. In doing so, the angular velocity and velocity expressions will involve only the independent time derivatives of coordinates. The second method, recursive_acc(), parses the Point/ReferenceFrame trees and forms the relative angular acceleration and acceleration. The third method, apply_gravitational_forces(), allows one to apply a gravitational forces to every particle and rigid body in the system with one command.

Currently, the NewtonianReferenceFrame class is acting like a container class for information about the system which is somewhat global in nature. For example, definitions of generalized coordinates, their time derivatives, and generalized speeds must all be assigned as attributes of this class in order for it to work correctly.

On the todo list is to write fr() and frstar(), which will also be recursive functions which will compute the equations of motion in one command. More work on the functions dealing with constraints is also needed, and there are at least two approaches that I could potentially implement for this. In particular, the approach taken with nonholonomic systems versus holonomic systems may necessitate different approaches. Currently, the work flow is as follows:
1) Form angular velocity and velocity expressions (linear in time derivatives of coordinates)
2) Form constraint equations arising from differentiated holonomic constraints and/or simple motion (rolling) constraints (m equations, linear in time derivatives of coordinates).
3) Define n-m generalized speeds as linear combinations of time derivatives of coordinates.
4) Solve generalized speed equations (item 3) for n-m ‘independent’ time derivatives of coordinates, in terms of the n-m generalized speeds and the m ‘dependent’ time derivatives of coordinates.
5) Solve constraint equations for m ‘dependent’ time derivatives of coordinates in terms of the ‘independent’ time derivatives of coordinates.
6) The complete set of kinematic differential equations will be the equations arising from item 5 and item 4.
7) Substitute these expressions into the expressions for angular velocity and linear velocity so that the constraints and definitions for generalized speeds are imposed, and so that the constrained partial velocities and partial angular velocities can be formed.

In taking this approach, the minimal set of n-m dynamic equations of motion will be formed, and they will obey all the constraints. These n-m equations, combined with the n kinematic differential equations comprise the complete set of motion equations that describe the system.

PyDy’s examples are being updated to reflect this new functionality, especially the rolling torus example. Additionally I have added a bicycle model which should soon have the nonlinear equations of motion for a very standard bicycle model.

Tags: , , , ,
Posted in PyDy, Sympy | No Comments »

GSoC developments

Monday, July 13th, 2009

In the last week, I finalized the locate method for the point class, as well as the .rel() and .vel() methods. Just as relative angular velocity between adjacent ReferenceFrame instances is stored, so is the relative position and relative velocity of points. By doing it in this fashion, the partial velocities can be formed incrementally, and there are some nice efficiencies to be realized when it comes time to form the equations of motion by avoiding the need to take redundant dot products. This will also translate into natural definitions for intermediate variables in the final form of the equations of motion.

One subtlety that is encountered in the kinematic chain of points is the possibility for closed kinematic chains. Typical examples are four and six bar linkages, or a bicycle with both wheels on the ground. These closed kinematic loops introduce holonomic constraint equations which are generally not solvable directly, but can be dealt with in differentiated form. PyDy now has a function called kinematic_chain() that deals with such loops and automatically forms the constraint equations and the differentiated holonomic constraint equations, and its use is illustrated in examples/fourbarlinkage.py in the git repository.

I have neared the completion of my updated sin,cos,tan,cot,sec,csc functions that handle arbitrary rational shifts by pi, and am in the process of testing and making things work with the hyperbolic functions that already exist in Sympy. This has been very challenging to get all the behavior that Mathematica does, but I think I have finally converged upon a solution that will be functional. Once this behavior is done, my next task will be to implement the Fu et al. paper so that trigsimp works. I will be putting PyDy development on hold until this is complete because PyDy needs to have some good trig simplification routines for many things to simplify as they should.

Tags: , , ,
Posted in PyDy, Sympy | No Comments »

Point class and doctest

Monday, July 6th, 2009

When I implemented the ReferenceFrame class with Ondrej, things were fairly straightfoward because the angular velocities of frames are 1) always described relative to another frame, and 2) because they add up in according to the angular velocity addition theorem. Points, on the other hand have proven a bit more complex to implement. When the user instantiates a new point using the locate method of the Point class, they provide the position Vector from the parent point to the new point. Optionally they may provide a ReferenceFrame which will determine how the velocity of the new point is determined (more on this later). Finding the position of one point relative to another is straightforward, and is implemented in the .rel() method of the Point class. What isn’t so straightforward is how to implement a velocity method to the point class.

When forming the velocity of a point, say point P1, one has to take the time derivative of a position vector of P1. This begs two questions:
1) Which relative position vector should be differentiated?
2) What frame should this position vector be differentiated in?
For forming the equations of motion, the answer to the first question is: the position from the inertial origin (or any other inertially fixed point) to the point in question. Similarly, the answer to the second question is: the inertial reference frame. This may not be what you want to do *ALL* of the time, but in every example I could think of, this is what is needed to form the equations of motion. The next question is how to do this efficiently.

The structure of the ReferenceFrame class is such that upon instantiation, ReferenceFrame objects have both an orientation and angular velocity that is relative to their parent frame. For the Point class, I have implemented the same thing — each new point has a position relative to it’s parent point. In both cases, this generates a tree structure with an inertial (Newtonian) frame at the top of the tree (in the case of the ReferenceFrame) and an inertial origin at the top of the tree (in the case of the Point class).

When the angular velocity of a frame is desired, one can use the .ang_vel() method of the ReferenceFrame class, and one needs to provide another ReferenceFrame whcih that angular velocity is computed relatve to (as in, the angular velocity of frame A relative to frame B). When the velocity of a point is desired on the other hand, what information should the user specify? My initial answer to this question is to make a method call .vel() with 2 default arguments that are None, and compute the velocity of the point with respect to the inertial origin unless otherwise specified. The issue with this is whether to pre-compute the time derivative (in the inertial frame) of each relative position vector, so that all one has to do is add the relative velocities of every point between the point in question and the inertial origin, so that the inertial velocity of that point is obtain. Additionally, when one uses the locate method, the option ReferenceFrame argument that can be specified will determine how this relative velocity is computed. Without specifying the optional argument direct time differentiation (in the inertial frame) of the position vector is performed, while specifying the optional argument results in the velocity being formed by omega (of the frame provided, relative to the inertial frame) crossed with the position vector. Ensuring that this subtlety is taken into account when implementing the .vel() method is a little tricky and that is what I’ve been working on for the last few days.

The Sympy documentation day last week was a great oppurtunity to write some documentation for how to subclass the StrPrinter, and it also exposed me to Sphinx and doctest. I started putting doctests in all of my docstrings, and I also started playing with Sphinx to figure out how to generate some nice documentation. I will be putting this up on docs.pydy.org in the next week.

I also spent a bunch of time learning how to manage a VPS. There is an excellent video blog series on www.guvnr.com where he goes through step by step (from zero to hero, as he puts it) how to setup a VPS to run the webserver nginx, a wordpress blog, etc. I’ve gotten the webserver up and running (now pydy.org is being redirected to the Google group page by my server instead of Ondrej’s), and I also got a nice secure ftp server up and running (using vsftpd in inetd mode with SSL). All in all very useful, but very time consuming!!! I’m exploring options for a wiki, probably either mediawiki or plone, but not sure yet. I managed to get plone up and running as well, but it looks like it might be more than is needed, plus I think it is very resource intensive.

Finally, over the last couple months I’ve been formally studying the book Programming in Python 3: A Complete Introduction to the Python Language by Mark Summerfield. I’ve learned about bunch of things I didn’t know, including named tuples, default dictionaries, and dictionary comprehensions (very handy for my Vector class), and I’m just about to dig into the modules chapter. Reading this book has really given me a better understanding of things like decorators (although they are still a little confusing), and also how to structure a module/package.

~Luke

Tags: , , ,
Posted in PyDy, Sympy | No Comments »

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.

Tags: ,
Posted in PyDy, Sympy | No Comments »

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.

Tags: , ,
Posted in Uncategorized | No Comments »