Vector Calculus module

[Prasoon Shukla]

Hello all,

I was going through the project ideas for GSoC 2013 and I saw an entry for implementation of vector calculus. As I have taken a semester's worth of vector calculus, this looked like something I could do.

I am not saying that I'll be going with this proposal a 100%, but I may.

Anyway, here are some of the questions that I have.

Q. Has it already been implemeneted to some extent?

I began with grep for keywords such as curl, divergence etc. Barely a few lines showed up. Then, I tried a google search and that led me to the mechanics module. Unfortunately, I am not familiar with the mechanics module. Therefore, I decided that it was best to ask the community about it. So, has vector calculus been implemented anywahere *inside* the mechanics module? That seems quite unlikely to me (because of the grep), still, I'm just checking whether anyone else has done some wrok on this that isn't on master.

Now, assuming that there is no implementation of this yet, I have the following questions:

Q. In case I decide to make my GSoC proposal for the addition of this new module, what functionality is expected from this new component other than the obvious?

Obviously, we would have to support general coordinate systems, tranformations, the usual operations for elements belonging to these vector fields and vector operations over scalar fields(curl, div on vectors, laplacian etc on scalars), integration over paths, surfaces and volumes, the usual theorems of vector calculus (Gauss, Stokes), support for vector functions etc.

There is no order/priority to this list. These are just the few things that came off the top of my mind. Now, other than this relatively obvious functionality, what are the other important things that we might have a relation with this module?

As an example, we can have the solution of the PDE for the magnetic vector potential. This would depend both on the new vector calculus module (to get the PDEs to solve) as well as on the PDE solvers (to get the solution of PDEs).

So, basically, I'm looking in both directions here. If the funcionality of some other part of SymPy increases because of this new module, then it's well and good. But, on the other hand, we are going to have dependencies on many other components as well (multiple integration is an example). So, if someone can point out some of these dependencies that might pose a problem later on, that would be helpful.

Also, as I asked before, what might be some of the non-obvious things that a user might require? I can't think of any such thing right now but it is certainly worth a look.

Q. Would this project be big enough for a summer?

My personal thought on this is yes. Though I have only vaguely considered the *how* yet - I am focussing more on the *what* for now -  still, I think that a well made VectorCalculus module may require quite a bit of work.

If this much work is enough for a summer, well and good. If not, then what more can be added to this project so as to make it big enough for a summer's worth of coding?

These are just some of questions that occur to me as of yet. Please state your opinions/suggestions/answers below. Once I am clear on these points, then, I can move ahead with figuring out a sketch for the implementation.

Thanks

Prasoon Shukla

[Jason Moore]

Prasoon,

In the mechanics module we've implemented the elements of vector calculus needed to do classical mechanics. This is a subset of what vector calculus is and it is also not necessarily written in a general form. We have a Vector class but it is limited to 3 dimensions and tied to the concept of a reference frame. If a vector calculus module existed we could potentially modify mechanics to use that functionality. So, understanding the mechanics module will give you a good example of one use case for a vector calculus module, but it may only give some ideas about how to implement a generic vector calculus module.

Jason
moorepants.info
+01 530-601-9791

Prasoon Shukla

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy?hl=en.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[Prasoon Shukla]

Thank you Jason, I'll look into the mechanics module and report back on this thread soon.

[Aaron Meurer]

Other places you might look for things that are already implemented
are the matrix and the differential geometry modules (including matrix
expressions).

There is also line_integrate, which uses Curve from the geometry
module. I believe that like the mechanics, the geometry module is
3-dimensional only.

By the way, it's not exactly vector calculus, though it's somewhat
related, but another great idea for a GSoC project is to expand the
matrix expressions module, especially with respect to derivatives.
See https://code.google.com/p/sympy/issues/detail?id=2759.

Aaron Meurer

[Prasoon Shukla]

@Aaron: Thanks for the inputs. I'll use this weekend to go through the modules you and Jason have mentioned and also brush up on the vector calculus. The matrix expressions module seems interesting at a first glance, too, but for the time being, I'll just focus on vector calculus. Let me experiment a bit and hopefully, I'll be able to decide on which way to go from here.

On a side note, I am using Calculus, Volume 2, by Tom Apostol (Wiley) as a reference here particularly because it also has elements of linear algebra prior to the vector calculus part.

[someone]

> I'll just focus on vector calculus. Let me experiment a bit and
> hopefully, I'll be able to decide on which way to go from here.

I's not vector calculus but rather abstract vector algebra
but you might want to look at:

"Rule-Based Simplification in Vector-Product Spaces"
by Songxin Liang and David J. Jeffrey.

I have an almost finished implementation/port for Axiom:

https://github.com/raoulb/fricas_code/tree/master/vecalg

[Stefan Krastanov]

I would like to elaborate on the remarks already made by Aaron and Jason.

There are two general approaches to this project:
1. taking care of the components of the vectors
2. working only on the algebraic properties of the vectors


Approach 1 is what is done obviously in the "matrix" module. Approach
2 is what is done in the "matrix expressions" module.

The quantum module, the diffgeom module and the mechanics module take
a mix of the two approaches, but in any case they seem to be closer to
approach 1. The new tensor canonicalization module that is still a
pull request (by Pernici) takes mostly approach 2.

I think that you can base your project on any of the following
modules: mechanics, diffgeom and quantum. However, given that there
are a lot of techniques that are specific to 3D euclidean space I
believe the best fit would be to abstract/extend the mechanics module.

I would love to see the vector manipulations in `mechanics` abstracted
outside of it (however the time-dependence stuff seems quite specific
to what they do, so it could be a problem). I think (the others may
disagree) that this would be the most valuable formulation of your
project.

If you follow this suggestion, it would still be very useful for you
to study the whole representation framework employed by the quantum
module. It has a very neat implementation of abstract hilbert spaces
that can be represented in any basis.

[Prasoon Shukla]

@Stefan: I have been looking at the implementation of vectors in the mechanics module and, just as you said, it would be better to extend the methods already used in the mechanics module. I am still in the process of looking at the code. I'll report back soon.

@Raoul: I have not yet taken a course in abstract algebra though I am currently enrolled in a course of linear algebra. I can read it up of course. Do you foresee any problems with not knowing abstract algebra just yet?

[someone]

> @Raoul: I have not yet taken a course in abstract algebra though I am
> currently enrolled in a course of linear algebra. I can read it up of
> course. Do you foresee any problems with not knowing abstract algebra
> just yet?

Oh, this one is not so much about "algebra" in the sense of either
polynomial or commutative algebra but the main point is that the
vectors are component free, this is what "algebra" refers to here.

The math behind this is fairly easy, mostly how to make the clever
use of about 10 axioms or "rules".

This is of course no full GSOC project. But thinking about it
might help you to make a good proposal for general vector calculus.

Even if this is beyond your scope and project but in the long term,
component-wise vectors calculus and component-free vector algebra
should go hand in hand. (The situation is similar to the one of
matrix calculus with explicit matrices/entries and matrix expressions
ignoring the entries at all.)

[Gilbert Gede]

This sounds like it could be a good project (or rather a decent part of a GSoC project). Vector/ReferenceFrame were indeed written for the needs of the mechanics module, and not for general SymPy use. This meant a slight focus on time (as Stefan mentioned) and a big focus on multiple reference frames (for multi-body dynamics).

That being said, one good thing is that when using the mechanics module the Vector/Dyadic classes are never directly instantiated by a user. So if they need to be "generalized", we can probably extend them within mechanics to get back what we need.

