Sympy.physics.mechanics and its approach

158 views
Skip to first unread message

Tim Lahey

unread,
Apr 6, 2012, 5:12:13 PM4/6/12
to sy...@googlegroups.com
Hi,

I now know why I was confused about SymPy supporting functions that
depend upon x,y,z,t. It's because the sympy.physics.mechanics doesn't
support it (since it requires this dynamicsymbol class). It assumes
that things are functions of time only.

Looking through the code for sympy.physics.mechanics, I've noticed
that it's very slanted towards the Kane and Levinson approach, even in
all the supporting functions and classes and as a result, it forces
one into something that they might not be comfortable with.

One example I'm thinking of is the approach of how it handles inertia.
If I have an inertia matrix in a defined frame, how am I supposed to
include it? If I know Ixx, Iyy, Izz and the cross-products of inertia
how do I use them? The example in the rigidbody.py file is,

>>> from sympy import Symbol
>>> from sympy.physics.mechanics import ReferenceFrame, Point, RigidBody
>>> from sympy.physics.mechanics import outer
>>> m = Symbol('m')
>>> A = ReferenceFrame('A')
>>> P = Point('P')
>>> I = outer (A.x, A.x)
>>> Inertia_tuple = (I, P)
>>> B = RigidBody('B', P, A, m, Inertia_tuple)
>>> # Or you could change them afterwards
>>> m2 = Symbol('m2')
>>> B.mass = m2

I don't see how to use the known info about I in this example. Looking
through functions.py, there's an inertia function that does what I
want, but this isn't clear if one is trying to define a rigid body
because they'd look at the rigid body class. Plus, printing out the
inertia matrix would be a nice option, rather than the given form.

How is parallel axis theorem carried out in this code? I can't find an
example. A parallel axis theorem matrix is defined in
sympy.physics.matrices, but I doubt it's used. What about rotation of
the inertia matrix?

Next, is the kinematic equations function. First, they seem to be hard
coded based upon the sequence of rotations. Second, what is meant by u
and qdot? Nowhere is this defined or explained, even by just giving a
reference. There's rotation matrices about the 1-, 2-, and 3-axes
available in sympy.matrices (if needed, but I'm guessing there's
another way you'd prefer).

Trying to puzzle out the examples in the documentation is difficult
since there's lots of Python code, but really no explanation as to
what the symbols are. What are u1, u2, u3, q1, q2, q3 and how do they
relate to the figure? It's mentioned that they're configuration and
speed variables but it's not stated which are which (it's stated only
at the end of the example). Plus the choice of names (configuration
and speed) is definitely a Kane thing. Most texts I've read would have
called them position and velocity variables.

What does orientnew do? It's not clear from the example. I'm assuming
it's a rotation about an axis (or presumably an angle). Plus, the
frames Y and L should be shown on the diagram (plus R should be moved
next to the specific frame), especially since you're defining angular
velocities in one of them.

What is Dmc (put this on the figure)? The r * L.z confused me a bit. I
figured out that you're saying a distance r in the L.z direction, but
I rarely write out things that way. I'm normally specifying the frame
and the complete vector. I get that's what you're doing, but could I
pass a matrix or a tuple and multiply by L? I'd probably create a
vector in the L frame and use that, but I get that r * L.z is a
shortcut. Nowhere is just defining a vector in a frame done in the
example on its own.

What does v2pt_theory and a2pt_theory do? (Also, why the theory
suffix?). I'm guessing that Dmc is the centre of mass of the disc, but
nowhere is that stated. I'm also guessing that v2pt_theory and
a2pt_theory are the velocity and acceleration of Dmc with respect to
point C.

Why compute the kinematic differential equations by hand? Why not show
how to do it with the module? It seems like a lost opportunity.

I'm coming at this from the perspective of one who has been involved
in the teaching of a number of physics courses. If you define a symbol
in your code, it should show up in the text and if at all possible, in
a corresponding figure. Also, it would be nice if the output was LaTeX
output to show the full utility of the package.

Thoughts?

Tim.

--
Tim Lahey
PhD Candidate, Systems Design Engineering
University of Waterloo
http://about.me/tjlahey

Jason Moore

