Penultimate week of the GSoC period and it has been busy busy busy.

In my last post I had spoken about opening a PR for LagrangesMethod and about cleaning up the PR on energy functions. (links to both are in last week’s post and I won’t repeat them here.) Much of this week was spent cleaning up both those PRs and quite extensive testing on LagrangesMethod. The testing has been mostly successful. I shall explaing why ‘mostly’ in the a bit. The PR for the renergy functions has been merged and I’m just waiting for approval from ‘the boss’ so that LagrangesMethod can be merged too.

I would like to direct the reader to my proposal. In it I had said that I would write this class only for unconstrained systems. The idea was to modify this to be a ‘complete’ class post the GSoC period. But as we got into the week to begin working on this class, Gilbert and I decided that we would make it a full fledged Lagrange class; that could handle any kind of constraint on mechanical systems. Constraints on mechanical systems are basically of 2 types – configuration constraints (or holonomic constraints) and velocity constraints (or non-holonomic constraints). Depending on the methods used (Newton-Euler or Lagrange or Kane’s method and so on) these constraint equations are accounted for differently In the case of Lagrange’s method, there are additional terms due to these constraints that result in the introduction of the Lagrange multipliers. So, basically, repeating myself for the sake of clairty, one can now obtain the equations of motion in sympy.physics.mechanics using LagrangesMethod for any kind of system. I would even like to go out on a limb (quite literally under my current circumstances :P) and claim that it could be made use of for more generic applications involving ‘The Method of Lagrange Multipliers’ provided the user has the Lagrangian and constraint equations. (The documentation will however be limited to the domain of mechanical systems but shouldn’t be too hard to translate into something more generic for a user). The class also handles nonconservative forces thus making it a complete class.

In terms of testing each of these functionalities, I feel the tests are pretty thorough. I have tested the nonconservative force handling of the class on a simple spring-mass-damper system. I have tested the handling of the constraint equations using the famous ‘disc on an inclined plane’ test problem. Tests on more complex 3d systems have been performed like the rolling disc (more on this later). And then tests on a multibody system have been verified for a double pendulum. All of these work correctly; results have been compared to solutions by ahnd.

So with all of this down, why did I say it was “mostly successful”? Well, as it turns out, the tests work perfectly well when limited to problems involving planar motion. The results match up to those obtained by hand. But the results from the class get extremely nasty when dealing with more complex systems; I have implemented the rolling disc in two separate cases. In one test, I use the minimal set of generalized coordinates and the correct eoms are generated. But in another case I tried to use the non-minimal set of GCs and the equations generated are near impossible to comprehend (or I haven’t found the best way to deal with them yet). A big contribution of this messiness is due to the way in which Lagranges approach requires the definition of generalized speeds. In his approach, it is erquired for the generlized speeds to be ‘simple’ i.e. the gen. speeds are just derivatives of the gen coords. This is different in Kane’s approach where the generalized speeds can but needn’t necessarily be ‘simple’. From my experience, Kane’s generalized speeds are defined in a manner which make physical sense. This definitely validates why most dynamicists today (or so I have heard) prefer to choose Kane’s method on complex multi-body systems. The only way I can think of circumventing this situation in teh ‘LagrangesMethod’ class right now is using the minimal set of GCs for well known systems like the rolling disc and hope for the best.

Having all the additional functionality in this class and also playing with the rolling disc in particular has definitely led to a lot of insight but also taken a good chunk of time away from a period I wanted to dedicate to the ‘code output section’ which I have been unable to get started on. It looks like I will be unable to meet that one goal by the ‘hard’ pencils down date as I complete and fine tune the documentation (pending final approval of the Lagrange PR, of course). But I do feel that the time spent on Lagrange has been for the good. The code, I personally feel, is easy to read and appears to be easy to use. With people’s comments I was able to weed out all the unnecessary stuff. It is also ‘complete’ like I previously highlighted. I will continue to work on ‘code output’ post the GSoC period though as it’s usefulness is undeniable and also because of a development of a general sense of interest in coding (surprise surprise!).

Anyhow, apart from this, the other stuff I got done this week- I wrote up minor functions to compute a Lagrangian, changed how the potential energy function behaves. That’s it for this week. See you next week, one last time.