I think the harder issue is determining an interface for a general use Vector/ReferenceFrame. This is probably of some relevance to Sachin Joglekar and his electrodynamics code as well as the other users of the physics module(s). I think that before you make a GSoC project proposal for this, that you get feedback from all relevant parties, and see if a common vector class will even be feasible (considering everyone's different needs). Making something that is all things to all people is probably impossible, but if you can get 80% of what people need and everyone can just extend that core into what their corner of physics needs, that would be a good thing. But if everyone's needs are too different, you might want to explore other possibilities for a GSoC project.

-Gilbert

[Sachin Joglekar]

Prasoon's idea does seem great. And as Gilbert pointed out, a nice interface for Vector and ReferenceFrame functions would go a long way in helping me build the module that I would propose in my project. Moreover, in an ideal case vector calculus is something that can be in sympy.physics, or in a module of its own, as it has various applications of its own, mechanics being one of the most important. Whether its a complete GSoC proposal or not is something that I cannot decide, but the idea is great nonetheless.

[Prasoon Shukla]

@Gilbert: I have been thinking along the lines of the mechanics module and as Stefan mentioned, I plan (at least right now) to take this idea forward with component wise operations on vectors. For component independent vector algebra, any numerical operation on any vector will inevitably get reduced to operations on scalars (individual components). Therefore, I believe, that a large part of component dependent vector operations will form a subset to a more complete component independent vector algebra. So definitely, in the long term, we will be able to extend this vector module to play nice with component free vector algebra, just as Raoul mentioned.

I will, at this stage, like to point out that my knowledge of mathematics is limited to vector calculus in R3; at least for the time being. The course I have done for vector calculus was only for vectors on R3. I have a working knowledge of vector spaces, linear transformation, some matrix theory but nothing from abstract algebra. I have started reading tensor analysis to have some idea about implementing vector operations on Rn (for example, the curl is defined in terms of permutation tensor). Right now, I will have no problems implementing the vector module on R3 but obviously, it would be great if it can be extended to Rn and, indeed, that is what I am trying to accomplish.

Also, and just as you mentioned, Vector in mechanics can be extended because the user does not directly interact with vector objects. This means that, at least for mechanics, we won't have too many problems as far as integration with the new module is concerned. Of course, this problem remains for other modules (quantum, electricity-magnetism modules proposed by Sachin). For those, I will interact with people in the community for their requirements and hopefully, decide on an idea that keeps everyone happy, or at least, which keeps most of the people happy.

Now, let us talk a bit about the requirement of reference frames. The reference frames will, I think, be required only for R3 (as in case of mechanics). I wish to abstract the idea of a reference frame as implemented in mechanics, albeit in a bit different way. The discussion on world based coordinate system is what I wish to point out here in this respect. Currently, we need to have a RefrenceFrame for any vector. World coordinates will allow vectors that just work and I think that it is a very important feature to have from a user's point of view. 

The current implementation that I have in mind will allow the user to directly instantiate a Vector. As an example, say I want to just initialize a vector and find out its curl. So, the user should be able do something like (just an example of how things could function):

>>> from vector import *

>>> from sympy.abc import x, y, z

>>> coord = CoordSystem('rect')

>>> v = Vector((-y, x, 0), (x, y, z), coord=coord)

>>> div(v)

0

>>> curl(v)

2*ez

Taking a clue from diffgeom, we are defining a coordinate system first (with 'rect' providing the basic cartesian coordinates, the definitions of which are inbuilt. We can have similar functionality for other commonly used coordinate systems). Then, we instantiate a vector and give it symbols to use for coordinates. The dimension of the vector can be calculated internally using the second n-tuple. The basis for any Rn space will be just the standard basis. Then, we can have component-wise operations that produce the result; in this case divergence and curl of the given vector field.

Also, for the n-dimension case, we can have something like:

>>> v = Vector(dim=7, symbol=y, coord=coord)

Here, the coord object will carry the relations with one of the standard coordinate systems and that will be used to define all the operations of vector calculus.( The coord object can be generated by performing operations on CoordSystem object of some standard coordinate system perhaps?)

This will generate a vector (unbound?) where the coordinates wrt to the standard basis (e1, ... , e7) are just (y1, y2, y3 ... y7).

Anyway, my ideas at this stage are sketchy at best - without knowing the requirements of the community, this speculation on the API would be little more than fruitless.

As should be obvious, this example represents nothing concrete; just a glimpse of some things I'd like to see.

To conclude, my point is that we should have direct interaction of the user with vectors, and not via reference frames. Also, the vector operations should be defined so that everything is just intuitive.

Again, as I am mentioned above, I am currently reading a book on tensor analysis. While I read the subject matter, I will like to gather opinions of others (specially people related with modules that will use/already use some form of vector calculus) on how best to implement this so that it will be sustainable in the long run, would be easily extensible, would be consistent with rest of SymPy codebase, and, so that it has a very generalized API that any other module can easily use.

One more thing. If there seems to be any mathematical inconsistency with my arguments, you can attribute that to my lack of knowledge of the requisite subject matter; as I mentioned before. Please feel free to correct any such inconsistencies you may find.

I shall be back again with some more reading. Until then, if you gentlemen have any comments/suggestions, please post them here.

Thanks.

[Stefan Krastanov]

First of all, I think it would be nice to start writing a wiki page
with the proposal as the general idea seems clear enough.

Second (and most important in this message), I would vote for
completely forgetting about anything outside 3D. Implementations of
vector calculus in Rn should be done (and are done to some extend) in
the geometrical algebra and differential geometry modules (and the
tensor canonicalization pull request by Pernici). The whole idea of a
vector calculus module, in my opinion, is that when working in R3 you
do not need all these generalized notions of, e.g, curl or vector
product in terms of completely antisymmetric tensors.

Obviously, we would love to have `rewrite_as` methods to couple it to
the diffgeom/ga/tensor modules but that does not mean it is a good
idea to use, for instance, the diffgeom module as a basis for your
project. It can be done and for puritans it will indeed be a deep
example of code deduplication/reuse and object oriented blah blah.
However, the code will be harder to understand and the obvious notions
of R3 calculus will be hidden behind their generalizations.

And simply from a practical standpoint, the mechanics module is
already based around R2/R3. It is more important to abstract that away
than to generalize it to cases that will not be used.

The rest of my comments are below, annotating your message.

> @Gilbert: I have been thinking along the lines of the mechanics module and
> as Stefan mentioned, I plan (at least right now) to take this idea forward
> with component wise operations on vectors. For component independent vector
> algebra, any numerical operation on any vector will inevitably get reduced
> to operations on scalars (individual components). Therefore, I believe, that
> a large part of component dependent vector operations will form a subset to
> a more complete component independent vector algebra. So definitely, in the
> long term, we will be able to extend this vector module to play nice with
> component free vector algebra, just as Raoul mentioned.

It is a bit unclear what you mean by "component dependent". Here are
the two ways this can be understood:

1. Simply using matrices without "metadata". I consider this a very
bad idea as a matrix is just a representation of the geometrical
entity vector and all that can be done this way is already in the
matrix module.
2. Using some symbolic object containing
- reference to the "world" in which we work
- reference to the reference/basis in which we represent the vector
- the coordinates in this basis, which is a matrix

EDIT: I see that you actually use the second approach. Great!

> Also, and just as you mentioned, Vector in mechanics can be extended because
> the user does not directly interact with vector objects. This means that, at
> least for mechanics, we won't have too many problems as far as integration
> with the new module is concerned. Of course, this problem remains for other
> modules (quantum, electricity-magnetism modules proposed by Sachin). For
> those, I will interact with people in the community for their requirements
> and hopefully, decide on an idea that keeps everyone happy, or at least,
> which keeps most of the people happy.

I would say, leave quantum away. Hilbert spaces in quantum mechanics
are sufficiently different from R3 to necessitate their own
implementation. On the other hand, I strongly suggest that this
project you are preparing be the basis for the electrodynamics project
by Sachin. You should definitely work together on this.

> Now, let us talk a bit about the requirement of reference frames. The
> reference frames will, I think, be required only for R3 (as in case of
> mechanics). I wish to abstract the idea of a reference frame as implemented
> in mechanics, albeit in a bit different way. The discussion on world based
> coordinate system is what I wish to point out here in this respect.
> Currently, we need to have a RefrenceFrame for any vector. World coordinates
> will allow vectors that just work and I think that it is a very important
> feature to have from a user's point of view.

Reference Frame makes sense outside mechanics (and for any dimension,
but as I said, I would like to see this project restricted to 3D done
right).

You might be confusing the following notions (or maybe you do not, but
in any case let's just spell them out):

- vector space vs affine space (or flat manifold): Physics and vector
analysis is __not__ done on vector spaces. Just look at your example
below: you have components of the vector (the vector space part) and
the coordinates of the point where this vector is (the manifold part).
- the basis is what you use to express the components of the vector
(the vector space part).
- the coordinate system is what you use to express the point in space
you are using.
- you can very well have two different cartesian coordinate systems
(centered at the same point but rotated wrt to each other or instead
parallel but with different origins). In the same way you can have two
different spherical coordinate systems.
- you may have noticed that the basis you use for the vectors can
change from point to point. For instance, the usual basis you use when
dealing with spherical coordinates is radial/longitudinal/latitudinal
vectors. However they change from point to point. This defines a
_frame_: having a frame is having a basis at each point of your space.
- importantly, there is a canonical frame for each coordinate system
(the one you obtain by having unit vectors in the directions in which
coordinates change).
- A reference frame can be (depending on textbook): the coordinate
system, the frame, or both taken together. I do not know which is in
the mechanics module. Just stick to it, and refactor it out into your
own module.

I think it would be good to refactor the above list, fix the
terminology you will be using and add it in your wiki page
application.

> The current implementation that I have in mind will allow the user to
> directly instantiate a Vector. As an example, say I want to just initialize
> a vector and find out its curl. So, the user should be able do something
> like (just an example of how things could function):
>
>>>> from vector import *
>>>> from sympy.abc import x, y, z
>>>> coord = CoordSystem('rect')
>>>> v = Vector((-y, x, 0), (x, y, z), coord=coord)
>>>> div(v)
> 0
>>>> curl(v)
> 2*ez

This seems sensible. You may want to abstract (x, y, z) out into a
Point class, as in diffgeom, however this might be excessive,
especially given that we do not mix coordinate systems with frames of
other coordinate systems.

> Taking a clue from diffgeom, we are defining a coordinate system first (with
> 'rect' providing the basic cartesian coordinates, the definitions of which
> are inbuilt. We can have similar functionality for other commonly used
> coordinate systems). Then, we instantiate a vector and give it symbols to
> use for coordinates. The dimension of the vector can be calculated
> internally using the second n-tuple. The basis for any Rn space will be just
> the standard basis. Then, we can have component-wise operations that produce
> the result; in this case divergence and curl of the given vector field.
> Also, for the n-dimension case, we can have something like:
>>>> v = Vector(dim=7, symbol=y, coord=coord)

Mostly agree, but drop dim!=3. And I have a lot of technical comments
concerning the the `symbol` argument, but they can wait for later.

>
> Here, the coord object will carry the relations with one of the standard
> coordinate systems and that will be used to define all the operations of
> vector calculus.( The coord object can be generated by performing operations
> on CoordSystem object of some standard coordinate system perhaps?)
>
> This will generate a vector (unbound?) where the coordinates wrt to the
> standard basis (e1, ... , e7) are just (y1, y2, y3 ... y7).

Concerning the very important "unbound?" question. In your first
example `Vector` is actually a vector field, i.e. bound vector for
physicists (there is a difference between vector field and a bound
vector but drop it for the moment). In your second example `Vector` is
actually a vector, not a vector field. Vector analysis cares about
vector fields, so just completely drop the unbounded vectors from your
proposal. They have place in linear algebra module, in matrix
expressions module and in hilbert spaces module (the quantum module),
not here.



PS

By the way, this question show an issue that must be discussed a lot
before we are clear about how to proceed. The mechanics module deals
with bounded vectors, not vector fields. The electrodynamics module
and your vector analysis module would deal with vector fields. How can
these two notions can be combined (is it possible/practical)? I would
like to see this discussed in your proposal.

Here is an example of the issue:

1. a field
some_vortex_field = Vector((-y, x, 0), (x, y, z))

2. a bounded vector
some_force_on_some_point = Vector((1,2,0), (p_x, p_y, p_z))

The first example is what one will do in electrodynamics do define
some field. You can calculate its curl for instance.

The second example is something from the mechanics module. The point
is for instance some joint that we are studying and the vector is the
force at that point.

Does it make sense to represent both these objects with the same class?

[Sachin Joglekar]

@Stefan, I agree with you on this. For this very reason, Gilbert, Jason and I had discussed the need for a CoordinateSystem class some time back, to supplement my electrodynamics work. Moreover, working in R3 seems natural for electromagnetic fields as of now. For now, if Prasoon and I start working on our projects simultaneously, I could initially use the existing infrastructure of sympy.physics.mechanics, but once his work gets merged, I could significantly reduce my code, or infact improve it largely. Gilbert can comment on this issue and about how Prasoon and I could combine our work.

[Stefan Krastanov]

This all seems clear. Hopefully SymPy will be able to accept these projects.

In any case, I can not imagine work done on an electrodynamics module
without having a vector calculus module.

[Prasoon Shukla]

First, sorry for replying so late. Had some assignments to complete (still do :P )

Second, thanks a lot for elucidating everything in such detail. I am constantly amazed by the wonderful community response here. Thank you all for that.

Now, getting back to the topic of discussion,

On Tuesday, March 12, 2013 12:43:28 AM UTC+5:30, Stefan Krastanov wrote:

First of all, I think it would be nice to start writing a wiki page 
with the proposal as the general idea seems clear enough. 

 

 Yes. I'll start with the proposal shortly; I still need to coordinate with Sachin for his proposed electrodynamics module.

Second (and most important in this message), I would vote for 
completely forgetting about anything outside 3D. Implementations of 
vector calculus in Rn should be done (and are done to some extend) in 
the geometrical algebra and differential geometry modules (and the 
tensor canonicalization pull request by Pernici). The whole idea of a 
vector calculus module, in my opinion, is that when working in R3 you 
do not need all these generalized notions of, e.g, curl or vector 
product in terms of completely antisymmetric tensors. 

Obviously, we would love to have `rewrite_as` methods to couple it to 
the diffgeom/ga/tensor modules but that does not mean it is a good 
idea to use, for instance, the diffgeom module as a basis for your 
project. It can be done and for puritans it will indeed be a deep 
example of code deduplication/reuse and object oriented blah blah. 
However, the code will be harder to understand and the obvious notions 
of R3 calculus will be hidden behind their generalizations. 

Implementing in R3 reduces quite a bit of my work - Now I don't have to read Tensor Analysis. And yes, keeping in line with my objectives, implementing in R3 will definitely be more accessible to the typical user. The rewrite_as method, as you mentioned, sounds really good to me for module interoperability. 
 

And simply from a practical standpoint, the mechanics module is 
already based around R2/R3. It is more important to abstract that away 
than to generalize it to cases that will not be used. 

The rest of my comments are below, annotating your message. 

> @Gilbert: I have been thinking along the lines of the mechanics module and 
> as Stefan mentioned, I plan (at least right now) to take this idea forward 
> with component wise operations on vectors. For component independent vector 
> algebra, any numerical operation on any vector will inevitably get reduced 
> to operations on scalars (individual components). Therefore, I believe, that 
> a large part of component dependent vector operations will form a subset to 
> a more complete component independent vector algebra. So definitely, in the 
> long term, we will be able to extend this vector module to play nice with 
> component free vector algebra, just as Raoul mentioned. 

It is a bit unclear what you mean by "component dependent". Here are 
the two ways this can be understood: 

1. Simply using matrices without "metadata". I consider this a very 
bad idea as a matrix is just a representation of the geometrical 
entity vector and all that can be done this way is already in the 
matrix module. 
2. Using some symbolic object containing 
   - reference to the "world" in which we work 
   - reference to the reference/basis in which we represent the vector 
   - the coordinates in this basis, which is a matrix 

EDIT: I see that you actually use the second approach. Great! 

 

As you said, I am using approach number 2.

> Also, and just as you mentioned, Vector in mechanics can be extended because 
> the user does not directly interact with vector objects. This means that, at 
> least for mechanics, we won't have too many problems as far as integration 
> with the new module is concerned. Of course, this problem remains for other 
> modules (quantum, electricity-magnetism modules proposed by Sachin). For 
> those, I will interact with people in the community for their requirements 
> and hopefully, decide on an idea that keeps everyone happy, or at least, 
> which keeps most of the people happy. 

I would say, leave quantum away. Hilbert spaces in quantum mechanics 
are sufficiently different from R3 to necessitate their own 
implementation. On the other hand, I strongly suggest that this 
project you are preparing be the basis for the electrodynamics project 
by Sachin. You should definitely work together on this. 

Alright then. In any case, I seem to be lacking the appropriate maths required for full-fledged treatment of quantum physics anyway.
 

Yes of course. As a matter of fact, I had, at one point, considered the idea of implementing a RefranceFrame by using an origin of coordinates w.r.t. to world coordinates and using a basis defined at that origin (my point might not have come across very clearly here, I admit).

Also, the fact that basis changes from point to point is taken care of by CoordSystem object passed to Vector (the way I see it as of yet). The operations (grad, curl, div etc) on a Vector will again be defined based on the CoordSystem object.

A Point class sound like a good idea to me. This way, we can define methods on the class which would help out a lot. As an example, we can have a method to rewrite the point coordinates for a different reference frame.

> Taking a clue from diffgeom, we are defining a coordinate system first (with 
> 'rect' providing the basic cartesian coordinates, the definitions of which 
> are inbuilt. We can have similar functionality for other commonly used 
> coordinate systems). Then, we instantiate a vector and give it symbols to 
> use for coordinates. The dimension of the vector can be calculated 
> internally using the second n-tuple. The basis for any Rn space will be just 
> the standard basis. Then, we can have component-wise operations that produce 
> the result; in this case divergence and curl of the given vector field. 
> Also, for the n-dimension case, we can have something like: 
>>>> v = Vector(dim=7, symbol=y, coord=coord) 

Mostly agree, but drop dim!=3. And I have a lot of technical comments 
concerning the  the `symbol` argument, but they can wait for later.

 

Okay. Since we are going ahead with R3, so, dim will be 3 by default. 

> 
> Here, the coord object will carry the relations with one of the standard 
> coordinate systems and that will be used to define all the operations of 
> vector calculus.( The coord object can be generated by performing operations 
> on CoordSystem object of some standard coordinate system perhaps?) 
> 
> This will generate a vector (unbound?) where the coordinates wrt to the 
> standard basis (e1, ... , e7) are just (y1, y2, y3 ... y7). 

Concerning the very important "unbound?" question. In your first 
example `Vector` is actually a vector field, i.e. bound vector for 
physicists (there is a difference between vector field and a bound 
vector but drop it for the moment). In your second example `Vector` is 
actually a vector, not a vector field. Vector analysis cares about 
vector fields, so just completely drop the unbounded vectors from your 
proposal. They have place in linear algebra module, in matrix 
expressions module and in hilbert spaces module (the quantum module), 
not here. 

Okay. I was just experimenting a bit so, yeah, will drop the idea of unbound vectors. 

 

PS 

By the way, this question show an issue that must be discussed a lot 
before we are clear about how to proceed. The mechanics module deals 
with bounded vectors, not vector fields. The electrodynamics module 
and your vector analysis module would deal with vector fields. How can 
these two notions can be combined (is it possible/practical)? I would 
like to see this discussed in your proposal. 

Here is an example of the issue: 

1. a field 
some_vortex_field = Vector((-y, x, 0), (x, y, z)) 

2. a bounded vector 
some_force_on_some_point = Vector((1,2,0), (p_x, p_y, p_z)) 

The first example is what one will do in electrodynamics do define 
some field. You can calculate its curl for instance. 

The second example is something from the mechanics module. The point 
is for instance some joint that we are studying and the vector is the 
force at that point. 

Does it make sense to represent both these objects with the same class? 

I haven't thought of this in detail yet. But, at a first glance, my answer is yes. We can just have a property Vector.is_field that can be used to differentiate between the two different types. Also, if it's just a constant vector, then curl, div are both zero.

Another way to go about this is inherent in the initialization of the vector. The second 3-tuple is actually used for getting the 'variables' that must be used in taking, say, the curl of the vector field. If we don't give out these variables explicitly, then the data members of Vector object that holds that 3-tuple is just going to be None. That way too, we can distinguish between vector fields and constant vectors.

Anyway, as I mentioned before, things are very sketchy right now. Let me first talk with Sachin and try out mechanics a bit as well. I plan to make refinements in layers, fine tuning things till they are right. So, I'll begin writing a proposal soon, but I won't necessarily put it on the wiki just yet. Once I have something concrete, then I'd post it on the wiki. In the meanwhile, let's keep this thread for all discussions concerning the vector module.

Again, thanks a ton for all the help.

[Gilbert Gede]

I think Stefan probably did a better job with the terminology than I could do, and covered a lot of points. Unfortunately, as a mechanical engineer, my terminology is a bit sloppy on these topics.

What ReferenceFrame really does, is provide a set of orthonormal basis vectors to define other vectors with; ReferenceFrame also provides functionality for keeping track of rotations between these orthonormal bases. The concept of an inertial frame is specified later, when equations of motion are constructed and it matters (I'm not sure if that is true for EM fields?). We didn't really allow for translation of ReferenceFrame's, but instead wrote it as if there was no knowledge of translation position/velocity/acceleration of one ReferenceFrame w.r.t. another. This places the burden on the end user to correctly define velocities of points in multiple frames.

In mechanics, we treated all vectors as vector functions. We didn't really make a distinction between bound and free vectors - there didn't seem to be a "clean" way to do it.

Perhaps you could change the ReferenceFrame of mechanics into something such as OrthonormalBasis? Other basis types could then be created. A ReferenceFrame would then need to be instantiated with a basis type. This could allow multiple bases to be embedded within one ReferenceFrame (all fixed relative to each other). While I've never attempted (nor needed to nor previously thought about) this, I would imagine that you could extend this to include other types of bases (e.g. spherical) and still allow for algebraic operations between vectors defined in this way.

Example:

If we had a vector written with an orthonormal basis: v = 1 * e_x + 2 * e_y + 3 * e_z

and another written with a spherical basis: u = 4 * e_r + 5 * e_theta + 6 * e_phi

that they could be added/subtracted/dotted, e.g. u + v = 1 * e_x + 2 * e_y + 3 * e_z + 4 * e_r + 5 * e_theta + 6 * e_phi

Does that sound right?

Maybe the code could look something like this:

A = ReferenceFrame('A', 'orthonormal')

A.attach_basis('spherical', _some_initial_orientation_info_)

v = A.x + A.y

u = A.r + A.theta

Maybe it's possible to extend this to include translational information also?

[Stefan Krastanov]

> What ReferenceFrame really does, is provide a set of orthonormal basis
> vectors to define other vectors with; ReferenceFrame also provides
> functionality for keeping track of rotations between these orthonormal
> bases. The concept of an inertial frame is specified later, when equations
> of motion are constructed and it matters (I'm not sure if that is true for
> EM fields?). We didn't really allow for translation of ReferenceFrame's, but
> instead wrote it as if there was no knowledge of translation
> position/velocity/acceleration of one ReferenceFrame w.r.t. another. This
> places the burden on the end user to correctly define velocities of points
> in multiple frames.

The way I understand this is that, because `mechanics` works only with
cartesian coordinate systems for now, it does not need to distinguish
between coordinate system / basis / frame, because the coordinates are
not curved and all this notions turn out to be the same.

> Perhaps you could change the ReferenceFrame of mechanics into something such
> as OrthonormalBasis? Other basis types could then be created. A
> ReferenceFrame would then need to be instantiated with a basis type. This
> could allow multiple bases to be embedded within one ReferenceFrame (all
> fixed relative to each other). While I've never attempted (nor needed to nor
> previously thought about) this, I would imagine that you could extend this
> to include other types of bases (e.g. spherical) and still allow for
> algebraic operations between vectors defined in this way.
> Example:
> If we had a vector written with an orthonormal basis: v = 1 * e_x + 2 * e_y
> + 3 * e_z
> and another written with a spherical basis: u = 4 * e_r + 5 * e_theta + 6 *
> e_phi
> that they could be added/subtracted/dotted, e.g. u + v = 1 * e_x + 2 * e_y +
> 3 * e_z + 4 * e_r + 5 * e_theta + 6 * e_phi
> Does that sound right?

I like the general idea.

> Maybe the code could look something like this:
> A = ReferenceFrame('A', 'orthonormal')
> A.attach_basis('spherical', _some_initial_orientation_info_)
> v = A.x + A.y
> u = A.r + A.theta
>
> Maybe it's possible to extend this to include translational information
> also?

Actually this is not too different from what was done for the diffgeom
module (and the Scheme codebase that inspired it).

By the way, here with the spherical example one can see a difference
between the basis/frame and the coordinate system
- the coordinate system is not of straight lines that can be
identified with scalar*unit_vector
- the basis is of unit vectors (however these vectors change
direction from point to point)

> In mechanics, we treated all vectors as vector functions. We didn't really
> make a distinction between bound and free vectors - there didn't seem to be
> a "clean" way to do it.

The fact that it works quite well for now makes me hopeful that it
will not be necessary to distinguish them (unbounded/bounded/field)
even in the proposed new module.

[Stefan Krastanov]

PS I was somewhat loud in suggesting that the work should be focused
only on R3. I do not have any authority on this, so before sticking
this in the application it might be a good idea to check whether the
community agrees.

[Prasoon Shukla]

On Wednesday, March 13, 2013 4:00:11 AM UTC+5:30, Gilbert Gede wrote:

I think Stefan probably did a better job with the terminology than I could do, and covered a lot of points. Unfortunately, as a mechanical engineer, my terminology is a bit sloppy on these topics.

What ReferenceFrame really does, is provide a set of orthonormal basis vectors to define other vectors with; ReferenceFrame also provides functionality for keeping track of rotations between these orthonormal bases. The concept of an inertial frame is specified later, when equations of motion are constructed and it matters (I'm not sure if that is true for EM fields?). We didn't really allow for translation of ReferenceFrame's, but instead wrote it as if there was no knowledge of translation position/velocity/acceleration of one ReferenceFrame w.r.t. another. This places the burden on the end user to correctly define velocities of points in multiple frames.

In mechanics, we treated all vectors as vector functions. We didn't really make a distinction between bound and free vectors - there didn't seem to be a "clean" way to do it.

Perhaps you could change the ReferenceFrame of mechanics into something such as OrthonormalBasis? Other basis types could then be created. A ReferenceFrame would then need to be instantiated with a basis type. This could allow multiple bases to be embedded within one ReferenceFrame (all fixed relative to each other). While I've never attempted (nor needed to nor previously thought about) this, I would imagine that you could extend this to include other types of bases (e.g. spherical) and still allow for algebraic operations between vectors defined in this way.

Example:

If we had a vector written with an orthonormal basis: v = 1 * e_x + 2 * e_y + 3 * e_z

and another written with a spherical basis: u = 4 * e_r + 5 * e_theta + 6 * e_phi

that they could be added/subtracted/dotted, e.g. u + v = 1 * e_x + 2 * e_y + 3 * e_z + 4 * e_r + 5 * e_theta + 6 * e_phi

Does that sound right?

Yes. That sounds perfectly alright to me. 

But, as Stefan pointed out, we need to be very clear in taking the correct transformation of the basis vectors. In this example, the vector u defined is ambiguous because we do not know the values of e_r, e_theta and e_phi (more precisely, we also need another point in space where this vector is defined).

Indeed, we can just take the default values of e_r, r_theta, e_phi to be at origin of coordinates but I think that isn't very explicit (python zen).

Again, I'd like to reiterate my point about the CoordSystem class. This class will contain all the information required for the transformations to other coordinate systems. So, adding v and u, or taking dot product will not be a problem.

 

Maybe the code could look something like this:

A = ReferenceFrame('A', 'orthonormal')

A.attach_basis('spherical', _some_initial_orientation_info_)

v = A.x + A.y

u = A.r + A.theta

Maybe it's possible to extend this to include translational information also?

That can be one way, but, as was mentioned before, this way is too mechanics oriented. What I'd like is a more direct interface to the vectors. Also, attaching a basis will be taken care of by the attaching another CoorSystem object which is basically a wrapper to all functionalities for any coordinate system.

@Stefan: Okay then.

So, let us have a poll here. Would it be better to go for R3 or Rn?

Certainly, R3 would be easier for me. Also the concepts and applications (use cases) will come more naturally in R3 . Also, there are several other points that Stefan has mentioned before.

On the other hand, Rn would be the more general approach. I'll have some problems with Rn though. First, I'll have to read some tensor analysis. That won't take too long though. The other problem would finding use cases (since I have to make the module so that the use cases will be covered in the implementation) and making my implementation to fit the stuff already made in mechanics.

That being said, I personally favour the R3 implementation.

So, please provide your views on this.

[Sachin Joglekar]

Looking at my proposed project, its obvious that my work would be based around R3. Moreover, if Prasoon's work is based around R3, the interface of his work which would be used by me would considerably simplify. However, I guess its the decision of the sympy.physics community ultimately.We could definitely modify our work for the 'greater good'.

[someone]

> Looking at my proposed project, its obvious that my work would be
> based around R3.

Why R3, why not RN?

Ok, no curved spaces, no manifolds, we have that in the diffgeom
module.

But R3 and RN are not that different mathematically. Hence, why
make this restriction? (We did not restrict to 3x3 matrices f.e.)

Can you please explain inside your proposal why to choose these
limitations in more detail?

Personally, I don't think the argument that physics (mechanics and
electrodynamics) uses mostly R3 is a strong enough a reason.

[Sachin Joglekar]

@somebody, I have to agree with you. I just favoured R3 because MOST of my work, as of my current concept, would revolve around R3. However, if Rn is decided as the default case, it wouldnt be much of a hassle for me or Prasoon  to change our implementation accordingly. As I pointed out, its for the sympy.physics people to decide.

About the proposal, I want to study parts of the Reference book I mentioned on the thread of my proposal, before I put up my concrete proposal on the wiki. I will do that most probably by the first week of April. From what the current situation seems to be, Prasoon and I would have to coordinate our efforts about quite a few aspects of vector field handling and discuss the implementations before GSOC begins. We have started corresponding for the same, and we hope to have a good, solid idea about the entire module as a whole by mid-April. If you guys have any inputs/ ideas, we would like them to discuss them and incorporate them in our proposals by that time.

Also, I want to know what you guys think about the CoordinateFrame concept to 'extend' the idea of ReferenceFrame.

[Prasoon Shukla]

@Raoul: I still haven't chosen to restrict anything yet. It's just that Stefan had mentioned that I should stay away from things other than R3. But, Rn would obviously be more desirable (just as your post mentioned). Hence the poll.

[Aaron Meurer]

Why not a middle ground? For some things it would be silly not to support n dimensions, like dot product, but for others, like cross product, the n-dimensional generalization is more complicated, and (if I understand correctly), not even technically a vector. For those cases, you could give an error if n > 3. 

Could someone write up a list if everything that would be tricky to do in n dimensions (rather than 3)? Someone mentioned creating a wiki page. That would be a good place for this. 

By the way, you say Rn, but is there a reason to not use Cn instead? Or maybe it won't actually matter for 99% of the code. 

Aaron Meurer

On Mar 13, 2013, at 6:59 PM, Prasoon Shukla <prasoon...@gmail.com> wrote:

@Raoul: I still haven't chosen to restrict anything yet. It's just that Stefan had mentioned that I should stay away from things other than R3. But, Rn would obviously be more desirable (just as your post mentioned). Hence the poll.

--

[someone]

> Why not a middle ground? For some things it would be silly not to
> support n dimensions, like dot product, but for others, like cross
> product, the n-dimensional generalization is more complicated, and
> (if I understand correctly), not even technically a vector. For those
> cases, you could give an error if n > 3.

Yes, I implicitly excluded the wedge/cross product from my wish
to have vectors in Rn. This is where we can make to border towards
geometric/Grassmann/exterior/Clifford Algebra things. (But make sure
one can enter that world from vector expressions)

Actually, here is another cross product in R7:

http://en.wikipedia.org/wiki/Seven-dimensional_cross_product

And for R1, R2 we can trivially pad vectors to R3.
Maybe one could even fill up R4,R5,R6 to R7? But that
I never tried.

> Could someone write up a list if everything that would be tricky to
> do in n dimensions (rather than 3)? Someone mentioned creating a wiki
> page. That would be a good place for this.

Actually, I the abstract algebra code & paper I mentioned, all basic
axioms of vector algebra are valid in Rn *except* the cross product.

> By the way, you say Rn, but is there a reason to not use Cn instead?
> Or maybe it won't actually matter for 99% of the code.

Oh sure! The base field should not matter. (Probably char 0, maybe
not even necessary to assume that.) Anyway, at least R and C should
be supported, they should not be hard-coded though.

[Alan Bromborsky]

You might want to look at the book "New Foundations for Classical
Mechanics (2nd edition)" by David Hestenes, Kluwer Academic Publishers.

Especially interesting is the spinor formulation of the two body problem
which gives an harmonic oscillator in the spinor variables. This is
especially important for a perturbation analysis of the n-body problem
since perturbed harmonic oscillators have nice convergence properties.

[Prasoon Shukla]

> Why not a middle ground?

Indeed. I like this idea a lot. +1

I think that Raoul's post takes care of the requirement of the list. Also, I will still take some time to write the proposal on the wiki. First, I want to be really clear on exactly what to implement and, more importantly, how to implement.

And sure, we can implement all things on a Cn over C vector space.

Raoul also mentioned the 7 dimensional cross product.

> Maybe one could even fill up R4,R5,R6 to R7? But that 

> I never tried.

Neither have I. So, let me. Though a first glance at Cayley's sample multiplication table (in the wikipedia page mentioned) reveals that the 7 dimensional cross product may generalize all the way down to 1 dimension. This I say because the cross product of ei, i=1, 2, 3, is same as for cross product in R3. So, this looks hopeful.

Anyway, in the meantime, I have been jotting down ideas as they come to me on a whiteboard and copying them down in a notebook. Right now, I am trying to deal with coordinate systems and how best to implement them in a class. Some ideas include user defined bases, conversion matrices between different coordinate systems, conversion between different bases etc. After I have a basic visualization of the structure of the coordinate system class, I'll post my ideas here. Then, next thing I will focus on is the concept of Reference Frames and a world based coordinate origin. Until then, please feel free to post any ideas/suggestions/opinions, however small, over here.

Thanks for the valuable inputs everyone. Much appreciated.

[Prasoon Shukla]

I won't be available till Tuesday (19th March); going on a field trip. In the meanwhile, please keep posting ideas/suggestions/problems if any.

Thanks.

[Sachin Joglekar]

@Prasoon, I recently did a bit of study and planning about our coordination on vector calculus as far as our GSoC proposals are concerned. I will post my total concept of the module by first week of April (@asmeurer, is this a good enough deadline? considering this project will need lots of inputs from the sympy.physics community) but I want to clear up a few things, to make your work (and mine) easier.

1. How do you plan to implement Gauss' theorem (divergence) and Stoke's fundamental theorem of curl? I am not talking about the actual implementation here, I would like to have an idea about the API of the functions/methods that you plan to implement. I know this is a bit early for such questions maybe, but I would appreciate it if you clarify this before I upload my proposal. Just a pointer this one.

2. I would also need a look at your API for vector field functions to integrate them with some functions I plan to implement for scalar potential fields. I am not clear about this concept myself (and whether implementing scalar fields is necessary), but having methods to return vector fields from potential field gradients would a good idea. I would like the community's views on this.

3. Do you plan to implement methods for polar and cylindrical coordinates? 

4. As your work would primarily deal with functions of the coordinates themselves, we would have to collaborate to integrate time-based fields into your work (maybe an auxiliary PR to your branch from my side, or extension of your classes as a part of my work). This is more of a query to our potential mentors here. What do you people think would be the best course of action in this case?

[Aaron Meurer]

On Mar 21, 2013, at 8:02 AM, Sachin Joglekar <srjogl...@gmail.com> wrote:

@Prasoon, I recently did a bit of study and planning about our coordination on vector calculus as far as our GSoC proposals are concerned. I will post my total concept of the module by first week of April (@asmeurer, is this a good enough deadline? considering this project will need lots of inputs from the sympy.physics community) but I want to clear up a few things, to make your work (and mine) easier.

Yes. See http://www.google-melange.com/gsoc/events/google/gsoc2013. 

1. How do you plan to implement Gauss' theorem (divergence) and Stoke's fundamental theorem of curl? I am not talking about the actual implementation here, I would like to have an idea about the API of the functions/methods that you plan to implement. I know this is a bit early for such questions maybe, but I would appreciate it if you clarify this before I upload my proposal. Just a pointer this one.

+1. Figuring out the APIs is more important than figuring out the implementation for these projects. 

2. I would also need a look at your API for vector field functions to integrate them with some functions I plan to implement for scalar potential fields. I am not clear about this concept myself (and whether implementing scalar fields is necessary), but having methods to return vector fields from potential field gradients would a good idea. I would like the community's views on this.

3. Do you plan to implement methods for polar and cylindrical coordinates? 

4. As your work would primarily deal with functions of the coordinates themselves, we would have to collaborate to integrate time-based fields into your work (maybe an auxiliary PR to your branch from my side, or extension of your classes as a part of my work). This is more of a query to our potential mentors here. What do you people think would be the best course of action in this case?

You may have to help one another, especially with bug fixes. Coordinating timelines is a great idea. We can also try to make sure that the mentors are coordinated. 

One thing: you shouldn't think "send a pull request to the other's branch". All completed work should be submitted for review and merged into master as early and often as possible, so you just sending pull requests fixing what is already in main repo. 

And another note, if your project does depend on another, you should put a contingency plan in your proposal in case only you are accepted. 

Aaron Meurer

[Prasoon Shukla]

I have been really busy with some other work this week. Nevertheless, I have made some progress with the ideas. I have finalized (almost) how to implement coordinate systems, and yes, we are going to have cylindrical as well as spherical coordinates. Right now, I am trying to decide how best to implement a Vector class. Next, I will be deciding on an implementation of reference frames.

Anyway, give me until tomorrow and I'll have answers to all your questions.

[someone]

> I have finalized (almost) how to implement coordinate systems, and yes,
> we are going to have cylindrical as well as spherical coordinates.

Great! But do not hard code the choice of available systems. There are
many many more useful coordinates:

http://mathworld.wolfram.com/OrthogonalCoordinateSystem.html

From the second-last paragraph:

"Three-dimensional orthogonal curvilinear coordinate systems of degree two or less include
bipolar cylindrical coordinates, bispherical coordinates, three-dimensional Cartesian coordinates,
confocal ellipsoidal coordinates, confocal paraboloidal coordinates, conical coordinates,
cyclidic coordinates, cylindrical coordinates, elliptic cylindrical coordinates,
oblate spheroidal coordinates, parabolic coordinates, parabolic cylindrical coordinates,
paraboloidal coordinates, prolate spheroidal coordinates, spherical coordinates,
and toroidal coordinates."

At least some of the most important elliptic and parabolic systems
for physics should be added too. (Maybe not as part of your GSOC but
later.) Its not that clear where to put to border towards diff geom now.

[Prasoon Shukla]

@Sachin, here are some of the answers you sought.

>1. How do you plan to implement Gauss' theorem (divergence) and Stoke's fundamental theorem of curl? I >am not talking about the actual implementation here, I would like to have an idea about the API of the >functions/methods that you plan to implement. I know this is a bit early for such questions maybe, but I would >appreciate it if you clarify this before I upload my proposal. Just a pointer this one.

I have not quite gotten to the point where I would think about integration of vector fields. Right now, I am still trying to finalize an idea for implementing a Vector class. I am pretty clear about how it will be implemented; just some more fine tuning remains. As I had mentioned before, I'll be taking up the idea for reference frames and world coordinate system next; only after that, I'll move on to standard vector calculus operations. However, in the meantime, I can think of a few things regarding your query.

Clearly, the integration for vector fields is going to be implemented in a fundamentally different way; in the Vector class at least. This class will convert the integrals into integrals over scalar functions and those can be further evaluated by SymPy calculus core. Anyway, I think that we will have an Integral class that we can use to represent integration over vector field type integrands. Now, in this Integral class, we can implement the basic theorems (Gauss, Stokes). So, we can have something like:

>>> x, y, z = symbols('x y z')

>>> c = CoordSys(basis_var='e', dim=3, type='rect')

>>> field = Vector( (y**2 + x, x**2 - y, x*y*z), (x, y, z), coord=c)

>>> integral = field.integrate( <values describing integral type - line, surface, vol etc and ranges>)

# this is very rough and sketchy yet

>>> integral.convert_gauss()

Vector.Integral

So, the last line returns a new Vector.Integral object with changed parameters - volume integral changers to surface for example. Again I'd like to emphasize this: I have not gotten to thinking about this part in detail. So, what I wrote above was just something that I though of quickly and there is no reason to believe that this is what things would finally look like.

And of course, since I haven't made this part concrete yet, I think I can accommodate your demands without much difficulty.

Coming to point #2, yes, we will have a gradient function that can be used to return the grad of a vector field. For the Vector API: I have given thought to how vector fields will work internally. But, there's no API design as of yet. Give me a two-three days and I'll post the APIs for CoordSys (coordinate system) and the Vector class (which depends on a coordinate system) here.

#3. Yes. I mentioned it in my last post. And as Raoul said, we can also have a few more useful coordinate systems implemented too.

#4. Time-dependent fields. Ah now that is something I haven't really thought about. But, as I see it, this shouldn't pose too many problems. We can just have one more variable, say t, that we can use in definition of vector fields. Just that this variable will not be differentiated with respect to the coordinates. So I think that the Vector API will be such that this will be automatically taken care of without needing to meddle with the Vector class itself.

I hope that those answers were satisfactory for the time being. Also, I again apologize for the slow pace; March has proven to be specially busy. We have some holidays next week here in India so hopefully, I'll make some good progress in that time.

[Sachin Joglekar]

Sure. As it is, I would be start finalising my idea and concepts next week after my quizzes end. This was more of a mental note to get these things clarified from your side. And about time dependent fields, I guess I will implement them as a part of my PR, which would extend your classes or use your API, cause time doesn't factor into vector calculus as such.

[Prasoon Shukla]

Okay, so I have been away a while due to a lot of reasons. Fortunately, I will have time from now on for about 3 weeks. As of now, I have a working API of the Coordinate System class, the Vector class, and, the Reference Frame class. The Reference Frame class is where I will try to implement all of the time dependent functionality. I will explain all of these in detail soon in my proposal on the wiki.

The next thing that I am trying to think of is the Vector.integrate() method. As is obvious, we will use it for integration over vector and scalar fields.

Now, the first thing is to have a VectorIntegral class that would represent an integral complete with the field, the type (line, area or volume) and the path/area/volume to integrate over. The problem I am having now is this: Right now, SymPy can solve integrals symbolically for the one-dimensional case. I need to implement

a) Line integrals - For this I need some notion of a path. Right now, the integrate function in SymPy takes only straight paths (that is we can only provide end points and the path is just one of the coordinate axes). But, I need support for three dimensional paths. One solution I thought of is to have a class called Path that would contain the parametric definitions for the path. The methods in this class can further be used as helpers to evaluate the integral. This can be done by reducing the line integral into simple integral in one variable which SymPy can already evaluate.

b) Surface Integrals - Again, I need some way to represent an area. So, perhaps a Surface class will do the job? Implementation will be similar as mentioned above.