unread,
Apr 6, 2012, 8:15:09 PM4/6/12
to sy...@googlegroups.com
Tim,

Responses intertwined below:


On Friday, April 6, 2012 2:12:13 PM UTC-7, Tim Lahey wrote:
Hi,

I now know why I was confused about SymPy supporting functions that
depend upon x,y,z,t. It's because the sympy.physics.mechanics doesn't
support it (since it requires this dynamicsymbol class). It assumes
that things are functions of time only.

Looking through the code for sympy.physics.mechanics, I've noticed
that it's very slanted towards the Kane and Levinson approach, even in
all the supporting functions and classes and as a result, it forces
one into something that they might not be comfortable with.

I would beg differ. The equation of motion derivation is completely separate from the kinematic and kinetic definitions. As is stands, any method of equation of motion derivation can be used to solve classical multi-body dynamics with the currently defined classes. We haven't written any examples for methods other than Kane's yet due to sparse resources to do so. We also haven't implemented automated classes for any other methods yet, for the same reason. One GoSC code applicant hopes to remedy both of those issues this summer.
 

One example I'm thinking of is the approach of how it handles inertia.
If I have an inertia matrix in a defined frame, how am I supposed to
include it? If I know Ixx, Iyy, Izz and the cross-products of inertia
how do I use them? The example in the rigidbody.py file is,

    >>> from sympy import Symbol
    >>> from sympy.physics.mechanics import ReferenceFrame, Point, RigidBody
    >>> from sympy.physics.mechanics import outer
    >>> m = Symbol('m')
    >>> A = ReferenceFrame('A')
    >>> P = Point('P')
    >>> I = outer (A.x, A.x)
    >>> Inertia_tuple = (I, P)
    >>> B = RigidBody('B', P, A, m, Inertia_tuple)
    >>> # Or you could change them afterwards
    >>> m2 = Symbol('m2')
    >>> B.mass = m2

I don't see how to use the known info about I in this example. Looking
through functions.py, there's an inertia function that does what I
want, but this isn't clear if one is trying to define a rigid body
because they'd look at the rigid body class. Plus, printing out the
inertia matrix would be a nice option, rather than the given form.

We have a dyadic class to handle inertia. Dyadics are typically not as popular as tensor formulations but they have the advantage of being basis independent (unlike tensors). We probably need more functionality to expose the dyadics as tensors so that they folks who don't want to learn about dyadics may have an easier time working with them. The current inertia function is one such function that does that. It allows you to define an inertia dyadic with a tensor like formulation.
 

How is parallel axis theorem carried out in this code? I can't find an
example. A parallel axis theorem matrix is defined in
sympy.physics.matrices, but I doubt it's used. What about rotation of
the inertia matrix?

We don't have an explicit parallel axis method or function. Here is an example of how you can write the parallel axis function with dyadic notation: https://github.com/angadhn/sympy/blob/a368bc0c871a75cbea0ffe1a2efc553fced1cbd6/sympy/physics/mechanics/kane.py#L471

We'd love for someone to formalize this into some functions that follow the more traditional ideas for ease of understanding.
 

Next, is the kinematic equations function. First, they seem to be hard
coded based upon the sequence of rotations. Second, what is meant by u
and qdot? Nowhere is this defined or explained, even by just giving a
reference. There's rotation matrices about the 1-, 2-, and 3-axes
available in sympy.matrices (if needed, but I'm guessing there's
another way you'd prefer).

As with all open source projects, the documentation is sometimes weak. Please submit issues with regards to missing documentation and we all can try to remedy it.
 

Trying to puzzle out the examples in the documentation is difficult
since there's lots of Python code, but really no explanation as to
what the symbols are. What are u1, u2, u3, q1, q2, q3 and how do they
relate to the figure? It's mentioned that they're configuration and
speed variables but it's not stated which are which (it's stated only
at the end of the example). Plus the choice of names (configuration
and speed) is definitely a Kane thing. Most texts I've read would have
called them position and velocity variables.


We'd love to have some examples written for different formulations and for the current documentation to be improved. We welcome any pull requests or issues that point these out specifically.
 

