Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP

Loading…

LagrangesMethod class for sympy.physics.mechanics. #1460

Merged
merged 6 commits into from

9 participants

@angadhn

The docstrings will likely need to be revamped.

@smichr smichr commented on the diff
sympy/physics/mechanics/lagrange.py
((13 lines not shown))
+ arguments. The Lagrangian multipliers are automatically generated and are
+ equal in number to the constraint equations.Similarly if there are any
+ non-conservative forces, they can be supplied in a list along with a
+ ReferenceFrame. This is discussed further in the __init__ method.
+
+ Attributes
+ ==========
+
+ mass_matrix : Matrix
+ The system's mass matrix
+
+ forcing : Matrix
+ The system's forcing vector
+
+ mass_matrix_full : Matrix
+ The "mass matrix" for the qdot's, qdoubledot's, and the
@smichr Collaborator
smichr added a note

in my experience, mixing single and double quotes in a docstring generates sphinx errors...let's see if this happens with your request.

@angadhn
angadhn added a note

This hasn't been an issue in the "Kane" class, so I presumed it wouldn't be a problem here but that could've been an erroneous assumption.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/physics/mechanics/lagrange.py
((64 lines not shown))
+
+ self._L = sympify(Lagrangian)
+ self.eom = None #initializing the eom Matrix
+ self._m_cd = Matrix([]) #Mass Matrix of differentiated coneqs
+ self._m_d = Matrix([]) #Mass Matrix of dynamic equations
+ self._f_cd = Matrix([]) #Forcing part of the diff coneqs
+ self._f_d = Matrix([]) #Forcing part of the dynamic equations
+ self.lam_coeffs = Matrix([]) #Initializing the coeffecients matrix of lams
+
+ self.forcelist = forcelist
+ self.inertial = frame
+
+ self.lam_vec = Matrix([])
+
+
+ #What used to be the coords method
@smichr Collaborator
smichr added a note

if it's not there anymore then this comment will have no meaning for anyone looking at this in the future, so it can either be deleted or changed to indicate what you are doing, perhaps?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/physics/mechanics/lagrange.py
((73 lines not shown))
+ self.forcelist = forcelist
+ self.inertial = frame
+
+ self.lam_vec = Matrix([])
+
+
+ #What used to be the coords method
+
+ q_list = list(q_list)
+ if not isinstance(q_list, list):
+ raise TypeError('Generalized coords. must be supplied in a list')
+ self._q = q_list
+ self._qdots = [diff(i, dynamicsymbols._t) for i in self._q]
+ self._qdoubledots = [diff(i, dynamicsymbols._t) for i in self._qdots]
+
+ #What used to be the constraints method
@smichr Collaborator
smichr added a note

ditto about what used to be

@angadhn
angadhn added a note

Thanks! I forgot to edit those lines. I had left those in as notes for myself.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/physics/mechanics/lagrange.py
((169 lines not shown))
+ raise TypeError('First entry in force pair is a point'
+ ' or frame')
+
+ else:
+ term4 = zeros(n,1)
+
+ self.eom = term1 - term2 -term3 -term4
+
+ #The mass matrix is generated by the following
+ self._m_d = (self.eom).jacobian(qdd)
+
+ return self.eom
+
+ @property
+ def mass_matrix(self):
+ # Returns the mass matrix, which is augmented by the Lagrange
@smichr Collaborator
smichr added a note

convert to docstring instead of comments

@angadhn
angadhn added a note

Once again, I was replicating the presentation of "Kane". I have changed this now though. Thanks.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/physics/mechanics/lagrange.py
((187 lines not shown))
+ # an n X n matrix is returned.
+ # If there are 'n' generalized coordinates and 'm' constraint equations
+ # have been supplied during initialization then an n X (n+m) matrix is
+ # returned. The (n + m - 1)th and (n + m)th columns contain the
+ # coeffecients of the lagrange multipliers.
+
+ if self.eom == None:
+ raise ValueError('Need to compute the equations of motion first')
+ if len(self.lam_coeffs) != 0:
+ return (self._m_d).row_join((self.lam_coeffs).transpose())
+ else:
+ return self._m_d
+
+ @property
+ def mass_matrix_full(self):
+ n = len(self._q)
@smichr Collaborator
smichr added a note

docstring?

@angadhn
angadhn added a note

The docstring I had written was a little contorted so I decided to leave it blank. I have added a docstring now though. Just waiting for the decision on the '*'/PEP 8 comments to push again.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@smichr smichr commented on the diff
sympy/physics/mechanics/lagrange.py
((205 lines not shown))
+ #THE FIRST TWO ROWS OF THE MATRIX
+ row1 = eye(n).row_join(zeros(n,n))
+ row2 = zeros(n,n).row_join(self.mass_matrix)
+ if self.coneqs != None:
+ m = len(self.coneqs)
+ I = eye(n).row_join(zeros(n,n+m))
+ below_eye = zeros(n+m,n)
+ A = (self.mass_matrix).col_join((self._m_cd).row_join(zeros(m,m)))
+ below_I = below_eye.row_join(A)
+ return I.col_join(below_I)
+ else:
+ A = row1.col_join(row2)
+ return A
+
+ @property
+ def forcing(self):
@smichr Collaborator
smichr added a note

docstring?

@angadhn Check the formatting for single line docstrings that shows up elsewhere and try to match that:
"""Some descriptive sentence."""

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@smichr smichr commented on the diff
sympy/physics/mechanics/lagrange.py
((227 lines not shown))
+
+ if self.coneqs != None:
+ lam = self.lam_vec
+ lamzero = dict(zip(lam, [0] * len(lam)))
+
+ #The forcing terms from the eoms
+ self._f_d = -((self.eom).subs(qddzero)).subs(lamzero)
+
+ else:
+ #The forcing terms from the eoms
+ self._f_d = -(self.eom).subs(qddzero)
+
+ return self._f_d
+
+ @property
+ def forcing_full(self):
@smichr Collaborator
smichr added a note

docstring?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/physics/mechanics/lagrange.py
((235 lines not shown))
+ else:
+ #The forcing terms from the eoms
+ self._f_d = -(self.eom).subs(qddzero)
+
+ return self._f_d
+
+ @property
+ def forcing_full(self):
+ if self.eom == None:
+ raise ValueError('Need to compute the equations of motion first')
+ if self.coneqs != None:
+ return (Matrix(self._qdots)).col_join((self.forcing).col_join(self._f_cd))
+ else:
+ return (Matrix(self._qdots)).col_join(self.forcing)
+
+ def rhs(self, method):
@smichr Collaborator
smichr added a note

convert comments to docstring

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@smichr smichr commented on the diff
sympy/physics/mechanics/lagrange.py
((211 lines not shown))
+ below_eye = zeros(n+m,n)
+ A = (self.mass_matrix).col_join((self._m_cd).row_join(zeros(m,m)))
+ below_I = below_eye.row_join(A)
+ return I.col_join(below_I)
+ else:
+ A = row1.col_join(row2)
+ return A
+
+ @property
+ def forcing(self):
+ # Returns the forcing vector
+ if self.eom == None:
+ raise ValueError('Need to compute the equations of motion first')
+
+ qdd = self._qdoubledots
+ qddzero = dict(zip(qdd, [0] * len(qdd)))
@smichr Collaborator
smichr added a note

PEP8, I believe, would suggest removing spaces around *

@moorepants Collaborator

Are you sure?

http://www.python.org/dev/peps/pep-0008/#other-recommendations

Seems like spaces around these operators is preferred.

@certik Owner
certik added a note

Yes, it seems to me as well, that PEP8 says you should have spaces around "*". In either case, I would use whatever looks more readable on the case by case basis, in this PR it seems that having spaces is more readable.

@asmeurer Owner
asmeurer added a note

We usually differ from PEP 8 in SymPy for * and ** because it makes reading expressions like 3*x**2 + 2*x**4 easier. Here I would say it's justified, though really either way looks fine.

@smichr Collaborator
smichr added a note

I reread PEP8 and see that I was wrong about the * -- I was running under the sympy-modifed version (that Aaron cites). So although we are doing exactly what PEP8 says not to do I think it looks better, too.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/physics/mechanics/lagrange.py
((158 lines not shown))
+ ' or tuples')
+ term4 = zeros(len(qd), 1)
+ for i,v in enumerate(qd):
+ for j,w in enumerate(forcelist):
+ if isinstance(w[0], ReferenceFrame):
+ speed = w[0].ang_vel_in(N)
+ term4[i] += speed.diff(v, N) & w[1]
+ if isinstance(w[0], Point):
+ speed = w[0].vel(N)
+ term4[i] += speed.diff(v, N) & w[1]
+ else:
+ raise TypeError('First entry in force pair is a point'
+ ' or frame')
+
+ else:
+ term4 = zeros(n,1)
@smichr Collaborator
smichr added a note