c) Volume Integrals - Needs a Volume class. Similar implementation.

The implementation of a Path, Surface and a Volume class seems to be the solution to this problem at this point. If anyone else has a better idea, please do tell. Otherwise, I think I'll go ahead with this.

Another problem that I have right now is this: How do I calculate the multiple integrals that will arise in evaluating the above three cases? I think that this can be implemented with some work in the integration module, which I'll have to do at a later stage.

But then again, SymPy is a symbolic maths library, and, the computations of the above mentioned integrals will be numerical in nature. So, I think that we may not need too many major changes to be able to implement this. For example, line integrals can just be reduced to integrals in one variable. Surface and volume integrals will need to be reduced to simple integrals in 2 and 3 dimensions and then be evaluated.

Anyway, these are just the ideas in my mind. And admittedly, they are in a crude form. Hopefully, in a day or two, I'll have more refined ideas to go on with. In the meantime, please make do tell any suggestions/changes that you gentlemen think should be there.

Thanks

[someone]

> a) Line integrals - For this I need some notion of a path. Right now,
> the integrate function in SymPy takes only straight paths (that is we
> can only provide end points and the path is just one of the
> coordinate axes). But, I need support for three dimensional paths.
> One solution I thought of is to have a class called Path that would
> contain the parametric definitions for the path. The methods in this
> class can further be used as helpers to evaluate the integral. This
> can be done by reducing the line integral into simple integral in one
> variable which SymPy can already evaluate.