What does orientnew do? It's not clear from the example. I'm assuming
it's a rotation about an axis (or presumably an angle). Plus, the
frames Y and L should be shown on the diagram (plus R should be moved
next to the specific frame), especially since you're defining angular
velocities in one of them.

Same as above.
 

What is Dmc (put this on the figure)? The r * L.z confused me a bit. I
figured out that you're saying a distance r in the L.z direction, but
I rarely write out things that way. I'm normally specifying the frame
and the complete vector. I get that's what you're doing, but could I
pass a matrix or a tuple and multiply by L? I'd probably create a
vector in the L frame and use that, but I get that r * L.z is a
shortcut. Nowhere is just defining a vector in a frame done in the
example on its own.

We have a Vector class that, by design, was not intended to be generated explicitly. The approach is to create a reference frame and specify vectors using the basis of the reference frame.

It would be relatively easy to add vector creation functions/methods as you suggest that interact with the Vector class, giving a different view on the creation and manipulation of vectors.

 

What does v2pt_theory and a2pt_theory do? (Also, why the theory
suffix?). I'm guessing that Dmc is the centre of mass of the disc, but
nowhere is that stated. I'm also guessing that v2pt_theory and
a2pt_theory are the velocity and acceleration of Dmc with respect to
point C.

v2pt and a2pt are helper methods to define velocities of points on the same rigid body. You don't have to use them if you don't want to.
 

Why compute the kinematic differential equations by hand? Why not show
how to do it with the module? It seems like a lost opportunity.

I'm coming at this from the perspective of one who has been involved
in the teaching of a number of physics courses. If you define a symbol
in your code, it should show up in the text and if at all possible, in
a corresponding figure. Also, it would be nice if the output was LaTeX
output to show the full utility of the package.

Thoughts?


We'd love to see some pull requests implementing features that you've suggested. We'd also gladly fix and improve and documentation errors or clarity issues that you find. The easy way to keep track of that would be submitting them as issues.

Thanks for your review!
 

Tim.

--
Tim Lahey
PhD Candidate, Systems Design Engineering
University of Waterloo
http://about.me/tjlahey


Tim Lahey

unread,
Apr 7, 2012, 12:56:15 AM4/7/12
to sy...@googlegroups.com
Jason,

Responses interleaved below,

On Fri, Apr 6, 2012 at 8:15 PM, Jason Moore <moore...@gmail.com> wrote:
>
> I would beg differ. The equation of motion derivation is completely separate
> from the kinematic and kinetic definitions. As is stands, any method of
> equation of motion derivation can be used to solve classical multi-body
> dynamics with the currently defined classes. We haven't written any examples
> for methods other than Kane's yet due to sparse resources to do so. We also
> haven't implemented automated classes for any other methods yet, for the
> same reason. One GoSC code applicant hopes to remedy both of those issues
> this summer.


True, you can use any approach for the kinematic definitions and still
derive the equations of motion using a different approach. That said,
what I'm saying is that the mechanics package's approach for the
kinematics is very different than what one would use if it was based
upon a different text. The approach to a vector illustrates this. The
approach I learnt for advanced dynamics embeds the frame as part of
the vector. Both approaches aren't wrong, just different.

The problem I have is that a lot of decisions have been fixed so
there's little flexibility. I'll probably have to ignore the whole
thing and write my own. For example, the decision representing the
inertia as a dyad rather than a matrix or a tensor. By choosing a
dyad, I can't just pull the matrix out and use matrix algebra tools
with it (e.g., matrix multiplication for rotations). The same problem
occurs with the use of "dynamicsymbol" which assumes that the symbol
is only a function of time. This precludes flexible bodies. This alone
means I can't use it.

>
> We have a dyadic class to handle inertia. Dyadics are typically not as
> popular as tensor formulations but they have the advantage of being basis
> independent (unlike tensors). We probably need more functionality to expose
> the dyadics as tensors so that they folks who don't want to learn about
> dyadics may have an easier time working with them. The current inertia
> function is one such function that does that. It allows you to define an
> inertia dyadic with a tensor like formulation.

I understand about inertia dyads, but they have limitations of their own.