So, in last week’s post I had said that the Lagrange class had neared completion. By mid-week I had it functional, so I opened a discussion on the mailing list asking for suggestions to improve the class. Several people in the group suggested that it’d be better to supply all the parameters on initialization. At first I was loath to make this change. This was for the following few reasons-

- My primary concern was that I felt it was taking away from the clarity in generating the equations of motion. I have spoken about that in some detail in the conversation on the group page about how I feel dynamicists would probably prefer typing an extra couple lines as long as the procedure to obtain the equations of motion feel similar to the work done by hand. I’m a a fan of the (possibly false) sense of security that comes from bookkeeping even whilst on a computer. Of course, I had to accede that it made no sense to have so many methods when everything could be passed at initialization. The argument that users would essentially have to memorize a sequence or procedure was a good one.
- The other thing was that I felt was that this was making the structure of the ‘Lagrange’ class tremendously different from the ‘Kane’ class. But after discussions with Luke, Jason and Gilbert it became pretty apparent that even they felt that it might not have been the best way to implement that class. At one point I too felt like it was little weird that I was creating so many little methods which could’ve been merged into the initialization of the class itself. But I was just concentrating on staying true to the current structure of things in the mechanics package.
- And the final reason for my apprehension was that I would have to revamp the whole method after having spent quite some time on it.

But it was very clear after a point that things had to change. So, I spent a good chunk of time making the required changes.

Also, we decided on, what I feel is, a more appropriate name for the class changing it from ‘Lagrange’ to ‘LagrangesMethod’. Most of the equation derivation techniques in dynamics have the ‘Method’ attached to the founder’s name. It made even more sense for this class because Lagrange’s contributions are numerous so just calling a class ‘Lagrange’ could lead to some ambiguity.

Rewriting the class also helped me hone my Python skills some more. I had come across the keyword arguments construct several times in my preparation for the summer of code but I was a little reluctant to use it. It was probably because it was something that felt so alien to me as I have never seen something like that in my fledgling programming career. But with great explanations on the groups message (linked above), things were made clearer as to how it should be done.

So having rewritten the class, I added the docstrings for it. I’m not too pleased with that part right now, but I’m confident it will get better with more input on the PR discussion.

So having opened that PR, I thought I would get back to working on the documentation as I had planned. But I ended up going off on a tangent with the discussion that was sparked on PR 1407 which is the stuff that I have added on the energy functions. I spent a good chunk of time going through that and almost completely changing the way the ‘kinetic_energy’ function works.

On that same PR, there was a discussion about how a more readable error should be generated if a user calls the ‘potential_energy’ property for either a Particle or RigidBody without having first set the potential energy using the ‘set_potential_energy’ method. What at the time seemed an innocuous thing to repair became a little interesting challenge for me. Without going into too many more details, I was pleased to have found a relatively simple fix with the ‘callable’ function in Python with the help of the online forum ‘Stack Overflow’.

Having handled most of the recommendations on PR1407, I decided to skip on the documentation for the time being and returned to “complete” what would be the most important part relatied to the “LagrangesMethod” class- the test! While writing the class, I had written a little dummytest to check for the little tihngs but I hadn’t subjected a real dynamical system to the ultimate test (pun intended). I decided to test the well known ‘disc rolling down an inclined plane’ problem. Not to generate any suspense, but I would like to point out that in my proposal I had said that I would only concentrate on unconstrained systems. But Gilbert and I spent a little more time to make the ‘LagrangesMethod’ class more useful and complete. The class should now be able to handle any system i.e. constrained or unconstrained. A lot of the credit goes to Gilbert for helping me through the numerous confusions I had with the implementation of the constrained systems. But back to the test. I picked that system because it has a configuration constraint and we handle configuration constraints a little unconventionally in this class. I was a little anxious about how the results for this would turn out but ti worked like a charm. With the one test that I have written, which I think is a pretty good system to test, it appears that the ‘LagrangesMethod’ class works like a charm.

Anyhow, it’s now time to get some shuteye and more importantly rest the leg as I have been a little cavalier with it in the last couple of days. Until next week.