I considered doing this for contour integrals in the complex plane
some long time ago. If you do it the "right" way, then the code could
maybe be reused for that case too. Although the typical paths in contour
integrals are somewhat specific combinations of circles and lines.
(I hoped that I could parametrize this all with the usual epsilons.)

> b) Surface Integrals - Again, I need some way to represent an area.
> So, perhaps a Surface class will do the job? Implementation will be
> similar as mentioned above.
>
> c) Volume Integrals - Needs a Volume class. Similar implementation.
>
> The implementation of a Path, Surface and a Volume class seems to be
> the solution to this problem at this point. If anyone else has a
> better idea, please do tell. Otherwise, I think I'll go ahead with
> this.

I have the hope that one could get away with just a single parametric
object depending on as many parameters as necessary:

L(u) := L(u_1, u_2, ..., u_n)

would then represent an integral over an N-dimensional hypervolume
(path: n=1, surface: n=2, volume: n=3 and so on) and the integral

Int(f , L) --> Int( f(L(u)) |J L| du
--> Int(...Int( f(L(u_1, ..., u_n)) |J L| du_1, ..., du_n

It is for sure not easy to do this all the correct way such
that all fits well together. But I think it's worth thinking
about it!

There is the most general version of the Stoke's theorem:

http://upload.wikimedia.org/math/8/e/4/8e4e3b1cd7d06acf2438e40509c4f4e7.png

in this notation independent on dimension of the "volume" Omega
and its border dOmega. From this we can get all others, and that
is the point why I think it should be possible to do this in a
mostly dimension agnostic way.

[Aaron Meurer]

On Mon, Apr 1, 2013 at 9:10 AM, Prasoon Shukla <prasoon...@gmail.com> wrote:

Okay, so I have been away a while due to a lot of reasons. Fortunately, I will have time from now on for about 3 weeks. As of now, I have a working API of the Coordinate System class, the Vector class, and, the Reference Frame class. The Reference Frame class is where I will try to implement all of the time dependent functionality. I will explain all of these in detail soon in my proposal on the wiki.

The next thing that I am trying to think of is the Vector.integrate() method. As is obvious, we will use it for integration over vector and scalar fields.

Now, the first thing is to have a VectorIntegral class that would represent an integral complete with the field, the type (line, area or volume) and the path/area/volume to integrate over. The problem I am having now is this: Right now, SymPy can solve integrals symbolically for the one-dimensional case. I need to implement

a) Line integrals - For this I need some notion of a path. Right now, the integrate function in SymPy takes only straight paths (that is we can only provide end points and the path is just one of the coordinate axes). But, I need support for three dimensional paths. One solution I thought of is to have a class called Path that would contain the parametric definitions for the path. The methods in this class can further be used as helpers to evaluate the integral. This can be done by reducing the line integral into simple integral in one variable which SymPy can already evaluate.

Take a look at the already existing line_integrate() as well.

 

b) Surface Integrals - Again, I need some way to represent an area. So, perhaps a Surface class will do the job? Implementation will be similar as mentioned above.

c) Volume Integrals - Needs a Volume class. Similar implementation.