>
> We don't have an explicit parallel axis method or function. Here is an
> example of how you can write the parallel axis function with dyadic
> notation:
> https://github.com/angadhn/sympy/blob/a368bc0c871a75cbea0ffe1a2efc553fced1cbd6/sympy/physics/mechanics/kane.py#L471
>
> We'd love for someone to formalize this into some functions that follow the
> more traditional ideas for ease of understanding.

This is one of the problems of using an inertia dyad. Most texts will
present an inertia matrix (or tensor) and give the parallel axis
theorem as a matrix. Plus, rotation of the inertia matrix can be
handled easily using rotation matrices. Given good matrix support,
this can be more efficient.

>
> As with all open source projects, the documentation is sometimes weak.
> Please submit issues with regards to missing documentation and we all can
> try to remedy it.

The problem is that there's a lot of problems with the documentation.
I only pointed out the problems with just the first example. I
shouldn't be the one to fix the issues since I only have a partial
understanding of the examples and I don't want to introduce errors.


>
> We have a Vector class that, by design, was not intended to be generated
> explicitly. The approach is to create a reference frame and specify vectors
> using the basis of the reference frame.
>
> It would be relatively easy to add vector creation functions/methods as you
> suggest that interact with the Vector class, giving a different view on the
> creation and manipulation of vectors.

I don't have a problem with creating a reference frame, it's just I
can't see why one couldn't just pass a reference frame to a vector
constructor along with components that multiply the basis vectors.

> v2pt and a2pt are helper methods to define velocities of points on the same
> rigid body. You don't have to use them if you don't want to.
>

I don't have a problem with them (except the names), but the example
doesn't explain what's happening when they're being called. Why the
_theory suffix, it seems unnecessary?

> We'd love to see some pull requests implementing features that you've
> suggested. We'd also gladly fix and improve and documentation errors or
> clarity issues that you find. The easy way to keep track of that would be
> submitting them as issues.

It will take some time, but I'll try and put together some issues to
help fix the documentation.

Cheers,

Gilbert gede

unread,
Apr 7, 2012, 1:58:54 AM4/7/12
to sympy
Tim,

I can try and explain why a few of these things are the way they are.

As we have learned to deal with inertias for multibody systems, it has
made the most sense to use dyadics, as they are reference frame
independent (as Jason stated).
I understand your confusion with the RigidBody example not showing the
inertia function, changing this probably could help to make things
more clear. I would also recommend that you look through the
documentation online: http://docs.sympy.org/dev/modules/physics/mechanics/index.html

The implementation of the Parallel Axis Theorem is indeed hidden
inside the Kane class; there are multiple reasons for this. One is
that the inertia of a rigid body will obviously be different around
different points, as the mr^2 term has to be added. When dealing with
a dyadic which represents the inertia, we don't have the mass of the
associated rigid body, and cannot compute the new inertia dyadic. By
putting this within Kane, this operation happens at a point in time
where we have both the mass and inertia dyadic for the rigid body. If
you want to do it outside of Kane, you can either repurpose the code
which is within Kane, or a function could be written. At the time, I
decided that if people were going to be doing operations on inertia
dyadics they would  be comfortable doing the dyadic manipulations in
their script, or that they would have a known inertia dyadic about
another point, did not want to have to think about this, and could
just supply things that way.

The kinematic_equations function is used to relate whatever
generalized speeds and generalized coordinates you are using which are
involved in the orientation of reference frames. They are indeed hard
coded for a variety of rotation types, mainly those using Euler angles
or Euler parameters. Obviously, more work could be done here;
Rodrigues parameters are obviously not included but could  be, or
rotations about an arbitrary axis. The point of this function was to
remove the need for the user to look up these kinematic differential
equations; more functionality could be added if desired.

dyanmicssymbols simply produces undefined functions of time. You could
certainly make your own undefined functions of time and coordinates.

I also have to admit I'm not sure which examples you have been looking
at; letting me know which ones have deficiencies allows me to improve
them.

I understand that the implementation of everything in this module
might not be to your taste.

I would, however, implore you to do two things to help with this, if
you would like to see improvements: submit issues for the problems
you've found, as Jason said, where things are not clear or need more
explanation; and please share (in electronic form) some references and
examples which illustrate the ways in which you do dynamics and in
which you have taught dynamics.

