Flux

I changed the code to the rolling torus example pretty significantly. Little fixes in how the arrows representing the body fixed unit vectors are now fixed. Additionally, I added some code that calculates the kinetic, potential and total energy of the system, and plots it. The scipy solver seems to do a pretty good (not perfect) job of maintaining constant energy. I think tightening up the error tolerances would improve this, but the simulation is well behaved at this point so I didn’t bother.

I just added an example of use for the derivation of the equations of a gyrostat using PyDy.  It is located in the examples/gyrostat/ subdirectory.

In case you aren’t familiar, here is the definition of a gyrostat:

A gyrostat is a mechanical system which is comprised of more than one body and yet has the rigid body property that its inertia components are time independent constants.

The reason I became interested in gyrostats is because in modeling the bicycle, the rear frame and the rear wheel can be treated as a gyrostat, and in doing so, the number of parameters that appear in the equations of motion will be reduced by two.  The same can be done for the front fork and handlebar assembly and front wheel.

In addition to the Python script that generates the equations, I took the time to do a complete write up of the model in LaTeX and generate a nice pdf of it.

I haven’t blogged in about two weeks, and there has been a lot of progress in PyDy. For starters, PyDy now ‘automatically’ derives the equations of motion for a rolling disc and a rigid body in space. I haven’t blogged in about two weeks, and there has been a lot of progress in PyDy. For starters, PyDy now ‘automatically’ derives the equations of motion for a rolling disc and a rigid body in space. One of the more useful methods I implemented is the output_eoms() method in the NewtonianReferenceFrame class. Once the equations of motion have been derived, this method creates python functions for them (in the form that scipy.odeint wants them in) and writes them to a file that you specify. Once you have this, you can import the function and use it for numerical integration. Additionally, this method also outputs a function that calculates the generalized speeds given the coordinates and their time derivatives, which is useful if you know the initial conditions of the coordinates time derivatives but need to determine the initial conditions of the generalized speeds when you want to perform a numerical integration. Finally, this method also optionally outputs an animate function, which takes a variable number of Points and (Axis, Angle) tuples and outputs a function which given the coordinates and parameters, will compute the inertial position and orientation of Points and/or Reference Frames. This is useful when you want to animate your simulation.
Both the rolling_disc.py and rigidbody.py files in the pydy/examples folder demonstrate the usage of these new methods.

One issue that has been repeatedly rearing its ugly head is the use of Function instances instead of Symbol instances. For everything except time differentiation, PyDy expressions could work fine with Symbols instead of Functions. I have become to realize that everything just works faster and more efficiently in Sympy with Symbols, so I will be working in the next week or so to replace the use of Functions in PyDy with Symbol instances instead. This should allow lots of things to work much better without having to call a bunch of helper methods to get expressions in their simplest form, or having to repeatedly substitute and back substitute with Symbols and Functions. For time differentiation, I will keep track of a list of symbols that represent the generalized coordinates, and then any expression that needs time differentiation will be examined to see if it has any occurrences of those coordinates. For expressions involving those coordinates, the partial derivative will be taken with respect to each coordinate, and then multiplied by the symbol that will be used to represent the time derivative of that coordinate, in essence a manual (as opposed to letting Sympy take care of it by using sympy Function) application of the chain rule. The user will be required to eclare these coordinate symbols in a special way, specifically using the .setcoords() method. Finally, I will implement method that will replace the symbols with Function instances (depending on time), so that if one wants to play with an expression interactively, all that they will have to do is call this function so that things like q1 = Symbol(’q1′) get replaced with q1 = Function(’q1′)(t), and then time differentiation will work just fine.

Last week I spent most of my time writing and testing new trigonometric functions for Sympy. Sympy already had sin, cos, tan, cot, asin, acos, atan, acot, but did not have sec, csc, asec, or acsc. So I added those four new functions and am in the process of writing very extensive tests for them.

The reason for writing these new functions is so that the trig functions have the same behavior as the ones in Mathematica. There are a few things that Mathematica does that is very nice, one of them being that any arguments that have a pi shift (as in sin(a + b*pi)), will always return another trig function whose pi shift is within the interval between 0 and pi/4. For example,

Getting all the behavior to work with the inverse trig functions is also implemented, as well as the special values of the trig functions for arguments which are integer multiples of pi/10. The next task is going to be to get all of this code working with the rest of Sympy so it can be put into place.

This behavior is essential to making the trigonometric simplification routine in Sympy to be more capable, and will greatly benefit both Sympy and PyDy.

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

Sympy 0.6.5 currently doesn’t allow for you to solve an expression, or differentiate with respect to, anything but a sympy Symbol.  This is problematic for PyDy because generally the generalized coordinates that are used are not just sympy Symbols, but instead are sympy Functions, implicitly dependent upon time.  So what I did was to implement the code that would allow solve to solve for not just Symbols, but also for Functions and Derivatives. So as a trivial example, the following code now works:

The patch that allows this behavior should be pushed into the main branch of sympy pretty soon, so if you clone the repo in the next few days, you should have this functionality as well.

Like I said, the example is trivial, but now, this will allow PyDy to solve the dynamic equations of motion (linear in the second time derivatives of the coordinates and/or first time derivatives of the generalized speeds). So as long as the linear solvers are up to the task of solving long, long, expressions, PyDy should be able to get the equations of motion in a form that will be directly ready for numerical integration, a key step to the task of generating animations. This is a good step forward.

The next task I’ll be tackling is to fix diff() so that it can differentiate with respect to objects other than the sympy Symbol, i.e., so that you can differentiate with respect to Functions and Derivative instances. This should be a very similar patch and should be done very soon.

Once that is done, I plan on spending some time on improving trigsimp(). Having this function work well will really help PyDy because with these long chains of rotations, the number of trigonometric terms becomes lengthy and having the most compact representation is ideal.