The implementation of a Path, Surface and a Volume class seems to be the solution to this problem at this point. If anyone else has a better idea, please do tell. Otherwise, I think I'll go ahead with this.

Another problem that I have right now is this: How do I calculate the multiple integrals that will arise in evaluating the above three cases? I think that this can be implemented with some work in the integration module, which I'll have to do at a later stage.

But then again, SymPy is a symbolic maths library, and, the computations of the above mentioned integrals will be numerical in nature. So, I think that we may not need too many major changes to be able to implement this. For example, line integrals can just be reduced to integrals in one variable. Surface and volume integrals will need to be reduced to simple integrals in 2 and 3 dimensions and then be evaluated.

How are they numerical in nature? Integrals are by their nature symbolic. They take symbolic functions and return a symbolic result. Numerics are useful if the answer is not expressible in terms of nice functions. 

 

Aaron Meurer

Anyway, these are just the ideas in my mind. And admittedly, they are in a crude form. Hopefully, in a day or two, I'll have more refined ideas to go on with. In the meantime, please make do tell any suggestions/changes that you gentlemen think should be there.

Thanks

--

You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.

Visit this group at http://groups.google.com/group/sympy?hl=en-US.

[Aaron Meurer]

On Mon, Apr 1, 2013 at 5:15 PM, someone <some...@bluewin.ch> wrote:

> a) Line integrals - For this I need some notion of a path. Right now,
> the integrate function in SymPy takes only straight paths (that is we
> can only provide end points and the path is just one of the
> coordinate axes). But, I need support for three dimensional paths.
> One solution I thought of is to have a class called Path that would
> contain the parametric definitions for the path. The methods in this
> class can further be used as helpers to evaluate the integral. This
> can be done by reducing the line integral into simple integral in one
> variable which SymPy can already evaluate.

I considered doing this for contour integrals in the complex plane
some long time ago. If you do it the "right" way, then the code could
maybe be reused for that case too. Although the typical paths in contour
integrals are somewhat specific combinations of circles and lines.
(I hoped that I could parametrize this all with the usual epsilons.)

That would be neat. Of course, the real work there would be getting the residue() function working better, so that you could actually compute such integrals.  But this is a GSoC idea onto itself.

 

> b) Surface Integrals - Again, I need some way to represent an area.
> So, perhaps a Surface class will do the job? Implementation will be
> similar as mentioned above.
>
> c) Volume Integrals - Needs a Volume class. Similar implementation.
>
> The implementation of a Path, Surface and a Volume class seems to be
> the solution to this problem at this point. If anyone else has a
> better idea, please do tell. Otherwise, I think I'll go ahead with
> this.

I have the hope that one could get away with just a single parametric
object depending on as many parameters as necessary:

  L(u) := L(u_1, u_2, ..., u_n)