I do not think simply ignoring everything in here and writing your own
approach will be very a good use of your time or will be beneficial
for everyone. It does not give us a chance to improve what is here so
it can be more accessible to more people.

-Gilbert




On Apr 6, 9:56 pm, Tim Lahey <tim.la...@gmail.com> wrote:
> Jason,
>
> Responses interleaved below,
>
> >https://github.com/angadhn/sympy/blob/a368bc0c871a75cbea0ffe1a2efc553...

Jason Moore

unread,
Apr 7, 2012, 2:00:19 AM4/7/12
to sy...@googlegroups.com
Tim,

More interspersed below:


On Friday, April 6, 2012 9:56:15 PM UTC-7, Tim Lahey wrote:
Jason,

Responses interleaved below,

On Fri, Apr 6, 2012 at 8:15 PM, Jason Moore <moore...@gmail.com> wrote:
>
> I would beg differ. The equation of motion derivation is completely separate
> from the kinematic and kinetic definitions. As is stands, any method of
> equation of motion derivation can be used to solve classical multi-body
> dynamics with the currently defined classes. We haven't written any examples
> for methods other than Kane's yet due to sparse resources to do so. We also
> haven't implemented automated classes for any other methods yet, for the
> same reason. One GoSC code applicant hopes to remedy both of those issues
> this summer.


True, you can use any approach for the kinematic definitions and still
derive the equations of motion using a different approach. That said,
what I'm saying is that the mechanics package's approach for the
kinematics is very different than what one would use if it was based
upon a different text. The approach to a vector illustrates this. The
approach I learnt for advanced dynamics embeds the frame as part of
the vector. Both approaches aren't wrong, just different.

The Vector class which is implemented does embed frame(s) as part of the vector. If you check out the Vector class design you will see that it contains the both scalars and reference frames. Keep in mind that Vectors can be defined with respect to multiple frames and the current formulation allows that.
 

The problem I have is that a lot of decisions have been fixed so
there's little flexibility. I'll probably have to ignore the whole
thing and write my own. For example, the decision representing the
inertia as a dyad rather than a matrix or a tensor. By choosing a
dyad, I can't just pull the matrix out and use matrix algebra tools
with it (e.g., matrix multiplication for rotations).

See below.
 

The same problem
occurs with the use of "dynamicsymbol" which assumes that the symbol
is only a function of time. This precludes flexible bodies. This alone
means I can't use it.

I feel like the dynamic symbol could surely be expanded to be a function of space and time, ultimately allowing for flexible bodies. Or some other spatial symbol introduced. I'm sure there is a solution that will work within the framework we've set up. But the reality is that it was not designed from the ground up for flexible bodies and that may be an issue. Probably not surmountable though.
 

>
> We have a dyadic class to handle inertia. Dyadics are typically not as
> popular as tensor formulations but they have the advantage of being basis
> independent (unlike tensors). We probably need more functionality to expose
> the dyadics as tensors so that they folks who don't want to learn about
> dyadics may have an easier time working with them. The current inertia
> function is one such function that does that. It allows you to define an
> inertia dyadic with a tensor like formulation.

I understand about inertia dyads, but they have limitations of their own.

>
> We don't have an explicit parallel axis method or function. Here is an
> example of how you can write the parallel axis function with dyadic
> notation:
> https://github.com/angadhn/sympy/blob/a368bc0c871a75cbea0ffe1a2efc553fced1cbd6/sympy/physics/mechanics/kane.py#L471
>
> We'd love for someone to formalize this into some functions that follow the
> more traditional ideas for ease of understanding.

This is one of the problems of using an inertia dyad. Most texts will
present an inertia matrix (or tensor) and give the parallel axis
theorem as a matrix. Plus, rotation of the inertia matrix can be
handled easily using rotation matrices. Given good matrix support,
this can be more efficient.

I don't think the the current inertia structure precludes this at all. There is nothing stopping anyone from writing a method which returns the inertia matrix from the Inertia object with respect to a reference frame. You can then do all the matrix operations on it that you want.

I could imagine some thing like:

I = inertia(N, 1, 2, 3) # create inertia dyadic
rotI = A.dcm(N).transpose() * I.matrix(N) * A.dcm(N) # rotI is just a matrix

where I.matrix(N) returns the matrix form of the inertia with respect to a given reference frame.

but why do that when you can already simply type:

I.express(A)

for the same result.

or

I.matrix(A)

if we add this new matrix method.

Another option would be to allow the add, subtract, multiply, etc operators for the inertia class interact with a matrix in the way you are asking. But careful attention will be needed by the user to keep track of the which reference frame the inertia is define with respect to as matrices do not contain this information.

Just because we have an inertia dyadic built in that is recognized by the RigidBody class and is inter-operable with reference frames and vectors doesn't necessarily mean you can't work with matrices instead.
 
>

> As with all open source projects, the documentation is sometimes weak.
> Please submit issues with regards to missing documentation and we all can
> try to remedy it.

The problem is that there's a lot of problems with the documentation.
I only pointed out the problems with just the first example. I
shouldn't be the one to fix the issues since I only have a partial
understanding of the examples and I don't want to introduce errors.

Adding new examples is also encouraged. Especially those from others' perspectives.
 


>
> We have a Vector class that, by design, was not intended to be generated
> explicitly. The approach is to create a reference frame and specify vectors
> using the basis of the reference frame.
>
> It would be relatively easy to add vector creation functions/methods as you
> suggest that interact with the Vector class, giving a different view on the
> creation and manipulation of vectors.

I don't have a problem with creating a reference frame, it's just I
can't see why one couldn't just pass a reference frame to a vector
constructor along with components that multiply the basis vectors.

This can be done, it would just be a simple case of the general vector construction that we have implemented. A vector function could be written to do that. For example:

v = vector(frameA, (a, b, c)) # would create a vector in frame A with a, b and c as the i, j, k measure numbers.

Keep in mind that our more general formulation allows one to create vectors from components defined in as many reference frames as you want.

v = a * N.x + b * A.y + c *H.z + d * N.y
 

> v2pt and a2pt are helper methods to define velocities of points on the same
> rigid body. You don't have to use them if you don't want to.
>

I don't have a problem with them (except the names), but the example
doesn't explain what's happening when they're being called. Why the
_theory suffix, it seems unnecessary?

No idea. At one point it didn't have the theory part on it. I'm not sure why it changed. You can git blame the code to see who changed it.
 

> We'd love to see some pull requests implementing features that you've
> suggested. We'd also gladly fix and improve and documentation errors or
> clarity issues that you find. The easy way to keep track of that would be
> submitting them as issues.

It will take some time, but I'll try and put together some issues to
help fix the documentation.

Awesome.
 

Tim Lahey

unread,
Apr 7, 2012, 2:30:03 AM4/7/12
to sy...@googlegroups.com
Hi,

I've responded below.

On Sat, Apr 7, 2012 at 1:58 AM, Gilbert gede <gilbe...@gmail.com> wrote:
> Tim,
>
> I can try and explain why a few of these things are the way they are.
>
> As we have learned to deal with inertias for multibody systems, it has
> made the most sense to use dyadics, as they are reference frame
> independent (as Jason stated).
> I understand your confusion with the RigidBody example not showing the
> inertia function, changing this probably could help to make things
> more clear. I would also recommend that you look through the
> documentation online: http://docs.sympy.org/dev/modules/physics/mechanics/index.html

I've been looking at the documentation, but part of the problem is
that documentation assumes a lot on the part of the reader.

>
> The implementation of the Parallel Axis Theorem is indeed hidden
> inside the Kane class; there are multiple reasons for this. One is
> that the inertia of a rigid body will obviously be different around
> different points, as the mr^2 term has to be added. When dealing with
> a dyadic which represents the inertia, we don't have the mass of the
> associated rigid body, and cannot compute the new inertia dyadic. By
> putting this within Kane, this operation happens at a point in time
> where we have both the mass and inertia dyadic for the rigid body. If
> you want to do it outside of Kane, you can either repurpose the code
> which is within Kane, or a function could be written. At the time, I
> decided that if people were going to be doing operations on inertia
> dyadics they would  be comfortable doing the dyadic manipulations in
> their script, or that they would have a known inertia dyadic about
> another point, did not want to have to think about this, and could
> just supply things that way.