In the week accompanying and the weeks following the evaluation, I have been working on the implementation of Lagrange’s method to the ‘mechanics’ armory. A lot can be written on the method but I will spare the details here for now. In terms of progress on the coding itself, I feel like I have made significant strides in completing it. I expect to open a pull request in the middle of the upcoming week on Lagrange’s method, once I have some more of the ‘kinks’ ironed out.

I’m simultaneously working on the documentation aspect mentioned on my proposal. Currently, I’m working on adding to the theoretical or ‘text-book’ documentation that mechanics has on sympy docs. Specifically, I’m working on the theoretical material of the helper functions that I have added. Ideally, I will be opening a PR for this at the same time as the PR for Lagrange but that may not happen.

I would also like to apologise for missing out on the blog post last week. There were certain extenuating circumstances with respect to my physical health that I had to deal with coupled with a failure of my internet connection (which my ISP finally resolved on Wednesday). With regards to my health- I have had a troublesome knee problem which needed surgery. I had put off the surgery until the end of August but unfortunately the situation took a turn for the worse in the week following the evaluation. Anyhow it ended with me requiring an ACL reconstruction and meniscus removal. Unexpected but I have tried hard to not let it hamper productivity as best as possible. Nonetheless I apologise for my tardiness.

In terms of the bigger picture with respect to my GSoC goals, it may look like I’m a little behind schedule but since I’m simultaneously working on two of my goals right now which will get knocked off, in the worst case scenario, in the next two weeks. I will still have a couple more weeks to work on my final goal which is the addition and improvement of the code output functionality of mechanics. At this time, it looks like I will meet my goals and in the worst case scenario, I will be mid way through my last goal.

Another week has gone by and another PR has been opened. Added the kinetic energy and potential energy methods to the Particle and RigidBody classes. Also added functions to determine the same for a system of bodies. Implementing the kinetic energy method was fairly straightforward. A slightly different approach was used when it came to potential energies. It was decided that the best thing to do would be to let the user specify the potential energy scalar. So there is a set_potential_energy method and a read-only property called potential_energy.

On a brighter note, the two PRs that I had opened for the momenta functions and partial velocity functions have been merged. Chuffed about that.

As things stand now, we are set up beautifully to being the implementation of Lagrange’s method. The next week is primarily going to revolve around reading as much literature on the topic as possible because an attempt is going to be made to deal with non-holonomic systems too, unlike what was stated in my proposal i.e. only holonomic systems. I have personally never used Lagrange’s method on non-holonomic systems; instructors in the past have said it’s a messy affair and have stressed on the superiority of Kane’s method especially for this class of systems. I don’t doubt the superiority of Kane’s method but unfortunately, it’s not something that every dynamicist is aware of. So, it is definitely going to be a useful addition; a natural progression in the learning of dynamical methods starts at Newton’s method and is typically followed by Lagrange’s method and then to Kane’s method.

With the partial velocities helper functions and energy methods and functions, we are in a very good place in terms of implementing this method. Until next week then! (It’s 7am and I haven’t gotten any sleep last night. Over caffeinated. This post may not have made for pretty reading consequently)

Another week and more things learned. A new functionality was added this week to mechanics- the partial velocities function. This function basically spits out a list of lists of partial velocities.

A little background on the use of partial velocities in determining equations of motion- Partial velocities are used in Kane’s method to determine, what Kane calls, generalized active force and generalized inertia force. A partial velocity is defined as the partial derivative of a velocity (or angular velocity) vector with respect to an independent generalized speed, which is a scalar. The outcome of these partial velocities is thus a vector. As many of these can be computed as may be relevant to the system. Then these partial velocities of the points are dotted with the appropriate active and inertial forces at those points.

As easy and straightforward as that may sound in terms of implementation, I found it a bit of a tedious task to decide what was the best way for implementing it. It’s one thing to write code that is for stand alone purposes and where you have creative reign but it is another when it comes to writing code that must replace chunks of existing code.

At first, when I wrote out this function, it was just spitting out a list of partial velocities. That wasn’t too hard. But then Jason in the comment on last week’s post suggested that changes be made into the Kane class so that the partial velocities function be incorporated into it. So as it stood, my function was useless for anything else but looking at partial velocities. So back to the drawing board!