would then represent an integral over an N-dimensional hypervolume
(path: n=1, surface: n=2, volume: n=3 and so on) and the integral

  Int(f , L)  -->  Int( f(L(u)) |J L| du
              -->  Int(...Int( f(L(u_1, ..., u_n)) |J L| du_1, ..., du_n

It is for sure not easy to do this all the correct way such
that all fits well together. But I think it's worth thinking
about it!

Manifolds in general have to be parameterized in patches (any manifold that is not diffeomorphic to R^n will by definition require more than one patch of differential parametrization to describe). The diffgeom module should take care of this sort of stuff, though.

 

Aaron Meurer 

There is the most general version of the Stoke's theorem:

http://upload.wikimedia.org/math/8/e/4/8e4e3b1cd7d06acf2438e40509c4f4e7.png

in this notation independent on dimension of the "volume" Omega
and its border dOmega. From this we can get all others, and that
is the point why I think it should be possible to do this in a
mostly dimension agnostic way.

[Prasoon Shukla]

I took a look at the line_integrate() function. The Curve class is something of a start. Though, I like very much Raoul's idea for a single class to represent everything (path/surface/volume). I think this is doable. But,

> would then represent an integral over an N-dimensional hypervolume

I have not had a relevant course that would acquaint me with 'N-dimensional hypervolumes'. Nevertheless, I am willing to read up on it. I hope that just implementing integrals in such hypervolumes wouldn't be too different from the 2 or 3 dimensional case.

> How are they numerical in nature? Integrals are by their nature symbolic. They take symbolic functions and   > return a symbolic result. Numerics are useful if the answer is not expressible in terms of nice functions. 

 

Indeed they are. I said numeric from the user's perspective. For example, if a user wants to calculate a line integral, say, then unlike integration in one variable where we could do indefinite integrals and expect the answer to be in terms of elementary functions, here, we might be just concerned with the numeric result of the integration. Of course, giving the output in terms of elementary functions would be much better but I am inclined to think that we will probably need rely heavily on some kind of a numeric backend.

[Aaron Meurer]

Don't concern yourself with whether or not integrate can give a
symbolic answer or not. That's somebody else's problem. If it can do
it, it will return the result. Otherwise, it will return an
unevaluated Integral, which can then be numerically evaluated.

Aaron Meurer

[someone]

> That would be neat. Of course, the real work there would be getting
> the residue() function working better, so that you could actually
> compute such integrals. But this is a GSoC idea onto itself.

Yes, sure. There was no focus on "residue", just on defining paths.
And this is the only point where these two projects meet.

> Manifolds in general have to be parameterized in patches (any
> manifold that is not diffeomorphic to R^n will by definition require
> more than one patch of differential parametrization to describe). The
> diffgeom module should take care of this sort of stuff, though.

Sure. And this makes the border towards diffgeo unclear.
For simple cases the method with one patch still works.

Maybe these integration functions should work with the objects
from the diffgeo module to specify hypervolumes ((closed) orientable
manifolds) for the complicated cases.

[someone]

> I took a look at the line_integrate() function. The Curve class is
> something of a start. Though, I like very much Raoul's idea for a
> single class to represent everything (path/surface/volume). I think
> this is doable.

The idea is to push this as far as possible before
switching to the heavy diffgeo stuff.

As mentioned, this breaks down when we need multiple
patches to cover the manifold.

> > would then represent an integral over an N-dimensional hypervolume
>
> I have not had a relevant course that would acquaint me with
> 'N-dimensional hypervolumes'. Nevertheless, I am willing to read up
> on it. I hope that just implementing integrals in such hypervolumes
> wouldn't be too different from the 2 or 3 dimensional case.

No, the principle is mostly the same, just there are more
variables. There is nothing fundamentally new in higher dimensions.

> > How are they numerical in nature? Integrals are by their nature
> > symbolic.
> They take symbolic functions and > return a symbolic result.
> Numerics are useful if the answer is not expressible in terms of nice
> functions.

Well, numerics is another topic. If we can not solve the integrals,
return them symbolically as unevaluated placeholder objects.
In the best case, the users then can call .evalf() and get a number.

I'd not make numerics part of your proposal. There is more than
enough work with the vector calc stuff.

[Prasoon Shukla]

I think I have a fair idea of what the class may behave like. Let's have a vector field like so (an example):

v(x, y, z) = f (x, y, z) i + g (x, y, z) j + h (x, y, z) k

Just and example in rectangular coordinates using variables x, y and z.

Now, we can have a class, say ParamSpace(). Let us define a spherical surface like so:

p = ParamRegion(x=5*sin(u)*cos(v), y=5*sin(u)*sin(v), z=5*cos(u), params=[u, v], bounds=[(u, 0, PI), (v, 0, 2*PI) ])

This object can then be passed to Vector.integrate() and can serve as the container for defining the region of integration.

Obviously, one problem is how to implement multiple patches. I am still thinking on it.

After this, the next step is the integration itself. This means I need to add functionality to the integration module so that it can handle multiple integrals too. (I couldn't find any such functions already existing in SymPy. Is there any implementation of it?).

Another problem is with vector calc theorems. For example, if you have a contour for a line integral, what area do you choose for applying the Stokes theorem (as the area can be any that has the contour as its boundary)? Similarly for Gauss theorem too. I am still thinking in this direction. Will post here when something occurs to me.

[Alan Bromborsky]

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy?hl=en-US.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

You might want to let i,j, and k be sympy noncommuting symbols then you automatically get vector addition, subtraction, and multiplication by a scalar and could write -

p = ParamRegion(r=5*sin(u)*cos(v)*i+5*sin(u)*sin(v)*j+5*cos(u)*k, params=[u, v], bounds=[(u, 0, PI), (v, 0, 2*PI) ])

change of coordinates becomes proper substitutions for i,j,k in terms of i',j',k'.

Practically speaking with stokes theorem a surface would be defined first and the closed contour would lie on the surface, not the other way around.

[Alan Bromborsky]

Another suggestion would be to overload the multiplication operator (__mul__) so that if a and b are vectors then a*b = a.b+a x b.  This could be accomplished with a
multiplication table for the basis vectors i,j,k and using the sympy member functions args_cnc() and subs(multiplication_table_dictionary) to perform the multiplications.  Thus for a*b the scalar part would always be the dot product and the non-commutative part the vector product.

[Aaron Meurer]

On Apr 4, 2013, at 4:27 AM, Prasoon Shukla <prasoon...@gmail.com> wrote:

I think I have a fair idea of what the class may behave like. Let's have a vector field like so (an example):

v(x, y, z) = f (x, y, z) i + g (x, y, z) j + h (x, y, z) k

Just and example in rectangular coordinates using variables x, y and z.

Now, we can have a class, say ParamSpace(). Let us define a spherical surface like so:

p = ParamRegion(x=5*sin(u)*cos(v), y=5*sin(u)*sin(v), z=5*cos(u), params=[u, v], bounds=[(u, 0, PI), (v, 0, 2*PI) ])

This object can then be passed to Vector.integrate() and can serve as the container for defining the region of integration.

Obviously, one problem is how to implement multiple patches. I am still thinking on it.

After this, the next step is the integration itself. This means I need to add functionality to the integration module so that it can handle multiple integrals too. (I couldn't find any such functions already existing in SymPy. Is there any implementation of it?).

Do you mean integrate(f(x, y), (x, a, b), (y, c, d))?

Aaron Meurer

 

Another problem is with vector calc theorems. For example, if you have a contour for a line integral, what area do you choose for applying the Stokes theorem (as the area can be any that has the contour as its boundary)? Similarly for Gauss theorem too. I am still thinking in this direction. Will post here when something occurs to me.

--

[Prasoon Shukla]

> Do you mean integrate(f(x, y), (x, a, b), (y, c, d))?

Yes, except that the limits may not be constant. I looked at the docstring of integrate() and didn't see any such examples. Perhaps I missed something?

[Aaron Meurer]

The limits can be anything. It doesn't even care if they depend on the
integration variable. It will just do the integrals iteratively until
it can't, substituting the limits as it goes. There are currently no
algorithms implemented to do any kind of special handling of multiple
integrals (like reordering the limits or rewriting them).

Aaron Meurer

[Prasoon Shukla]

Alright, I think most of the ideas are clear at this point. With that, I'll start writing the proposal and hopefully post it by Monday. Then, we can continue refinements based on suggestions from the community. And again, I'd like to say a big thanks to everyone who has helped me out here. Thanks guys!

[Prasoon Shukla]

I am still in the process of writing the proposal. I just want to clarify all points in it and hence, it is taking somewhat longer.

@Aaron: In the meantime, can you tell me who might be my mentor should my application get accepted? I would like to get in touch with them on the IRC.

[Aaron Meurer]

I really can't say. It will be whoever volunteers to do it.

Aaron Meurer

[Prasoon Shukla]

So, after spending a couple of days writing the proposal, I've uploaded on to the wiki.

https://github.com/sympy/sympy/wiki/GSoC-2013-Application-Prasoon-Shukla:-Vector-Calculus-Module

@All community members: Please give it a read. This is the first draft of the proposal so it's bound to be rough-ish. Please point out things that you don't like or would like to see improved. Also, please suggest any additions that you'd like to see.

Thank you.

[Alan Bromborsky]

> --
> You received this message because you are subscribed to the Google
> Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy?hl=en-US.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

Analogous to term pythonic code I would use sympythonic code and make
the following suggestions -

Start by defining you basis vectors as noncommutative symbols -

(e1,e2,e3) = symbols('e_1 e_2 e_3',commutative=False)

Then if a1, a2, and a3 are commutative sympy expressions (symbols) any
vector a is -

a = a1*e1+a2*e2+a3*e3

Then you automatically get vector addition, subtraction, and scalar
multiplication (if c, b1, b2, and b3 are scalars)

b = b1*e1+b2*e2+b3*e3

and

a+b = (a1+b1)*e1+(a2+b2)*e2+(a3+b3)*e3

a-b = (a1-b1)*e1+(a2-b2)*e2+(a3-b3)*e3

c*a = c*a1*e1+c*a2*e2+c*a3*e3

Then if you define dictionaries for the dot and cross products, dot product

dot_dict =
{e1**2:1,e1*e2:0,e1*e3:0,e2*e1:0,e2**2:1,e2*e3:0,e3*e1:0,e3*e2:0,e3**2:1}

then -

dot(a,b) = (a*b).subs(dot_dict)

cross product

cross_dict =
{e1**2:0,e1*e2:e3,e1*e3:-e2,e2*e1:-e3,e2**2:0,e2*e3:e1,e3*e1:e2,e3*e2:-e1,e3**2:0}

cross(a,b) = (a*b).subs(cross_dict)

Like wise for coordinate transformations -

Let the bases be e1,e2,e3 ang g1,g2,g3 be related by g1 = f1(e1,e2,e3),
g1 = f2(e1,e2,e3), g1 = f2(e1,e2,e3)
where f1, f2, and f3 could also be functions of the coordinates then

g_to_e_dict = {g1:f1,g2:f2,g3:f3}

a = a1*g1+a2*g2+a3*g3

a_in_terms_of_e = a.subs(g_to_e_dict)

etc.

[Aaron Meurer]

Some of your examples in the interface section don't make any sense.
For example:

>>> w = 5*v
w

>>> f = w - 4*v
>>> f is v
True

Aaron Meurer

[Prasoon Shukla]

On Saturday, April 13, 2013 5:26:16 AM UTC+5:30, Aaron Meurer wrote:

Some of your examples in the interface section don't make any sense.
For example:

>>> w = 5*v
w

>>> f = w - 4*v
>>> f is v
True

Just trying to show vector addition and comparison. f = 5v-4v = v (mathematically). Using the 'is' was wrong of course, I admit. We can overload the == operator for the check. So,

>>> f == v

True

Any other inconsistencies that you found?

[Prasoon Shukla]

Thanks a lot. I'll make the changes soon. 

[Aaron Meurer]

On Fri, Apr 12, 2013 at 8:07 PM, Prasoon Shukla
<prasoon...@gmail.com> wrote:
>
>
> On Saturday, April 13, 2013 5:26:16 AM UTC+5:30, Aaron Meurer wrote:
>>
>> Some of your examples in the interface section don't make any sense.
>> For example:
>>
>> >>> w = 5*v
>> w

I actually meant this statement too. Assignment in Python does not
print the name of the variable assigned to.

Aaron Meurer

>>
>> >>> f = w - 4*v
>> >>> f is v
>> True
>
>
> Just trying to show vector addition and comparison. f = 5v-4v = v
> (mathematically). Using the 'is' was wrong of course, I admit. We can
> overload the == operator for the check. So,
>
>>>> f == v
> True
>
> Any other inconsistencies that you found?
>

[someone]

> Calculating Gradient
>
> A grad() function will be present in the namespace of Vector module.

Make sure this generalizes to the Jacobian Matrix for functions R^n -> R^M.

And do not forget about Hessians too.

[Stefan Krastanov]

>> Calculating Gradient
>>
>> A grad() function will be present in the namespace of Vector module.
>
> Make sure this generalizes to the Jacobian Matrix for functions R^n -> R^M.

More like a Jacobian operator. I believe the module is about vectors
as geometrical entities, not matrices.

> And do not forget about Hessians too.

I think that Hessians will be of very little use in this module, as
they are very clumsy to represent geometrically. Nonetheless they
would be an amusing addition at the later stages.

[Prasoon Shukla]

Thanks for all your responses. My semester-end exams will be starting soon and so, I won't be able to devote quite as much time as I would like to this proposal. Nevertheless, I'll make the suggested changes.

Also, I will take a look into the vector analysis packages of other CAS, Mathematica in particular. Accordingly, I'll modify the proposal for whatever functionality/implementation details that I have missed so that my proposal is more complete.

I'll report back here once the changes are done.

[Prasoon Shukla]

On Saturday, April 13, 2013 2:40:52 AM UTC+5:30, brombo wrote:

First, I thank you for the suggestion sir.

Certainly, this methods will work. But in our case, a Vector object is not merely a sum of c*e like terms. Indeed, my whole project hinges on the segregation of the two elements required to represent a vector : The coordinate system providing a way to represent points in space and also providing bases and other relevant relations, and, the vector object that will hold the components of a vector field. Also, this project is based on a component based approach to vectors. Therefore, even though this method, namely, substituting from a dictionary, works, still it seems at this point that I will need something different.

Also we will need to implement the vector calculus operations (div, curl), which wouldn't work in such a manner. What I think is that we should implement a procedure that will take care of both the dot, cross products and the div, curl. (Similar to what happens when we use the del operator). But perhaps that is a long shot. And also, this operation needs to work on individual components; the method you provide, however elegant it might be, still works on the vector as a whole. So, IMO, we should probably take this problem on with a different approach. But that's just what I think.

> Make sure this generalizes to the Jacobian Matrix for functions R^n -> R^M.

I have made some changes under the Gradient heading. I hope that takes care of the Jacobian matrix. Hessian is on hold for a bit, according to what Stefan said. But if I have some time towards the end, then certainly I shall try to implement the Hessian matrix as well.

@Stefan : A couple of days back, I asked on this thread who might be the mentor for this project, should it be accepted. From our interactions on this thread, I gather that you have a pretty strong mathematical knowledge as concerns my proposal. Also, I feel that you can be a very good guide; both for Python related things, and, the project implementation in general. You also seem to understand my level of mathematical and programming competence fairly well. So, if my project idea does get accepted, I was hoping that you would be willing to mentor for me. Would that be alright?

I understand that this is a pretty big request but I believe that you'd be a very good mentor and hence, I just had to ask.

[Stefan Krastanov]

You will need a smart way to represent dual vectors and scalar
products to represent Hessians. Even if some of us consider Hessians
less important, dual vectors and scalar products should be very high
on your list (and btw Hessians become trivial when you have the
aforementioned objects).

>
> @Stefan : A couple of days back, I asked on this thread who might be the
> mentor for this project, should it be accepted. From our interactions on
> this thread, I gather that you have a pretty strong mathematical knowledge
> as concerns my proposal. Also, I feel that you can be a very good guide;
> both for Python related things, and, the project implementation in general.
> You also seem to understand my level of mathematical and programming
> competence fairly well. So, if my project idea does get accepted, I was
> hoping that you would be willing to mentor for me. Would that be alright?
>
> I understand that this is a pretty big request but I believe that you'd be a
> very good mentor and hence, I just had to ask.

It depends on a lot of things (and no body knows how many slots will
be provided by Google), but indeed this is one of the projects in
which I might have enough math/compsci background to volunteer to
mentor (if accepted, etc, etc).

[Prasoon Shukla]

Thanks for clarifying that. And sorry for the late reply. I am buried to the neck in academics and will stay that way till about 10th May so forgive me if I cannot reply soon.

About dual vectors; I haven't read about dual vectors as of yet. So I looked around on Google and found an introductory link, http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/other/dualVectors/ . Now, if only I had some time on my hands, I could have gone to the library, found a book and read up more on the topic. Unfortunately, I do not have any time till the 10th of May and hence I won't be able to modify my application to accommodate dual vectors and how to represent dual vectors.

Nevertheless, as soon as I get free, I will definitely study the topic. Since the final results for GSoC will be announced on 27th May, so I am hoping that I can use the time before that date to read up on dual vectors and then I'll be more clear on how to implement them in code.

Also, I will submit the application today on the melange gsoc portal.

[Stefan Krastanov]

> About dual vectors; I haven't read about dual vectors as of yet.

If you wish you can think of the dual of a vector (a one-form) as a
one row matrix (and the vector as a one column matrix).

If you have a scalar product operation defined there is a canonical
way to get a dual of a vector.

For the vector v and scalar product <,>, the oneform for the vector is
v' such that v'(v) = <v, v>

[Prasoon Shukla]

Okay so let me spell out what I have understood from your comment above, just to be clear.

If you wish you can think of the dual of a vector (a one-form) as a
one row matrix (and the vector as a one column matrix).

So, if I have column vector v = [ a1, a2, ... , an ] '  (using ' for transpose), then the dual of v is then a row vector. 

 

If you have a scalar product operation defined there is a canonical
way to get a dual of a vector.

For the vector v and scalar product <,>, the oneform for the vector is
v' such that v'(v) = <v, v>

So, let's say H is a nxn positive definite Hermitian matrix, and let us define the inner product < > of vectors a and b as of dimension n:

<a, b> = Y* H X   (where Y and X coordinate matrices of vectors a and b in some ordered basis)

Then, from your comment I understand that the dual of a vector v is v' = <v, v> = V* H V (where V is the coordinate matrix). So it's just an inner product then?

But, before you mentioned that the dual of a vector is row matrix, but as far as I understand, the inner product will just be a scalar, not a row vector. So, where is it that I am not getting you correctly? Also, I am assuming that there is a standard Hermitian matrix that is used here, like maybe the identity matrix perhaps?

Also, I tried searching a bit more on Google and I found a lot of mention of the dual space, L(V, F) [ V is a vector space over the field F] which is why I am forced to believe that there might be some relation between these two (dual of a vector and dual space). Can you elaborate on this?

Again I'd like to say that I am really busy with academics for next 3 weeks. So, if my replies seem to reflect thoughtlessness on my part, please excuse me for that. I just am really badly occupied.

[Stefan Krastanov]

Don't worry, we understand that many of the applicants have exams at the moment.

For simplicity I will talk about vector spaces only in a finite number
of dimensions.

Consider a vector space V and a vector v in that space. Representing v
as a matrix in some basis results in a one column matrix.

Each vector space has a dual space V' which is the space of all
entities that take a vector from V and return a scalar. Consider v' an
element of V'. v' acting on v is written as v'(v) - a scalar.
Representing v' as a matrix in some basis results in a one row matrix
and v'(v) is just the matrix product of the two.

I will write repr(v') and repr(v) for the representations of v' and v
in the basis in question.

To say that v' is the dual of v with respect to the scalar product <,>
means that v'(v)=<v,v>.

Writing this down with matrices would be
repr(v')*repr(v)=repr(v).T*H*repr(v) where H is the hermitian matrix
corresponding to the scalar product <,>.

Indeed often the basis is chosen such that H is the identity matrix,
which means that repr(v).T=repr(v') (i.e. transposing the matrix
representing v results in the matrix representing the dual of v)

Hopefully I did not introduce any mistakes above.

[Prasoon Shukla]

I understand now. Thanks for the clear explanation.

So, let's talk about functionality that we want to give the user. In the module I have proposed, the CoordSystem class provides a method to initialize a user-defined basis, let us say B.

Now, let us define a vector v in this coordinate system. In code, we will just pass to the Vector class this CoordSystem object. So, now we have a vector v defined whose components in the basis B are [ a1, a2, ... , an] transpose.

Now, as you said, we need a way to represent the inner product. For this, we can have a class, InnerProduct. To initialize this class, the user will have to provide a positive definite Hermitian matrix which will be used to implement the inner product. The other way we can do this is for user to provide a definition for <v1, v2>. Then we can check whether this definition conforms to the standard rules of the inner product.

Once we have a well defined InnerProduct object with us, then we can move on to finding the dual of v represented in the basis B. So, our objective is to find out the row vector v'. Now, since

repr(v')*repr(v) = <v, v> = repr(v).T * H *repr(v)

So, if we assume repr(v') = [ b1, b2, .. , bn], then the LHS will be just:

[ a1 b1 + a2 b2 + .. + an bn]

and RHS is, for a given H and v, just a scalar.

So, as I see it, we have n-1 free variables and one dependent variable ( b1 through bn are all the variables here).

Am I right thus far? If I am, then we can just set the n-1 variables to be unity ( in fact, zero seems a simpler choice) and the remaining variable can be found from the relation I mentioned above.

Is this what you wanted? Or am I off base?

[Prasoon Shukla]

Also, I have added the proposal to the melange portal. Here:

http://www.google-melange.com/gsoc/proposal/review/google/gsoc2013/prasoon2211/1

[Aaron Meurer]

 Since you are taking about inner products, you might look at generalizing some of the code in the quantum module so that it just uses the vector module. Some notes:

- the quantum module uses infinite dimensional vector spaces, so it will need to work more symbolically than just rewriting everything as matrices. But you can still abstract away a lot of code this way. I haven't read your proposal in depth, so this may not be within its scope, but I would encourage you to consider it. 

- physicists represent inner products different from mathematicians. They use bras and kets, and more importantly, the conjugation convention is reversed. So the base classes would need to be ignorant of things like this, so that the main SymPy subclass can use the standard math convention and the physics subclass can use the physics convention. 

By the way, real vector spaces and complex vector spaces are fundamentally different. So far you seem to be using only real, since I haven't noticed any mention of conjugation. In general, you can actually use any field, but that is probably too algebraicmformwhat you are trying to implement.

Aaron Meurer

On Friday, April 26, 2013, Prasoon Shukla wrote:

Also, I have added the proposal to the melange portal. Here:

http://www.google-melange.com/gsoc/proposal/review/google/gsoc2013/prasoon2211/1

--

[Prasoon Shukla]

> Since you are taking about inner products, you might look at generalizing some of the code in the quantum module so     > that it  just uses the vector module.

Hmm .. I think it can be done but thinking it through and going through the code in quantum module will take some time. My final exams start in 2 days and will continue past 3rd May. Is it okay if I discuss those ideas here on the mailing list after 9th of May (which will be after the application deadline obviously) and still get them considered as part of my proposal? The final allotment will be on 27th of May so I suppose before this time, I can think more about this and solidify whatever ideas I get.

> By the way, real vector spaces and complex vector spaces are fundamentally different. So far you seem to be using only    > real, since I haven't noticed any mention of conjugation. In general, you can actually use any field, but that is probably too   > algebraic form what you are trying to implement.

Fundamentally different? How? I thought that just the field over which the vector space is defined changes in this case. Obviously C is algebraically closed so that might mean something here but I am probably missing the point.

[Aaron Meurer]

Symbolically, you can usually just consider R as a special case of C, but note that the dual spaces are different for the different fields. Actually, I'll have to think if it will make a difference for the symbolic manipulation case, or if it just matters in proofs of theorems. 

Aaron Meurer

[someone]

Hi,

> Symbolically, you can usually just consider R as a special case of C,
> but note that the dual spaces are different for the different fields.

You also have to be very careful when adding inner products.

For R: <u,v> = transpose(u) dot v
For C: <u,v> = conjugate(transpose(u)) dot v # physicists
For C: <u,v> = u dot conjugate(transpose(u)) # mathematicians

[Prasoon Shukla]

I had a 2 day break between exams so I began work on the CoordSystem class. Hopefully, I'll be able to make a WIP PR before 15th.

One problem I am facing is this: In this class, there is a lot of initialization going on. A part of that is generating non-commutative symbols to represent the basis elements. Since we planning to support n dimensional coordinate system for rectangular coordinates, I would need n basis elements.

So, I basically need to do initialization n times like so:

self.e1 = Symbol('e1') ... upto n times.

Obviously this is not feasible since we don't know n.

I can get a tuple of symbols using the symbols('e1: '+ str(n+1)) but then again is the problem of initializing the object properties.

The only way I can think of doing this is to use the exec statement.

sym_tuple = symbols('e1:' + str(n+1))

for i in range(1, n+1):

    cmd = 'self.e' + str(i) + ' = sym_tuple[ ' + str(i-1) + ']'

    exec cmd

But then again, I know that using exec is frowned upon in general. So, is there a better way to go about this?

[Aaron Meurer]

Try using numbered_symbols().

Aaron Meurer

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.

[Prasoon Shukla]

I have been reading the documentation of the mechanics module and looking at the code to understand how the vector module will fit in with `mechanics`. Most of the things are clear to me at this point.

Now I'd like to know a little more about using dyadics. I have read about dyadics from some web resources. A Dyadic class is already implemented in mechanics. Now, I have two options: Either to modify the Dyadic class present in mechanics to play nice with the vector module, or, to make a separate Dyadic class as a part of the vector module.

If dyadics can be used elsewhere in SymPy, then I think that the latter option would be better. If not, and if dyadics are being used only in mechanics, then we can just modify the existing Dyadic class.

So, which option would be better?

[Aaron Meurer]

I don't think dyadics are used elsewhere, but the real question is if they are useful outside of physics or not. If they are, then they should be defined outside of physics, and just used in mechanics. If there are certain aspects that only make sense in physics, then they should be defined outside of physics and subclassed in physics (I say physics instead of mechanics because this sort of thing actually applies to all classes in the physics module).  

Aaron Meurer

--

You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.

[Prasoon Shukla]

But that is exactly my question - where else can dyadics be used?
Perhaps someone from the physics community can comment?

Thanks