Having it as part of the Kane class is a bad idea since it leads to
duplication of code. I don't think it makes sense to have an inertia
outside of a rigid body. If you contain it in a rigid body, then you
have all the information you need. Plus, if you associate it with a
frame and not just a point, you can handle the rotation of the inertia
to a new frame or translation to another frame by telling it to switch
the inertia to a different frame. Based upon the relationship between
the frames, you'd apply parallel axis theorem or rotations as
appropriate.

>
> The kinematic_equations function is used to relate whatever
> generalized speeds and generalized coordinates you are using which are
> involved in the orientation of reference frames. They are indeed hard
> coded for a variety of rotation types, mainly those using Euler angles
> or Euler parameters. Obviously, more work could be done here;
> Rodrigues parameters are obviously not included but could  be, or
> rotations about an arbitrary axis. The point of this function was to
> remove the need for the user to look up these kinematic differential
> equations; more functionality could be added if desired.

I was wondering about the hard-coding aspect rather than deriving them
directly. It just seemed like a lot of unnecessary work over a
derivation approach. Plus, since the symbols aren't explained anywhere
in the docstring, it's a bit hard to read.

>
> dyanmicssymbols simply produces undefined functions of time. You could
> certainly make your own undefined functions of time and coordinates.

It's not clear how to do that since dynamicsymbol is kind of
throughout the code.

>
> I also have to admit I'm not sure which examples you have been looking
> at; letting me know which ones have deficiencies allows me to improve
> them.

I was looking at the Rolling Disc first example.

>
> I understand that the implementation of everything in this module
> might not be to your taste.
>
> I would, however, implore you to do two things to help with this, if
> you would like to see improvements: submit issues for the problems
> you've found, as Jason said, where things are not clear or need more
> explanation; and please share (in electronic form) some references and
> examples which illustrate the ways in which you do dynamics and in
> which you have taught dynamics.

The text I learnt dynamics from doesn't exist in electronic (or
printed form) yet. My supervisor and his co-author are working on it
for Cambridge University Press (I think). They originally learnt it at
UTIAS (University of Toronto Institute for Aerospace Studies). It's
well suited to computational implementation and has nice compact
notation. Another grad student I know implemented it using Matlab's
symbolic toolbox. That said, the approach the mechanics module is
using is different from a lot of texts I've read.

>
> I do not think simply ignoring everything in here and writing your own
> approach will be very a good use of your time or will be beneficial
> for everyone. It does not give us a chance to improve what is here so
> it can be more accessible to more people.

The problem is that a lot of assumptions are built-in (e.g., the
dynamics symbol business and how inertia is handled) and it's likely
faster for me to just write my own rather than work around them. I'm
not going to do it in the short term since I already have a working
package in Maple that meets my needs. It's not as automatic, but I
prefer that so I don't have to worry as much about possible errors in
a code black box. I want to peek inside to make sure everything is
correct.

Tim Lahey

unread,
Apr 7, 2012, 2:48:31 AM4/7/12
to sy...@googlegroups.com
Jason,

I've responded below.

On Sat, Apr 7, 2012 at 2:00 AM, Jason Moore <moore...@gmail.com> wrote:
>
> The Vector class which is implemented does embed frame(s) as part of the
> vector. If you check out the Vector class design you will see that it
> contains the both scalars and reference frames. Keep in mind that Vectors
> can be defined with respect to multiple frames and the current formulation
> allows that.
>

Defining a vector with respect to multiple frames seems like a bad
idea and could lead to errors. I could see how one might want to do
that, though. I'd do it as a sum of two vectors each defined in the
different frames though.

>
> I feel like the dynamic symbol could surely be expanded to be a function of
> space and time, ultimately allowing for flexible bodies. Or some other
> spatial symbol introduced. I'm sure there is a solution that will work
> within the framework we've set up. But the reality is that it was not
> designed from the ground up for flexible bodies and that may be an issue.
> Probably not surmountable though.
>

It might be possible, by the dynamic symbol is kind of throughout the
code so I expect that any solution will require a fair bit of work to
extend it.