In the mechanics package, the scalar equations that are computed are stored in a matrix. So, I thought what would be ideal is to get a matrix of partial velocities that would mimic the partial velocities table that we get when we try to derive equations of motion by hand. So I went about doing that. Finally when I tested my code, I was getting an error message I couldn’t figure out. Now it was time to debug the code but with an error message that I couldn’t really grasp. I guess I went about trying to decipher the error in an unconventional way, where I went and rewrote each line of my code into the test file. But this was good because I finally figured out the issue with the function- The Matrix package in sympy doesn’t accept objects which are of the Vector class type that we typically use to represent vectors in mechanics.

Speaking about this with Gilbert finally led to the conclusion that nested lists would be best in this case. So after having written the same function several times, I was finally able to rewrite the function successfully and a pertinent PR opened. (I should add that in all this hoopla of writing and rewriting this function, I completely forgot about python list indexing beginning at 0. In hindsight it was hilarious to see myslef struggle with debugging my function one last time when I forgot one of Python’s quirks. I was up until the wee hours of the morning trying to figure that one out and then it hit me! Hah!). Currently I have been looking at the Kane.py module to see how to incorporate this function into it. In theory it shouldn’t be hard because the math that needs to be done fits right in, but it could be a bit more tedious than I recognize. But I’m confident that the function itself will not have to be tweaked to be used which is a good sign.

So this upcoming week, I will find a way to use this functionality in the Kane class. I will also be meeting with Gilbert and we will determine a plan in terms of how to implement a Lagrange class which is one of the major goals for this summer.

This has been a good week in terms of progress and clarity. In terms of completed code, I now have the entire suite of momenta functions working perfectly. The pull request that I linked to last week now has a series of commits. Several helpful comments were made in the PR; most of it had to with improving the readability of the docstrings and I’m extremely thankful for the time taken by everyone to review my work.

Onto the details of the commits now- I implemented methods to determine linear and angular momenta of both individual particles and rigidbodies. These are working well. As I said last week, I was having trouble figuring out how to deal with a system of bodies. Consequently, I was thrown off track for a bit where I tried to figure out how I would go about doing this without having to restructure the existing code. Sadly (when I think of the time I lost doing this but simultaneously glad that I figured it out eventually), the answer had been staring me in the face. Whenever anyone uses ‘mechanics’ to compute the equations of motion, they are required to set up a body list so all I had to do was tap into that. I ultimately wrote up a little procedure that computes both the momenta for a system of rigid bodies and/or particles. I then polished up the previous tests I had written and wrote some additional tests for the newer functionality. Overall, it was very satisfying to see the body of work that has been contributed, It will be a pretty useful tool when the time comes to implement the Newton-Euler method to determine the equations of motion, though I may not be the one to do it (at least not probably within this summer) as I plan to attack Lagrange’s method for equation derivation.

To that end, I spent some time this week trying to see how I would go about including Lagrange’s method. In my proposal, I have stated that the intention is to have Lagrange’s method for holonomic systems, but this is the more trivial case and is pretty well known to most dynamicists. I’m hoping that I can push the boundaries a little further and not be limited to just holonomic systems. The problem with the Lagrange approach is that it requires the determination of Lagrange multipliers whilst dealing with nonholonomic systems and this is not a method that is well documented in the classical textbooks. Finding material on this is a little tedious and I’m still uncertain about implementing this portion but I’m hopeful that it can be done.

So apart from the completion of code and a literature review on Lagrange, I have begun on a function to determine partial velocities and another for kinetic energies of particles and rigidbodies. These aren’t particularly hard but I’m constantly thinking of how I would go about using these functions when it comes to Lagrange’s equations of motion. (Another issue here is that there a couple ways to determine the equations of motion of holonomic systems even in Lagrange’s approach. One involves determining a Lagrangian and performing some operations on the Lagrangian and another approach involves determining just kinetic energies and generalized forces.) The coming week will lead to some clarity on that, it appears. I ill also be talking to my grad school advisor about this matter as I’m confident he will have more insight on this.