space afater comma

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/physics/mechanics/lagrange.py
((119 lines not shown))
+ if not isinstance(coneqs, list):
+ raise TypeError('Enter the constraint equations in a list')
+
+ o = len(coneqs)
+
+ #Creating the multipliers
+ self.lam_vec = Matrix(dynamicsymbols('lam1:' + str(o+1)))
+
+ #Extracting the coeffecients of the multipliers
+ coneqs_mat = Matrix(coneqs)
+ qd = self._qdots
+ self.lam_coeffs = -coneqs_mat.jacobian(qd)
+
+ #Determining the third term in Lagrange's EOM
+ #term3 = ((self.lam_vec).transpose() * self.lam_coeffs).transpose()
+ term3 = self.lam_coeffs.transpose() * self.lam_vec
@smichr Collaborator
smichr added a note

remove here (and throughout) spaces around * as per PEP8

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/physics/mechanics/lagrange.py
((110 lines not shown))
+ #Determining the second term in Lagrange's EOM
+ term2 = (L.jacobian(q)).transpose()
+
+ #term1 and term2 will be there no matter what so leave them as they are
+
+ if self.coneqs != None:
+ coneqs = self.coneqs
+ #If there are coneqs supplied then the following will be created
+ coneqs = list(coneqs)
+ if not isinstance(coneqs, list):
+ raise TypeError('Enter the constraint equations in a list')
+
+ o = len(coneqs)
+
+ #Creating the multipliers
+ self.lam_vec = Matrix(dynamicsymbols('lam1:' + str(o+1)))
@smichr Collaborator
smichr added a note

add space around + in o+1

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@travisbot

This pull request fails (merged b3a58b0a into 9625918).

@Krastanov
Collaborator

SymPy Bot Summary: :eight_spoked_asterisk: All tests have passed.

Test command: setup.py test
master hash: 9625918
branch hash: b3a58b0a0b6efe5025feedf214575df54ef9b002

Interpreter 1: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/local/bin/python2.5 (2.5.6-final-0)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYtLYiDA

Interpreter 2: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python2.7 (2.7.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYxqMjDA

Interpreter 3: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python3.2 (3.2.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYs7YiDA

Build HTML Docs: :eight_spoked_asterisk: All tests have passed.

Docs build command: make html-errors
Sphinx version: 1.1.3

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYq48iDA

Automatic review by SymPy Bot.

@certik
Owner

Thanks for opening the pull request. So the first priority is regular tests. Either add them to some file in sympy/physics/mechanics/tests, or create a new one. Then some doctests showing how to use it in the docstring of the various methods and also what Chris mentioned above.

The regular tests are the most important.

@moorepants moorepants commented on the diff
sympy/physics/mechanics/lagrange.py
@@ -0,0 +1,253 @@
+__all__ = ['LagrangesMethod']
+
+from sympy import diff, zeros, Matrix, eye, sympify
+from sympy.physics.mechanics import (dynamicsymbols, ReferenceFrame, Point)
+
+class LagrangesMethod(object):
@moorepants Collaborator

the Kane method class isn't called KanesMethod, so why do that here? It's just more to type. There are no other naming conflicts in this module are there?

@angadhn
angadhn added a note

So, @hazelnusse, @gilbertgede and I talked about this in the lab yesterday. We talked about how other methods of deriving EOMs are referred to i.e. 'Kane's method', 'Hamilton's method', etc and how maybe the "Kane" class could be renamed to "KanesMethod" possibly. Also, just calling it "Lagrange" might be ambiguous.

@moorepants Collaborator

Ok great, we should change the name of the Kane Class too then.

@moorepants I will take care of that in another PR that I am going to open soon.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@travisbot

This pull request fails (merged 8c26b504 into 9625918).

@Krastanov
Collaborator

SymPy Bot Summary: :red_circle: There were test failures.

@angadhn: Please fix the test failures.

Test command: setup.py test
master hash: 7d92368
branch hash: 8c26b50426c7150788f999ca2b17a2db9cca4c24

Interpreter 1: :red_circle: There were test failures.

Interpreter: /usr/local/bin/python2.5 (2.5.6-final-0)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYsKsjDA

Interpreter 2: :red_circle: There were test failures.

Interpreter: /usr/bin/python2.7 (2.7.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYsd0iDA

Interpreter 3: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python3.2 (3.2.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYr6sjDA

Build HTML Docs: :red_circle: There were test failures.

Docs build command: make html-errors
Sphinx version: 1.1.3

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY580iDA

Automatic review by SymPy Bot.

@travisbot

This pull request fails (merged 3092f5e1 into 7d92368).

@Krastanov
Collaborator

SymPy Bot Summary: :red_circle: There were test failures.

@angadhn: Please fix the test failures.

Test command: setup.py test
master hash: 1627b32
branch hash: 3092f5e14af92de40c7c14496902ebfe46f15641

Interpreter 1: :red_circle: There were test failures.

Interpreter: /usr/local/bin/python2.5 (2.5.6-final-0)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYz6MjDA

Interpreter 2: :red_circle: There were test failures.

Interpreter: /usr/bin/python2.7 (2.7.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYle0iDA

Interpreter 3: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python3.2 (3.2.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYoIwjDA

Build HTML Docs: :eight_spoked_asterisk: All tests have passed.

Docs build command: make html-errors
Sphinx version: 1.1.3

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYjvAhDA

Automatic review by SymPy Bot.

@travisbot

This pull request fails (merged 619dcc6d into 65b6582).

@travisbot

This pull request fails (merged 84a5ea6f into 65b6582).

@Krastanov
Collaborator

SymPy Bot Summary: :red_circle: There were test failures.

@angadhn: Please fix the test failures.

Test command: setup.py test
master hash: 98cc80f
branch hash: 84a5ea6f2899b4dbce6df52b3c0b6f01ad4497f9

Interpreter 1: :red_circle: There were test failures.

Interpreter: /usr/local/bin/python2.5 (2.5.6-final-0)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY8P8hDA

Interpreter 2: :red_circle: There were test failures.

Interpreter: /usr/bin/python2.7 (2.7.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYjbMjDA

Interpreter 3: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python3.2 (3.2.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY7_8hDA

Build HTML Docs: :eight_spoked_asterisk: All tests have passed.

Docs build command: make html-errors
Sphinx version: 1.1.3

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYu6sjDA

Automatic review by SymPy Bot.

sympy/physics/mechanics/tests/test_lagrange.py
((102 lines not shown))
+ # of the string fixed frame is redefined as 'LagrangesMethod' doesn't require
+ # generalized speeds, per se. (Lagrangian mechanics requires 'simple'
+ # generalized speeds)
+ A.orientnew('A', 'Axis', [q, N.z])
+ A.set_ang_vel(N, qd *A.z)
+ P.v2pt_theory(O,N,A)
+ T = 1/2.0 * m * P.vel(N) & P.vel(N) # T is the kinetic energy of the system
+ V = - m * g * l * cos(q) # V is the potential energy of the system
+ L = T - V # L is the Lagrangian
+
+ # The 'LagrangesMethod' class is invoked and the equations of motion are generated.
+ l = LagrangesMethod(L, [q])
+ l.lagranges_equations()
+ l.rhs("GE")
+
+ # Finally, we the results of both methods and it's seen that they are
@certik Owner
certik added a note

the word "compare" is missing?

@angadhn
angadhn added a note

thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/physics/mechanics/lagrange.py
((22 lines not shown))
+ The system's mass matrix
+
+ forcing : Matrix
+ The system's forcing vector
+
+ mass_matrix_full : Matrix
+ The "mass matrix" for the qdot's, qdoubledot's, and the
+ lagrange multipliers (lam)
+
+ forcing_full : Matrix
+ The forcing vector for the qdot's, qdoubledot's and
+ lagrange multipliers (lam)
+
+ rhs : Matrix
+ Solves for the states (i.e. q's, qdot's, and multipliers)
+
@certik Owner
certik added a note

Can you put a simple doctest here showing how to use the class?

@angadhn
angadhn added a note

I've added it now.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@certik
Owner

You seem to have some trailing whitespace problems there: http://travis-ci.org/#!/sympy/sympy/jobs/2056371

sympy/physics/mechanics/lagrange.py
((75 lines not shown))
+
+ self.lam_vec = Matrix([])
+
+
+ # Creating the qs, qdots and qdoubledots
+
+ q_list = list(q_list)
+ if not isinstance(q_list, list):
+ raise TypeError('Generalized coords. must be supplied in a list')
+ self._q = q_list
+ self._qdots = [diff(i, dynamicsymbols._t) for i in self._q]
+ self._qdoubledots = [diff(i, dynamicsymbols._t) for i in self._qdots]
+
+ self.coneqs = coneqs
+
+ def lagranges_equations(self):

Maybe form_lagranges_equations instead?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
sympy/physics/mechanics/lagrange.py
((130 lines not shown))
+ self.lam_coeffs = -coneqs_mat.jacobian(qd)
+
+ #Determining the third term in Lagrange's EOM
+ #term3 = ((self.lam_vec).transpose() * self.lam_coeffs).transpose()
+ term3 = self.lam_coeffs.transpose() * self.lam_vec
+
+ #Taking the time derivative of the constraint equations
+ diffconeqs = [diff(i, dynamicsymbols._t) for i in coneqs]
+
+ #Extracting the coeffecients of the qdds from the diff coneqs
+ diffconeqs_mat = Matrix(diffconeqs)
+ qdd = self._qdoubledots
+ self._m_cd = diffconeqs_mat.jacobian(qdd)
+
+ #The remaining terms i.e. the 'forcing' terms in diff coneqs
+ qddzero = dict(zip(qdd, [0] * len(qdd)))

Isn't len(qdd) equal to n, which you use elsewhere?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@angadhn

So it appears that I have mistakenly rebased a couple of my branches and so my commits (including the most recent one on this page) aren't showing up correctly. Any way I can fix this?

@angadhn

Actually, it's fine here. Never mind.

@travisbot

This pull request fails (merged 758c380b into 0d88b25).

@certik
Owner

You still have a trailing whitespace in lagrange.py line 49, see the Travis test results:

AssertionError: File contains trailing whitespace: /home/vagrant/virtualenv/python2.6/lib/python2.6/site-packages/sympy/physics/mechanics/lagrange.py, line 49.
@Krastanov
Collaborator

SymPy Bot Summary: :red_circle: There were test failures.

@angadhn: Please fix the test failures.

Test command: setup.py test
master hash: 0d88b25
branch hash: 758c380b3b2d871776782ceb5da985b936e58fb4

Interpreter 1: :red_circle: There were test failures.

Interpreter: /usr/local/bin/python2.5 (2.5.6-final-0)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYwqsjDA

Interpreter 2: :red_circle: There were test failures.

Interpreter: /usr/bin/python2.7 (2.7.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYou0iDA

Interpreter 3: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python3.2 (3.2.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY-JsjDA

Build HTML Docs: :eight_spoked_asterisk: All tests have passed.

Docs build command: make html-errors
Sphinx version: 1.1.3

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYkZQjDA

Automatic review by SymPy Bot.

@travisbot

This pull request passes (merged 9d7d8b3b into 0d88b25).

@travisbot

This pull request passes (merged e4d9026f into 0d88b25).

@moorepants
Collaborator

@certik @Krastanov

So this pull request may take a while to get in. It implements a complicated method and the testing is still weak. @angadhn pulled earlier than usual so that we can comment on the code easily and discuss some of the issues. Is there anyway not to get bombarded by the bots constantly? I personally find it stressful. We probably only need the bot testing when the PR is closer to being mergable.

@Krastanov
Collaborator

SymPy Bot Summary: :eight_spoked_asterisk: All tests have passed.

Test command: setup.py test
master hash: 0d88b25
branch hash: e4d9026f4f3ad212df580540d3018c28e81fdfda

Interpreter 1: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/local/bin/python2.5 (2.5.6-final-0)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYoPAhDA

Interpreter 2: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python2.7 (2.7.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYjcYiDA

Interpreter 3: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python3.2 (3.2.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY3KMjDA

Build HTML Docs: :eight_spoked_asterisk: All tests have passed.

Docs build command: make html-errors
Sphinx version: 1.1.3

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY5tUiDA

Automatic review by SymPy Bot.

@angadhn

So the docs commit came into Lagrange because I thought that the Lagrange class was complete w.r.t. the tests; they are identical in number to those in Kane (apart from that I will also be adding the rolling disc). The two degree of freedom test in my case is the double pendulum. Like Kane, I also have the simple pendulum. Lagrange's equations don't need the test for auxiliary equations and ditto for the parallel axis theorem.

@travisbot

This pull request passes (merged e5d912a2 into 15c2c89).

@moorepants
Collaborator

Angadh, your current tests don't really cover the use cases of the class. An ideal set of tests would probably include a test for each problem case individually and at least one test case problem that covers a problem that has all dynamics features. This way you know that your class plays well with whatever the user may throw at it. It will also prove that your class gives the correct answer for a robust dynamics problem and if it does that then we would all be convinced that it can probably handle any lesser problems. You'd hate to send out a class into the world that incorrectly derives equations of certain types of problems.

I think you need at least a 3D problem with multiple degrees of freedom, both holonomic and non-holonomic constraints, and analytically unsolvable kinematic loop.

The nice thing is that Gilbert already has problems like this in his Kane test suite. The only thing you have to do is to write out the Lagrangian for those problems and then see if your Lagrange class gives the same answer as Gilbert's. In fact, both Kane and Lagrange classes should be tested by the same problem test suite as they are both just different ways to crack the coconut and should give the exact same answer.

sympy/physics/mechanics/tests/test_lagrange.py
((68 lines not shown))
+ N = ReferenceFrame('N')
+ A = N.orientnew('A', 'Axis', [q, N.z])
+ A.set_ang_vel(N, qd * N.z)
+
+ # Next, we create the point O and fix it in the inertial frame. We then
+ # locate the point P to which the bob is attached. Its corresponding
+ # velocity is then determined by the 'two point formula'.
+ O = Point('O')
+ O.set_vel(N, 0)
+ P = O.locatenew('P', l * A.x)
+ P.v2pt_theory(O, N, A)
+
+ # The 'Particle' which represents the bob is then created.
+ Pa = Particle('Pa', P, m)
+
+ T = 1/2.0 * m * P.vel(N) & P.vel(N) # T is the kinetic energy of the system
@moorepants Collaborator

Whenever you use the & and ^ for the dot and cross products it is probably good practice to always enclose the operation in parentheses because these two operators have special order of operation rules that don't necessarily act like you think they do. In this case it is fine, but I think all example code that a user may possibly look at to figure out how to write their own problems should have this explicit. New users may spend forever trying to figure out the bug if they type something like

A & B ^ C ^ D

and expect the order of operations to match they way they think about it on paper.

Look back on the sympy mailing list, the online docs, or old pull requests to see how these operators work (or you may already know). http://docs.sympy.org/dev/modules/physics/mechanics/vectors.html#vector-algebra-in-mechanics only has a small warning about this.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@travisbot

This pull request passes (merged 85ac0417 into 15c2c89).

sympy/physics/mechanics/lagrange.py
((21 lines not shown))
+
+ mass_matrix : Matrix
+ The system's mass matrix
+
+ forcing : Matrix
+ The system's forcing vector
+
+ mass_matrix_full : Matrix
+ The "mass matrix" for the qdot's, qdoubledot's, and the
+ lagrange multipliers (lam)
+
+ forcing_full : Matrix
+ The forcing vector for the qdot's, qdoubledot's and
+ lagrange multipliers (lam)
+
+ rhs : Matrix
@moorepants Collaborator

This is a method, so you should have a Methods section to accompany your Attributes section. Also, don't forget the form_equations method.

Jason, I don't think that Methods are usually listed in the class docstring, although rhs should be taken out of here.

@moorepants Collaborator

Ok, here is some potential reference material for that:

https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt#documenting-classes

Although, I don't know the standard in SymPy is.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@travisbot

This pull request fails (merged cdaf29ba into 15c2c89).

@gilbertgede

@moorepants I have to disagree with you on some of these points Jason.

We just realized that the Kane class actually returns a negative mass matrix and a negative forcing vector. In all calculations that have been done, these cancel out and it goes unnoticed. Angadh and I discovered it today though. I suggested that he ensure that his results check out against hand/reference calculations rather than what the Kane class gives (also, if any other changes ever break Kane, Lagrange won't fail it's tests).

Not all of the test cases go between the 2 classes either - if the generalized speeds are not just the derivatives of the coordinates, things get messy. That being said, storing some of these reference results that could be used for both cases in 1 place seems like a good idea.

Also, Angadh is working to get reference results for a few more tests. Note that the Kane class does not have a test within SymPy for holonomic and nonholonomic constraints.

@travisbot

This pull request fails (merged 0fc013fe into 15c2c89).

@moorepants
Collaborator

I'm not sure what you are disagreeing with, seems more like agreeing. I'm glad you found the bug in the Kane class.

Yes, results should always check against well known results that are in literature as the de facto test. I'm suggesting that both classes can test to the same benchmark problems from literature. i.e. setup up a benchmark problem with known canonical results and then run the two classes to see if you get the same answers as literature. There seems to be little reason to make different test problems for both of the classes. These test cases certainly need to be designed such that the resulting equations of motion use the same coordinates and such, but they certainly should both be able to solve the same problems and give the same EoMs.

Secondly, maintaining different sets of benchmark tests for many automated EoM derivation methods will ultimately be a nightmare. I can imagine that all methods test to the same benchmarks, but then there may be some special test cases for each Method class that test the particularities and advantages of one method or an another.

I thought you have the Bicycle example as a test for a "complete" dynamics problem. This test could probably be run in the test suite even though it is slower, especially since the tests are mostly run by the bots in the background.

Also a simple four bar linkage or a 3D version of one could be a nice holonomic/kinematic loop test. I think the problem in Kane's online dynamics in Chapter Zero is a good one for that.

Lastly, we should find a problem that is "complete" but which is simpler than the bicycle for computation time reasons. We've talked about this in the past and I vaguely remember suggesting one. I'll dig around for that.

@Krastanov
Collaborator

SymPy Bot Summary: :eight_spoked_asterisk: All tests have passed.

Test command: setup.py test
master hash: 15c2c89
branch hash: 0fc013fe86fdfd38d84a04d1b1d611e4250acac2

Interpreter 1: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/local/bin/python2.5 (2.5.6-final-0)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYxt0iDA

Interpreter 2: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python2.7 (2.7.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYgpwjDA

Interpreter 3: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python3.2 (3.2.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY_P8hDA

Build HTML Docs: :eight_spoked_asterisk: All tests have passed.

Docs build command: make html-errors
Sphinx version: 1.1.3

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYgZwjDA

Automatic review by SymPy Bot.

@moorepants
Collaborator

I asked Mont if he knew of a simpler problem than the bicycle which was full featured and he suggested the unicycle with the riders leg(s). The leg and crank form a four bar linkage. He doubted than anyone had the solution to this problem written down, but it may be possible. We could compute the EoM's using Autolev and then compare the results to what the classes in SymPy give if writing the solution by hand is a mess.

I was think of something similar today too that added a four bar linkage to the the rolling disc, but it doesn't have as nice of a physical interpretation.

I'm willing to make this test up, but can't do it till after he 21st (my dissertation due date).

@hazelnusse
@certik
Owner

@gilbertgede, @moorepants, are you ok with pusing this in? It looks good to me.

@angadhn, for the future, see here how to write good commit messages: https://plus.google.com/u/0/104039945248245758823/posts/GcJmpxCVDfg
So for example instead of "More edits." it's better to explain what exactly you did and why. See the blog post, it can be just one line, or more paragraphs if needed.

@moorepants
Collaborator
@angadhn

@moorepants Yeah. I'm redoing the rolling disc.

@gilbertgede

I agree with @moorepants , we should add at least 1 test which is more complex, but I am not sure how complex we should make it. @angadhn is working on a test right now though, so something should be up soon.

@certik
Owner

Excellent. Please ping me when it's ready for review.

@angadhn

Everyone, I just wanted to give a heads up that I will be rebasing this branch; I need to have access to the energy functions that just got merged into master. as I will need to use this for the 'LagrangesMethod' tests. And I won't amend the commits this time around ;)

angadhn added some commits
@angadhn angadhn Revised6: LagrangesMethod class for sympy.physics.mechanics. a054386
@angadhn angadhn Minor edits. 9f5999f
@angadhn angadhn ADDED ROLLING DISC TEST AND MINOR EDITS.
The rolling disc test has been included.

Two minor edits have been included:
- addition of the energy function to determine the Lagrangian.
- slight modification of the 'rhs' method in the LagrangesMethod class.
1171a1c
@angadhn

@certik I have added the test for the rolling disc.

@hazelnusse @moorepants I tried to run the rolling disc with the non-minimal coordinates and so did @gilbertgede. It took forever to run. 12 mins into running the test, there was still no output so we decided to stick with the three generalized coordinates. I guess having trivial generalized (speeds as defined by Lagrange) slows things down drastically when we have 6 generalized coordinates AND three constraint equations.

@travisbot

This pull request passes (merged 1171a1c into 2797590).

sympy/physics/mechanics/tests/test_lagrange.py
((178 lines not shown))
+ C.set_vel(N, 0)
+ Dmc = C.locatenew('Dmc', r * L.z)
+ Dmc.v2pt_theory(C, N, R)
+
+ # Forming the inertia dyadic.
+ I = inertia(L, m / 4 * r**2, m / 2 * r**2, m / 4 * r**2)
+ BodyD = RigidBody('BodyD', Dmc, R, m, (I, Dmc))
+
+ # Finally we form the equations of motion, using the same steps we did
+ # before. Supply the Lagrangian, the generalized speeds.
+ T = kinetic_energy(N, BodyD)
+ BodyD.set_potential_energy(- m * g * r * cos(q2))
+ V = potential_energy(BodyD)
+ Lag = T - V
+ q = [q1, q2, q3]
+ l = LagrangesMethod(Lag, q)
@moorepants Collaborator

I'm wondering if this would flow better like:

BodyD.set_potential_energy(- m * g * r * cos(q2))
q = [q1, q2, q3]
l = LagrangesMethod(BodyD, (q1, q2, q3))

For 99% of problems the Lagrangian is T - V. I'm not sure why it would be otherwise. So why force all this extra typing? The above will mean that the kinetic energy and potential energy calls happen inside the Lagrange class.

So for that 1% where you might want to define a different Lagrangian, why not just have either an optional keyword argument or check the first object to see if it is a Body, Particle, or an expression.

So you can either type:

LagrangesMethod(body1, body2, particle1, particle2, (q1, q2, q3))

or

LagrangesMethod(L, (q1, q2, q3))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@moorepants
Collaborator

I fine with the tests as is. Nice job on this!

I've only the one comment above that I'm curious about the design of the code and why the Lagrangian can't be formed in the class.

@angadhn

How about adding a 'Lagrangian' function in the functions.py instead. That way users can feed in a Lagrangian however they desire.
I'm thinking of something like this in 'functions.py'-

def Lagrangian(frame, *body):
    "compute Lagrangian and return it here"

The frame would still be needed to determine the kinetic energy. What do you think?

@angadhn

Also a git question-
After having rebased and committed once (without amending), will it be okay to 'git commit --amend'?

@Krastanov
Collaborator

SymPy Bot Summary: :red_circle: There were test failures.

@angadhn: Please fix the test failures.

Test command: setup.py test
master hash: 2797590
branch hash: 1171a1c

Interpreter 1: :red_circle: There were test failures.

Interpreter: /usr/local/bin/python2.5 (2.5.6-final-0)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYjvUiDA

Interpreter 2: :red_circle: There were test failures.

Interpreter: /usr/bin/python2.7 (2.7.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY1Y8iDA

Interpreter 3: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python3.2 (3.2.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY2bYiDA

Build HTML Docs: :eight_spoked_asterisk: All tests have passed.

Docs build command: make html-errors
Sphinx version: 1.1.3

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYr_AhDA

Automatic review by SymPy Bot.

@hazelnusse

@angadhn @gilbertgede @moorepants @certik

I'm concerned about the test taking 12 minutes. Where was the majority of the time spent? Changing back to the minimal coordinates and speeds seems to imply you won't be testing functionality that is really important to this class. Maybe I'm missing something but to me, 12 minutes seems to be an indication that your algorithm needs improvement and/or the way you are using core sympy functionality is less than optimal. Can you explain why it was taking so long?

@asmeurer
Owner

git commit --amend is essentially the same as git rebase HEAD~1, as far as changing the history goes.

@angadhn

@hazelnusse Things get slow at the mass matrix inversion. The 'form_lagranges_equations' method doesn't stall so I'm wondering if maybe the issue lies in inverting a huge mass matrix, in this case, it was a 15 by 15 matrix.

@smichr
Collaborator

If you use git commit -a (a=all, not ammend) it will do what you are thinking, I believe: just add the changed files and put you in the editor to add a message. If you do git commit -a -m "a title" for a title-only commit and it will bypass the editing step.

@hazelnusse

@angadhn
Is the mass matrix inversion for solving for the q double dot's and the lagrange multipliers? Is this done by the class, or in user code after the mass matrix and forcing function has been obtained? If the former, think this is a pretty clear example of why coupling your class with a linear system solver is dicey. Also, what technique were you using to do the matrix inverse? I think in my experience, doing something like M.adjugate() / M.det(method='...') was pretty fast, even for relatively big matrices.

@angadhn

@hazelnusse yes. and also solves for the qdots (which is of course trivial.) .
And yes, this is done within clsas in the 'rhs' method. The inversion of the mass matrix MM is done by

MM.inv("GE", try_block_diag = True)

I will try it out the way you have suggested and see how it works.

@angadhn

@hazelnusse no dice. it seems to be taking forever to compute the determinant as well as the adjoint of the mass matrix.

@angadhn

I stand corrected. so the issue seems to be with printing the adjugate and not calculating it. As for just using the the method buil in to the class, that does take forever.

@hazelnusse
@angadhn

the method == rhs method to solve for the states.

@angadhn angadhn ADDED LAGRANGIAN AND SOME OTHER CORRECTIONS
A Lagrangian function has been added to 'functions.py' and appropriately
used now in the tests.

Corrected an error in the docstring of 'kinetic_energy'.

Modified the 'potential_energy' methods for both Particle and RigidBody;
the default value is now zero thus preventing the user from having to
set potential energy when there is none.
ba38ebe
@angadhn

@smichr @asmeurer Thanks for the information. Very helpful.

@hazelnusse Thanks for showing the path to faster inversions. Very insightful. Nonetheless, I have spent a good few hours since we met trying to decipher the gargantuan 'rhs' terms and I didn't find a feasible way to make those terms any simpler. It's quite horrendous as you may have seen.

@hazelnusse @moorepants @gilbertgede I will try to find a "simpler" test case for nonholnomicity later if you are okay with it as I personally feel this doesn't fit the bill in the case of using Lagrange's method.

@travisbot

This pull request passes (merged ba38ebe into 2797590).

@Krastanov
Collaborator

SymPy Bot Summary: :eight_spoked_asterisk: All tests have passed.

Test command: setup.py test
master hash: 3b77c94
branch hash: ba38ebe

Interpreter 1: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/local/bin/python2.5 (2.5.6-final-0)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY148iDA

Interpreter 2: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python2.7 (2.7.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY2rYiDA

Interpreter 3: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python3.2 (3.2.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYiYAiDA

Build HTML Docs: :eight_spoked_asterisk: All tests have passed.

Docs build command: make html-errors
Sphinx version: 1.1.3

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY7qMjDA

Automatic review by SymPy Bot.

@certik
Owner

It's definitely ok with me to merge this now, if you want. Or should I wait?

@angadhn angadhn ADDED WHITESPACE AFTER A COMMA.
I had forgotten to address one of Ondrej's comments on PR 1407 so I have
taken care of that here.
ce17e7b
@travisbot

This pull request passes (merged ce17e7b into 2797590).

@Krastanov
Collaborator

SymPy Bot Summary: :eight_spoked_asterisk: All tests have passed.

Test command: setup.py test
master hash: b0a9866
branch hash: ce17e7b

Interpreter 1: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/local/bin/python2.5 (2.5.6-final-0)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYvowjDA

Interpreter 2: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python2.7 (2.7.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYocYiDA

Interpreter 3: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python3.2 (3.2.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY8aMjDA

Build HTML Docs: :eight_spoked_asterisk: All tests have passed.

Docs build command: make html-errors
Sphinx version: 1.1.3

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY0t0iDA

Automatic review by SymPy Bot.

@gilbertgede

@certik I'm OK with this going in - the tests aren't as complex as they could be, but I believe that all the functionality is covered.

@angadhn

@certik With that last commit I'm OK with this going in too.

@travisbot

This pull request passes (merged d04c599 into 2797590).

@Krastanov
Collaborator

SymPy Bot Summary: :eight_spoked_asterisk: All tests have passed.

Test command: setup.py test
master hash: b0a9866
branch hash: d04c599

Interpreter 1: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/local/bin/python2.5 (2.5.6-final-0)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY_dUiDA

Interpreter 2: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python2.7 (2.7.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYwYwjDA

Interpreter 3: :eight_spoked_asterisk: All tests have passed.

Interpreter: /usr/bin/python3.2 (3.2.3-candidate-2)
Architecture: Linux (64-bit)
Cache: yes

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY_NUiDA

Build HTML Docs: :eight_spoked_asterisk: All tests have passed.

Docs build command: make html-errors
Sphinx version: 1.1.3

Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYrbMjDA

Automatic review by SymPy Bot.

@certik certik merged commit 86696a5 into sympy:master

1 check passed

Details default The Travis build passed
@certik
Owner

Excellent, it's in!

@angadhn

Everyone- thanks a lot for all the input.

@asmeurer
Owner

This has led to some test failures in master. See #1370 (comment).

@asmeurer
Owner

Actually, I can't reproduce those. They might be related to the copy_py3k_sympy feature of sympy-bot.

@gilbertgede

Yeah, it appears to only be with the python3 tests.

Is this just a random issue, or will this code keep popping up as an issue?

@asmeurer
Owner

Hopefully it will go away. I think it was an issue of py3k-sympy not being updated correctly.

@certik
Owner

The Python 3.2 tests are well tested by travis:

http://travis-ci.org/#!/sympy/sympy/builds

you can see that there are no failures in master, except some bugs in Travis itself.

@angadhn

That's a relief! Thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Commits on Aug 10, 2012
  1. @angadhn
  2. @angadhn

    Minor edits.

    angadhn authored
  3. @angadhn

    ADDED ROLLING DISC TEST AND MINOR EDITS.

    angadhn authored
    The rolling disc test has been included.
    
    Two minor edits have been included:
    - addition of the energy function to determine the Lagrangian.
    - slight modification of the 'rhs' method in the LagrangesMethod class.
Commits on Aug 11, 2012
  1. @angadhn

    ADDED LAGRANGIAN AND SOME OTHER CORRECTIONS

    angadhn authored
    A Lagrangian function has been added to 'functions.py' and appropriately
    used now in the tests.
    
    Corrected an error in the docstring of 'kinetic_energy'.
    
    Modified the 'potential_energy' methods for both Particle and RigidBody;
    the default value is now zero thus preventing the user from having to
    set potential energy when there is none.
  2. @angadhn

    ADDED WHITESPACE AFTER A COMMA.

    angadhn authored
    I had forgotten to address one of Ondrej's comments on PR 1407 so I have
    taken care of that here.
Commits on Aug 12, 2012
  1. @angadhn
This page is out of date. Refresh to see the latest.
View
4 doc/src/modules/physics/mechanics/examples.rst
@@ -68,7 +68,7 @@ equations in the end. ::
Kinematic differential equations; how the generalized coordinate time
derivatives relate to generalized speeds. Here these were computed by hand. ::
- >>> kd = [q1d - u3/cos(q3), q2d - u1, q3d - u2 + u3 * tan(q2)]
+ >>> kd = [q1d - u3/cos(q2), q2d - u1, q3d - u2 + u3 * tan(q2)]
Creation of the force list; it is the gravitational force at the center of mass of
the disc. Then we create the disc by assigning a Point to the center of mass
@@ -145,7 +145,7 @@ and along the ground in an perpendicular direction. ::
>>> vel = Dmc.v2pt_theory(C, N, R)
>>> acc = Dmc.a2pt_theory(C, N, R)
>>> I = inertia(L, m / 4 * r**2, m / 2 * r**2, m / 4 * r**2)
- >>> kd = [q1d - u3/cos(q3), q2d - u1, q3d - u2 + u3 * tan(q2)]
+ >>> kd = [q1d - u3/cos(q2), q2d - u1, q3d - u2 + u3 * tan(q2)]
Just as we previously introduced three speeds as part of this process, we also
introduce three forces; they are in the same direction as the speeds, and
View
4 sympy/physics/mechanics/__init__.py
@@ -34,3 +34,7 @@
import essential
from essential import *
__all__.extend(essential.__all__)
+
+import lagrange
+from lagrange import *
+__all__.extend(lagrange.__all__)
View
63 sympy/physics/mechanics/functions.py
@@ -13,7 +13,8 @@
'linear_momentum',
'angular_momentum',
'kinetic_energy',
- 'potential_energy']
+ 'potential_energy',
+ 'Lagrangian']
from sympy.physics.mechanics.essential import (Vector, Dyadic, ReferenceFrame,
MechanicsStrPrinter,
@@ -577,7 +578,8 @@ def kinetic_energy(frame, *body):
==========
frame : ReferenceFrame
- The frame in which angular momentum is desired.
+ The frame in which the velocity or angular velocity of the body is
+ defined.
body1, body2, body3... : Particle and/or RigidBody
The body (or bodies) whose kinetic energy is required.
@@ -650,7 +652,6 @@ def potential_energy(*body):
>>> a = ReferenceFrame('a')
>>> I = outer(N.z, N.z)
>>> A = RigidBody('A', Ac, a, M, (I, Ac))
- >>> BL = [Pa, A]
>>> Pa.set_potential_energy(m * g * h)
>>> A.set_potential_energy(M * g * h)
>>> potential_energy(Pa, A)
@@ -665,3 +666,59 @@ def potential_energy(*body):
else:
raise TypeError('*body must have only Particle or RigidBody')
return pe_sys
+
+def Lagrangian(frame, *body):
+ """Lagrangian of a multibody system.
+
+ This function returns the Lagrangian of a system of Particle's and/or
+ RigidBody's. The Lagrangian of such a system is equal to the difference
+ between the kinetic energies and potential energies of its constituents. If
+ T and V are the kinetic and potential energies of a system then it's
+ Lagrangian, L, is defined as
+
+ L = T - V
+
+ The Lagrangian is a scalar.
+
+ Parameters
+ ==========
+
+ frame : ReferenceFrame
+ The frame in which the velocity or angular velocity of the body is
+ defined to determine the kinetic energy.
+
+ body1, body2, body3... : Particle and/or RigidBody
+ The body (or bodies) whose kinetic energy is required.
+
+ Examples
+ ========
+
+ >>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame
+ >>> from sympy.physics.mechanics import RigidBody, outer, Lagrangian
+ >>> from sympy import symbols
+ >>> M, m, g, h = symbols('M m g h')
+ >>> N = ReferenceFrame('N')
+ >>> O = Point('O')
+ >>> O.set_vel(N, 0 * N.x)
+ >>> P = O.locatenew('P', 1 * N.x)
+ >>> P.set_vel(N, 10 * N.x)
+ >>> Pa = Particle('Pa', P, 1)
+ >>> Ac = O.locatenew('Ac', 2 * N.y)
+ >>> Ac.set_vel(N, 5 * N.y)
+ >>> a = ReferenceFrame('a')
+ >>> a.set_ang_vel(N, 10 * N.z)
+ >>> I = outer(N.z, N.z)
+ >>> A = RigidBody('A', Ac, a, 20, (I, Ac))
+ >>> Pa.set_potential_energy(m * g * h)
+ >>> A.set_potential_energy(M * g * h)
+ >>> Lagrangian(N, Pa, A)
+ -M*g*h - g*h*m + 350
+
+ """
+
+ if not isinstance(frame, ReferenceFrame):
+ raise TypeError('Please supply a valid ReferenceFrame')
+ for e in body:
+ if not isinstance(e, (RigidBody, Particle)):
+ raise TypeError('*body must have only Particle or RigidBody')
+ return kinetic_energy(frame, *body) - potential_energy(*body)
View
332 sympy/physics/mechanics/lagrange.py
@@ -0,0 +1,332 @@
+__all__ = ['LagrangesMethod']
+
+from sympy import diff, zeros, Matrix, eye, sympify
+from sympy.physics.mechanics import (dynamicsymbols, ReferenceFrame, Point)
+
+class LagrangesMethod(object):
@moorepants Collaborator

the Kane method class isn't called KanesMethod, so why do that here? It's just more to type. There are no other naming conflicts in this module are there?

@angadhn
angadhn added a note

So, @hazelnusse, @gilbertgede and I talked about this in the lab yesterday. We talked about how other methods of deriving EOMs are referred to i.e. 'Kane's method', 'Hamilton's method', etc and how maybe the "Kane" class could be renamed to "KanesMethod" possibly. Also, just calling it "Lagrange" might be ambiguous.

@moorepants Collaborator

Ok great, we should change the name of the Kane Class too then.

@moorepants I will take care of that in another PR that I am going to open soon.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+ """Lagrange's method object.
+
+ This object generates the equations of motion in a two step procedure. The
+ first step involves the initialization of LagrangesMethod by supplying the
+ Lagrangian and a list of the generalized coordinates, at the bare minimum.
+ If there are any constraint equations, they can be supplied as keyword
+ arguments. The Lagrangian multipliers are automatically generated and are
+ equal in number to the constraint equations.Similarly any non-conservative
+ forces can be supplied in a list (as described below and also shown in the
+ example) along with a ReferenceFrame. This is also discussed further in the
+ __init__ method.
+
+ Attributes
+ ==========
+
+ mass_matrix : Matrix
+ The system's mass matrix
+
+ forcing : Matrix
+ The system's forcing vector
+
+ mass_matrix_full : Matrix
+ The "mass matrix" for the qdot's, qdoubledot's, and the
@smichr Collaborator
smichr added a note

in my experience, mixing single and double quotes in a docstring generates sphinx errors...let's see if this happens with your request.

@angadhn
angadhn added a note

This hasn't been an issue in the "Kane" class, so I presumed it wouldn't be a problem here but that could've been an erroneous assumption.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+ lagrange multipliers (lam)
+
+ forcing_full : Matrix
+ The forcing vector for the qdot's, qdoubledot's and
+ lagrange multipliers (lam)
+
+ Examples
+ ========
+
+ This is a simple example for a one degree of freedom translational
+ spring-mass-damper.
+
+ In this example, we first need to do the kinematics.$
+ This involves creating generalized coordinates and its derivative.
+ Then we create a point and set its velocity in a frame::
+
+ >>> from sympy.physics.mechanics import LagrangesMethod, Lagrangian
+ >>> from sympy.physics.mechanics import ReferenceFrame, Particle, Point
+ >>> from sympy.physics.mechanics import dynamicsymbols, kinetic_energy
+ >>> from sympy import symbols
+ >>> q = dynamicsymbols('q')
+ >>> qd = dynamicsymbols('q', 1)
+ >>> m, k, b = symbols('m k b')
+ >>> N = ReferenceFrame('N')
+ >>> P = Point('P')
+ >>> P.set_vel(N, qd * N.x)
+
+ We need to then prepare the information as required by LagrangesMethod to
+ generate equations of motion.
+ First we create the Particle, which has a point attached to it.
+ Following this the lagrangian is created from the kinetic and potential
+ energies.
+ Then, a list of nonconservative forces/torques must be constructed, where
+ each entry in is a (Point, Vector) or (ReferenceFrame, Vector) tuple, where
+ the Vectors represent the nonconservative force or torque.
+
+ >>> Pa = Particle('Pa', P, m)
+ >>> Pa.set_potential_energy(k * q**2 / 2.0)
+ >>> L = Lagrangian(N, Pa)
+ >>> fl = [(P, -b * qd * N.x)]
+
+ Finally we can generate the equations of motion.
+ First we create the LagrangesMethod object.To do this one must supply an
+ the Lagrangian, the list of generalized coordinates. Also supplied are the
+ constraint equations, the forcelist and the inertial frame, if relevant.
+ Next we generate Lagrange's equations of motion, such that:
+ Lagrange's equations of motion = 0.
+ We have the equations of motion at this point.
+
+ >>> l = LagrangesMethod(L, [q], forcelist = fl, frame = N)
+ >>> print l.form_lagranges_equations()
+ [b*Derivative(q(t), t) + 1.0*k*q(t) + m*Derivative(q(t), t, t)]
+
+ We can also solve for the states using the 'rhs' method.
+
+ >>> print l.rhs()
+ [ Derivative(q(t), t)]
+ [(-b*Derivative(q(t), t) - 1.0*k*q(t))/m]
+
+ Please refer to the docstrings on each method for more details.
+
+ """
+
+ def __init__(self, Lagrangian, q_list, coneqs = None, forcelist = None, frame = None):
+ """Supply the following for the initialization of LagrangesMethod
+
+ Lagrangian : Sympifyable
+
+ q_list : list
+ A list of the generalized coordinates
+
+ coneqs : list
+ A list of the holonomic and non-holonomic constraint equations.
+ VERY IMPORTANT NOTE- The holonomic constraints must be
+ differentiated with respect to time and then included in coneqs.
+
+ forcelist : list
+ Takes a list of (Point, Vector) or (ReferenceFrame, Vector) tuples
+ which represent the force at a point or torque on a frame. This
+ feature is primarily to account for the nonconservative forces
+ amd/or moments.
+
+ frame : ReferenceFrame
+ Supply the inertial frame. This is used to determine the
+ generalized forces due to non-sonservative forces.
+
+ """
+
+ self._L = sympify(Lagrangian)
+ self.eom = None #initializing the eom Matrix
+ self._m_cd = Matrix([]) #Mass Matrix of differentiated coneqs
+ self._m_d = Matrix([]) #Mass Matrix of dynamic equations
+ self._f_cd = Matrix([]) #Forcing part of the diff coneqs
+ self._f_d = Matrix([]) #Forcing part of the dynamic equations
+ self.lam_coeffs = Matrix([]) #Initializing the coeffecients of lams
+
+ self.forcelist = forcelist
+ self.inertial = frame
+
+ self.lam_vec = Matrix([])
+
+ self._term1 = Matrix([])
+ self._term2 = Matrix([])
+ self._term3 = Matrix([])
+ self._term4 = Matrix([])
+
+ # Creating the qs, qdots and qdoubledots
+
+ q_list = list(q_list)
+ if not isinstance(q_list, list):
+ raise TypeError('Generalized coords. must be supplied in a list')
+ self._q = q_list
+ self._qdots = [diff(i, dynamicsymbols._t) for i in self._q]
+ self._qdoubledots = [diff(i, dynamicsymbols._t) for i in self._qdots]
+
+ self.coneqs = coneqs
+
+ def form_lagranges_equations(self):
+ """Method to form Lagrange's equations of motion.
+
+ Returns a vector of equations of motion using Lagrange's equations of
+ the second kind.
+
+ """
+
+ q = self._q
+ qd = self._qdots
+ qdd = self._qdoubledots
+ n = len(q)
+
+ #Putting the Lagrangian in a Matrix
+ L = Matrix([self._L])
+
+ #Determining the first term in Lagrange's EOM
+ self._term1 = L.jacobian(qd)
+ self._term1 = ((self._term1).diff(dynamicsymbols._t)).transpose()
+
+ #Determining the second term in Lagrange's EOM
+ self._term2 = (L.jacobian(q)).transpose()
+
+ #term1 and term2 will be there no matter what so leave them as they are
+
+ if self.coneqs != None:
+ coneqs = self.coneqs
+ #If there are coneqs supplied then the following will be created
+ coneqs = list(coneqs)
+ if not isinstance(coneqs, list):
+ raise TypeError('Enter the constraint equations in a list')
+
+ o = len(coneqs)
+
+ #Creating the multipliers
+ self.lam_vec = Matrix(dynamicsymbols('lam1:' + str(o + 1)))
+
+ #Extracting the coeffecients of the multipliers
+ coneqs_mat = Matrix(coneqs)
+ qd = self._qdots
+ self.lam_coeffs = -coneqs_mat.jacobian(qd)
+
+ #Determining the third term in Lagrange's EOM
+ #term3 = ((self.lam_vec).transpose() * self.lam_coeffs).transpose()
+ self._term3 = self.lam_coeffs.transpose() * self.lam_vec
+
+ #Taking the time derivative of the constraint equations
+ diffconeqs = [diff(i, dynamicsymbols._t) for i in coneqs]
+
+ #Extracting the coeffecients of the qdds from the diff coneqs
+ diffconeqs_mat = Matrix(diffconeqs)
+ qdd = self._qdoubledots
+ self._m_cd = diffconeqs_mat.jacobian(qdd)
+
+ #The remaining terms i.e. the 'forcing' terms in diff coneqs
+ qddzero = dict(zip(qdd, [0] * n))
+ self._f_cd = -diffconeqs_mat.subs(qddzero)
+
+ else:
+ self._term3 = zeros(n, 1)
+
+ if self.forcelist != None:
+ forcelist = self.forcelist
+ N = self.inertial
+ if not isinstance(N, ReferenceFrame):
+ raise TypeError('Enter a valid ReferenceFrame')
+ if not isinstance(forcelist, (list, tuple)):
+ raise TypeError('Forces must be supplied in a list of: lists'
+ ' or tuples')
+ self._term4 = zeros(n, 1)
+ for i,v in enumerate(qd):
+ for j,w in enumerate(forcelist):
+ if isinstance(w[0], ReferenceFrame):
+ speed = w[0].ang_vel_in(N)
+ self._term4[i] += speed.diff(v, N) & w[1]
+ if isinstance(w[0], Point):
+ speed = w[0].vel(N)
+ self._term4[i] += speed.diff(v, N) & w[1]
+ else:
+ raise TypeError('First entry in force pair is a point'
+ ' or frame')
+
+ else:
+ self._term4 = zeros(n, 1)
+
+ self.eom = self._term1 - self._term2 - self._term3 - self._term4
+
+ return self.eom
+
+ @property
+ def mass_matrix(self):
+ """ Returns the mass matrix, which is augmented by the Lagrange
+ multipliers, if necessary.
+
+ If the system is described by 'n' generalized coordinates and there are
+ no constraint equations then an n X n matrix is returned.
+
+ If there are 'n' generalized coordinates and 'm' constraint equations
+ have been supplied during initialization then an n X (n+m) matrix is
+ returned. The (n + m - 1)th and (n + m)th columns contain the
+ coefficients of the Lagrange multipliers.
+
+ """
+
+ if self.eom == None:
+ raise ValueError('Need to compute the equations of motion first')
+
+ #The 'dynamic' mass matrix is generated by the following
+ self._m_d = (self.eom).jacobian(self._qdoubledots)
+
+ if len(self.lam_coeffs) != 0:
+ return (self._m_d).row_join((self.lam_coeffs).transpose())
+ else:
+ return self._m_d
+
+ @property
+ def mass_matrix_full(self):
+ """ Augments the coefficients of qdots to the mass_matrix. """
+
+ n = len(self._q)
+ if self.eom == None:
+ raise ValueError('Need to compute the equations of motion first')
+ #THE FIRST TWO ROWS OF THE MATRIX
+ row1 = eye(n).row_join(zeros(n,n))
+ row2 = zeros(n,n).row_join(self.mass_matrix)
+ if self.coneqs != None:
+ m = len(self.coneqs)
+ I = eye(n).row_join(zeros(n,n+m))
+ below_eye = zeros(n+m,n)
+ A = (self.mass_matrix).col_join((self._m_cd).row_join(zeros(m,m)))
+ below_I = below_eye.row_join(A)
+ return I.col_join(below_I)
+ else:
+ A = row1.col_join(row2)
+ return A
+
+ @property
+ def forcing(self):
@smichr Collaborator
smichr added a note

docstring?

@angadhn Check the formatting for single line docstrings that shows up elsewhere and try to match that:
"""Some descriptive sentence."""

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+ """ Returns the forcing vector from 'lagranges_equations' method. """
+
+ if self.eom == None:
+ raise ValueError('Need to compute the equations of motion first')
+
+ qdd = self._qdoubledots
+ qddzero = dict(zip(qdd, [0] * len(qdd)))
@smichr Collaborator
smichr added a note

PEP8, I believe, would suggest removing spaces around *

@moorepants Collaborator

Are you sure?

http://www.python.org/dev/peps/pep-0008/#other-recommendations

Seems like spaces around these operators is preferred.

@certik Owner
certik added a note

Yes, it seems to me as well, that PEP8 says you should have spaces around "*". In either case, I would use whatever looks more readable on the case by case basis, in this PR it seems that having spaces is more readable.

@asmeurer Owner
asmeurer added a note

We usually differ from PEP 8 in SymPy for * and ** because it makes reading expressions like 3*x**2 + 2*x**4 easier. Here I would say it's justified, though really either way looks fine.

@smichr Collaborator
smichr added a note

I reread PEP8 and see that I was wrong about the * -- I was running under the sympy-modifed version (that Aaron cites). So although we are doing exactly what PEP8 says not to do I think it looks better, too.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+
+ if self.coneqs != None:
+ lam = self.lam_vec
+ lamzero = dict(zip(lam, [0] * len(lam)))
+
+ #The forcing terms from the eoms
+ self._f_d = -((self.eom).subs(qddzero)).subs(lamzero)
+
+ else:
+ #The forcing terms from the eoms
+ self._f_d = -(self.eom).subs(qddzero)
+
+ return self._f_d
+
+ @property
+ def forcing_full(self):
@smichr Collaborator
smichr added a note

docstring?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+ """ Augments qdots to the forcing vector above. """
+
+ if self.eom == None:
+ raise ValueError('Need to compute the equations of motion first')
+ if self.coneqs != None:
+ return (Matrix(self._qdots)).col_join((self.forcing).col_join(self._f_cd))
+ else:
+ return (Matrix(self._qdots)).col_join(self.forcing)
+
+ def rhs(self, method = "GE"):
+ """ Returns equations that can be solved numerically
+
+ Parameters
+ ==========
+
+ method : string
+ The method by which matrix inversion of mass_matrix_full must be
+ performed such as Gauss Elimination or LU decomposition.
+
+ """
+
+ # TODO- should probably use the matinvmul method from Kane
+
+ return ((self.mass_matrix_full).inv(method, try_block_diag = True) *
+ self.forcing_full)
View
9 sympy/physics/mechanics/particle.py
@@ -41,7 +41,7 @@ def __init__(self, name, point, mass):
self._name = name
self.set_mass(mass)
self.set_point(point)
- self._pe = None
+ self._pe = sympify(0)
def __str__(self):
return self._name
@@ -144,7 +144,7 @@ def angular_momentum(self, point, frame):
def kinetic_energy(self, frame):
"""Kinetic energy of the particle
- The kinetic energy, T, of a particle,P, is given by
+ The kinetic energy, T, of a particle, P, is given by
'T = 1/2 m v^2'
@@ -218,7 +218,4 @@ def potential_energy(self):
"""
- if callable(self._pe) == True:
- return self._pe
- else:
- raise ValueError('Please set the potential energy of the Particle')
+ return self._pe
View
6 sympy/physics/mechanics/rigidbody.py
@@ -53,6 +53,7 @@ def __init__(self, name, masscenter, frame, mass, inertia):
self.set_mass(mass)
self.set_frame(frame)
self.set_inertia(inertia)
+ self._pe = sympify(0)
def __str__(self):
return self._name
@@ -276,7 +277,4 @@ def potential_energy(self):
"""
- if callable(self._pe) == True:
- return self._pe
- else:
- raise ValueError('Please set the potential energy of the RigidBody')
+ return self._pe
View
185 sympy/physics/mechanics/tests/test_lagrange.py
@@ -0,0 +1,185 @@
+from sympy.physics.mechanics import (dynamicsymbols, ReferenceFrame, Point,
+ RigidBody, LagrangesMethod, Particle,
+ kinetic_energy, dynamicsymbols, inertia,
+ potential_energy, Lagrangian)
+from sympy import symbols, pi, sin, cos, simplify, expand
+
+def test_disc_on_an_incline_plane():
+ # Disc rolling on an inclined plane
+ # First the generalized coordinates are created. The mass center of the
+ # disc is located from top vertex of the inclined plane by the generalized
+ # coordinate 'y'. The orientation of the disc is defined by the angle
+ # 'theta'. The mass of the disc is 'm' and its radius is 'R'. The length of
+ # the inclined path is 'l', the angle of inclination is 'alpha'. 'g' is the
+ # gravitational constant.
+ y, theta = dynamicsymbols('y theta')
+ yd, thetad = dynamicsymbols('y theta', 1)
+ m, g, R, l, alpha = symbols('m g R l alpha')
+
+ # Next, we create the inertial reference frame 'N'. A reference frame 'A'
+ # is attached to the inclined plane. Finally a frame is created which is attached to the disk.
+ N = ReferenceFrame('N')
+ A = N.orientnew('A', 'Axis', [pi/2 -alpha, N.z])
+ B = A.orientnew('B', 'Axis', [-theta, A.z])
+
+ # Creating the disc 'D'; we create the point that represents the mass
+ # center of the disc and set its velocity. The inertia dyadic of the disc
+ # is created. Finally, we create the disc.
+ Do = Point('Do')
+ Do.set_vel(N, yd * A.x)
+ I = m * R**2 / 2 * B.z | B.z
+ D = RigidBody('D', Do, B, m, (I, Do))
+
+ # To construct the Lagrangian, 'L', of the disc, we determine its kinetic
+ # and potential energies, T and U, respectively. L is defined as the
+ # difference between T and U.
+ D.set_potential_energy(m * g * (l - y) * sin(alpha))
+ L = Lagrangian(N, D)
+
+ # We then create the list of generalized coordinates and constraint
+ # equations. The constraint arises due to the disc rolling without slip on
+ # on the inclined path. Also, the constraint is holonomic but we supply the
+ # differentiated holonomic equation as the 'LagrangesMethod' class requires
+ # that. We then invoke the 'LagrangesMethod' class and supply it the
+ # necessary arguments and generate the equations of motion. The'rhs' method
+ # solves for the q_double_dots (i.e. the second derivative with respect to
+ # time of the generalized coordinates and the lagrange multiplers.
+ q = [y, theta]
+ coneq = [yd - R * thetad]
+ m = LagrangesMethod(L, q, coneq)
+ m.form_lagranges_equations()
+ assert m.rhs()[2] == 2*g*sin(alpha)/3
+
+def test_simp_pen():
+ # This tests that the equations generated by LagrangesMethod are identical
+ # to those obtained by hand calculations. The system under consideration is
+ # the simple pendulum.
+ # We begin by creating the generalized coordinates as per the requirements
+ # of LagrangesMethod. Also we created the associate symbols
+ # that characterize the system: 'm' is the mass of the bob, l is the length
+ # of the massless rigid rod connecting the bob to a point O fixed in the
+ # inertial frame.
+ q, u = dynamicsymbols('q u')
+ qd, ud = dynamicsymbols('q u ', 1)
+ l, m, g = symbols('l m g')
+
+ # We then create the inertial frame and a frame attached to the massless
+ # string following which we define the inertial angular velocity of the
+ # string.
+ N = ReferenceFrame('N')
+ A = N.orientnew('A', 'Axis', [q, N.z])
+ A.set_ang_vel(N, qd * N.z)
+
+ # Next, we create the point O and fix it in the inertial frame. We then
+ # locate the point P to which the bob is attached. Its corresponding
+ # velocity is then determined by the 'two point formula'.
+ O = Point('O')
+ O.set_vel(N, 0)
+ P = O.locatenew('P', l * A.x)
+ P.v2pt_theory(O, N, A)
+
+ # The 'Particle' which represents the bob is then created and its
+ # Lagrangian generated.
+ Pa = Particle('Pa', P, m)
+ Pa.set_potential_energy(- m * g * l * cos(q))
+ L = Lagrangian(N, Pa)
+
+ # The 'LagrangesMethod' class is invoked to obtain equations of motion.
+ lm = LagrangesMethod(L, [q])
+ lm.form_lagranges_equations()
+ RHS = lm.rhs()
+ assert RHS[1] == -g*sin(q)/l
+
+def test_dub_pen():
+
+ # The system considered is the double pendulum. Like in the
+ # test of the simple pendulum above, we begin by creating the generalized
+ # coordinates and the simple generalized speeds and accelerations which
+ # will be used later. Following this we create frames and points necessary
+ # for the kinematics. The procedure isn't explicitly explained as this is
+ # similar to the simple pendulum. Also this is documented on the pydy.org
+ # website.
+ q1, q2 = dynamicsymbols('q1 q2')
+ q1d, q2d = dynamicsymbols('q1 q2', 1)
+ q1dd, q2dd = dynamicsymbols('q1 q2', 2)
+ u1, u2 = dynamicsymbols('u1 u2')
+ u1d, u2d = dynamicsymbols('u1 u2', 1)
+ l, m, g = symbols('l m g')
+
+ N = ReferenceFrame('N')
+ A = N.orientnew('A', 'Axis', [q1, N.z])
+ B = N.orientnew('B', 'Axis', [q2, N.z])
+
+ A.set_ang_vel(N, q1d * A.z)
+ B.set_ang_vel(N, q2d * A.z)
+
+ O = Point('O')
+ P = O.locatenew('P', l * A.x)
+ R = P.locatenew('R', l * B.x)
+
+ O.set_vel(N, 0)
+ P.v2pt_theory(O, N, A)
+ R.v2pt_theory(P, N, B)
+
+ ParP = Particle('ParP', P, m)
+ ParR = Particle('ParR', R, m)
+
+ ParP.set_potential_energy(- m * g * l * cos(q1))
+ ParR.set_potential_energy(- m * g * l * cos(q1) - m * g * l * cos(q2))
+ L = Lagrangian(N, ParP, ParR)
+ lm = LagrangesMethod(L, [q1, q2])
+ lm.form_lagranges_equations()
+
+ assert expand(l*m*(2*g*sin(q1) + l*sin(q1)*sin(q2)*q2dd
+ + l*sin(q1)*cos(q2)*q2d**2 - l*sin(q2)*cos(q1)*q2d**2
+ + l*cos(q1)*cos(q2)*q2dd + 2*l*q1dd) - (simplify(lm.eom[0]))) == 0
+ assert expand((l*m*(g*sin(q2) + l*sin(q1)*sin(q2)*q1dd
+ - l*sin(q1)*cos(q2)*q1d**2 + l*sin(q2)*cos(q1)*q1d**2
+ + l*cos(q1)*cos(q2)*q1dd + l*q2dd)) - (simplify(lm.eom[1]))) == 0
+
+def test_rolling_disc():
+ # Rolling Disc Example
+ # Here the rolling disc is formed from the contact point up, removing the
+ # need to introduce generalized speeds. Only 3 configuration and 3
+ # speed variables are need to describe this system, along with the
+ # disc's mass and radius, and the local gravity.
+ q1, q2, q3 = dynamicsymbols('q1 q2 q3')
+ q1d, q2d, q3d = dynamicsymbols('q1 q2 q3', 1)
+ r, m, g = symbols('r m g')
+
+ # The kinematics are formed by a series of simple rotations. Each simple
+ # rotation creates a new frame, and the next rotation is defined by the new
+ # frame's basis vectors. This example uses a 3-1-2 series of rotations, or
+ # Z, X, Y series of rotations. Angular velocity for this is defined using
+ # the second frame's basis (the lean frame).
+ N = ReferenceFrame('N')
+ Y = N.orientnew('Y', 'Axis', [q1, N.z])
+ L = Y.orientnew('L', 'Axis', [q2, Y.x])
+ R = L.orientnew('R', 'Axis', [q3, L.y])
+
+ # This is the translational kinematics. We create a point with no velocity
+ # in N; this is the contact point between the disc and ground. Next we form
+ # the position vector from the contact point to the disc's center of mass.
+ # Finally we form the velocity and acceleration of the disc.
+ C = Point('C')
+ C.set_vel(N, 0)
+ Dmc = C.locatenew('Dmc', r * L.z)
+ Dmc.v2pt_theory(C, N, R)
+
+ # Forming the inertia dyadic.
+ I = inertia(L, m / 4 * r**2, m / 2 * r**2, m / 4 * r**2)
+ BodyD = RigidBody('BodyD', Dmc, R, m, (I, Dmc))
+
+ # Finally we form the equations of motion, using the same steps we did
+ # before. Supply the Lagrangian, the generalized speeds.
+ BodyD.set_potential_energy(- m * g * r * cos(q2))
+ Lag = Lagrangian(N, BodyD)
+ q = [q1, q2, q3]
+ l = LagrangesMethod(Lag, q)
+ l.form_lagranges_equations()
+ RHS = l.rhs()
+ RHS.simplify()
+ assert (l.mass_matrix[3:6] == [0, 5*m*r**2/4, 0])
+ assert (RHS[4] == (-4*g*sin(q2) + 5*r*sin(q2)*cos(q2)*q1d**2
+ + 6*r*cos(q2)*q1d*q3d)/(5*r))
+ assert RHS[5] == (5*sin(q2)**2*q1d + 6*sin(q2)*q3d - q1d)*q2d/cos(q2)
Something went wrong with that request. Please try again.