>
> I don't think the the current inertia structure precludes this at all. There
> is nothing stopping anyone from writing a method which returns the inertia
> matrix from the Inertia object with respect to a reference frame. You can
> then do all the matrix operations on it that you want.
>
> I could imagine some thing like:
>
> I = inertia(N, 1, 2, 3) # create inertia dyadic
> rotI = A.dcm(N).transpose() * I.matrix(N) * A.dcm(N) # rotI is just a matrix
>
> where I.matrix(N) returns the matrix form of the inertia with respect to a
> given reference frame.
>
> but why do that when you can already simply type:
>
> I.express(A)
>
> for the same result.
>
> or
>
> I.matrix(A)
>
> if we add this new matrix method.

And if I want to rotate the inertia matrix? How does one accomplish
this with the current setup?

>
> Another option would be to allow the add, subtract, multiply, etc operators
> for the inertia class interact with a matrix in the way you are asking. But
> careful attention will be needed by the user to keep track of the which
> reference frame the inertia is define with respect to as matrices do not
> contain this information.
>

True, matrices don't contain the reference frames. But matrices work
well for handling the parallel axis theorem or rotation to a different
frame behind the scenes as well. Personally, I do keep track of what
frame each thing is in. One of the problems I have with the current
module is that it's hard to see what's happening since the module is
designed to hide the complexity of the derivation. I'd like to peek at
the details to verify things are correct.

>
> Adding new examples is also encouraged. Especially those from others'
> perspectives.

I can't exactly add new examples when I don't know what the module is
doing or how to do various things.

>
> This can be done, it would just be a simple case of the general vector
> construction that we have implemented. A vector function could be written to
> do that. For example:
>
> v = vector(frameA, (a, b, c)) # would create a vector in frame A with a, b
> and c as the i, j, k measure numbers.
>
> Keep in mind that our more general formulation allows one to create vectors
> from components defined in as many reference frames as you want.
>
> v = a * N.x + b * A.y + c *H.z + d * N.y
>

I'd consider that to be the sum of three vectors in different
reference frames (one in frame N, one in frame A and another in frame
H). It's a big mess as it's written there and I'm sure I'd make a
mistake. I certainly work with things like this, but I always try to
convert to a base frame first. I'm assuming it's happening behind the
scenes here.

Cheers,

Tim.
--
Tim Lahey
tjl...@cgl.uwaterloo.ca

Gilbert gede

unread,
Apr 7, 2012, 5:47:17 AM4/7/12
to sympy
Inertia dyadics are fixed to the frame they are originally defined in;
you cannot rotate a fixed vector/dyadic in its frame, but you can
rotate the frame they are fixed in relative to another frame. What you
can then do, for both dyadics and vectors, is express them in this
different rotated frame.

This shows how to get the inertia matrix for an inertia dyadic,
rotated into whatever frame you want.

In [1]: from sympy import *
In [2]: from sympy.physics.mechanics import *
In [3]: q1 = dynamicsymbols('q1')
In [4]: Ixx, Iyy, Izz = symbols('Ixx Iyy Izz')
In [5]: N = ReferenceFrame('N')
In [6]: A = N.orientnew('A', 'Axis', [q1, N.x])
In [7]: I = inertia(N, Ixx, Iyy, Izz)
In [8]: rotI = lambda I, f: Matrix([j & I & i for i in f for j in
f]).reshape(3,3)
In [9]: rotI(I, N)
Out[9]:
[Ixx, 0, 0]
[ 0, Iyy, 0]
[ 0, 0, Izz]
In [10]: rotI(I, A)
Out[10]:
[Ixx,
0, 0]
[ 0, Iyy*cos(q1(t))**2 + Izz*sin(q1(t))**2, -
Iyy*sin(q1(t))*cos(q1(t)) + Izz*sin(q1(t))*cos(q1(t))]
[ 0, -Iyy*sin(q1(t))*cos(q1(t)) +
Izz*sin(q1(t))*cos(q1(t)), Iyy*sin(q1(t))**2 +
Izz*cos(q1(t))**2]
Reply all
Reply to author
Forward
0 new messages