One thing that I have failed to talk about in my previous posts is the difficulty (read as ‘frustrations’) with using git. For someone with a pretty poor background in programming, the amount of jargon in python alone is hard to deal with. I remember when I was accustoming myself to git via the progit online book, it all seemed so easy. In the last few weeks, I have realized how much more harder it is to do the right things in the right order when you’re working on a real project and not on an example. Somwhow theory is always easier than practice. The frustrations that I had had lesser to with git and more to do with the fact that it was taking so long to master. But I took some solace after conversations with the other people in my lab who have worked on it when they said that it took them a while to get to grips with it. But nonetheless, I can see why it such an important tool. Over the last couple of weeks, I have gotten very familiar with not only the trivial aspects of making new branches and rebasing but also with the more powerful things like branching off of the right branches correctly and just in general better planning of projects. It was also also useful to learn that I could pull work from other people and collaborate on their work too.

One of the lessons that I have learned this week though is that I shouldn’t be worried about being a bit of a slow coder as everyone takes a while to get fast at things. What’s more important is to have ‘usable’ code. BUT I will be picking up the pace from here on, I believe. Python has definitely got me interested in programming like never before. All in all. a pretty good week, personally.

This has been a very interesting week for me. I have learned a lot more about the intricacies of Python and how sympy.physics.mechanics has been written. I should admit that as an engineer, my programming background isn’t as strong as that of others so maybe I’m not asking the right questions. But I’m also with several questions after writing my first ‘formal’ code.

I have opened a pull request with my functions for angular momenta for both particles and rigid body. I don’t intend to go into the details of the mechanics of it but one thing that I have asked myself upon writing my first set of code is whether it is useful in and of itself. I’m not convinced about it being the best implementation. The code definitely spits out the (linear or angular) momentum of any individual component of a system but I haven’t found a suitable way to determine the momentum of an entire system that may compromise a multi-body system i.e. one that has more than one kind of object, in this case ‘Particles’ and ‘RigidBodies’. I believe it is fairly trivial to deal with just one of these object types when they comprise a system but it’s a lot more intriguing to think about how to incorporate a system with multiple objects. The way I’m thinking of it right now is that ‘Particles’ and ‘RigidBodies’ must be derived classes of a more general class, somethig like ‘class System(object):’ with the rest of the code following that. But I also feel, at the same time, that I’m not approaching this problem the right way. I gave myself a couple days to figure this out and I was really hoping that I could figure this out on my own but I guess I will need more guidance on this issue. This has also made me question the simple code that I have to determine the energies, kinetic and potential, of either system. It is useful as it stands for individual components of multibody systems as it is trivial to have a method to determine either energy of a particle or a rigid body but adding this functionality keeping in mind that it has to be used to determine the Lagrangian of a more complex multibody system is also proving to be a bit of barrier. I’m uncertain but I feel that the answer to this also lies in what I have spoken about prior- the issues with angualr momentum of a system.

This has been primarily what has dominated my time since I have worked on my momenta functions.

In terms of tangibles, I have the momenta code up and running. I have also tested the better checking methods for vectors and dyadics as per PR 1269. (Note- That particular pull request has proven to be extremely crucial to my work. When angular momentum vector of a rigid body is to be found about the mass center, the generic formula for angular momentum that I implemented in ‘mechanics’ falls apart because in the math that follows, a cross product between a zero vector and a non zero vector needs to be performed. It turns out that the methods currently available in ‘mechanics’ aren’t up to task to cope with this (even though cross products between a vector and a zero vector are performed correctly.)) After having spent almost 2 days trying to examine why this generic formula wasn’t working, I ended up writing some patchy code with some conditional statements for the angular momentum to work with the current setup. After this I went through Gilbert’s updated code (PR 1269) and saw that it handles the math being performed in my angular momentum function perfectly. I hope that both that PR and mine get merged in that order. And as I am also waiting on that merge, I have been furhter examining the pydy.org website. I have jotted down several recommendations that I have to make for that page. This will prove to be important because that is the page that students in UC Davis use to learn sympy.physics,mechanics. As this has already turned into a bit of a prosaic blog post, I shall end this here. But, yes this has been a long week and a lot has been learned and I’m pretty satisfied with what I’ve learned and am even more invigorated for all that is in store in the coming weeks.