Vectors EM framework - sympy/sympy GitHub Wiki

Note: See this google doc for updated changes to the API. Link. This is just bare bones API to get started with coding - more will be added as we code. Please comment on the ML if you want any changes. Currently Sachin and Prasoon have edit access.

For the moment the page is separated in a proposal by Sachin and a proposal by Prasoon.

Sachin

The inspiration for the RefFrame/ CoordSys framework has been taken from this article.

Important points to be noted -

"An observational frame (such as an inertial frame or non-inertial frame of reference) is a physical concept related to state of motion."

"An observer in an observational frame of reference can choose to employ any coordinate system (Cartesian, polar, curvilinear, generalized, …) to describe observations made from that frame of reference. A change in the choice of this coordinate system does not change an observer's state of motion, and so does not entail a change in the observer's observational frame of reference."

"The idea of a reference frame is really quite different from that of a coordinate system. Frames differ just when they define different spaces (sets of rest points) or times (sets of simultaneous events). So the ideas of a space, a time, of rest and simultaneity, go inextricably together with that of frame. However, a mere shift of origin, or a purely spatial rotation of space coordinates results in a new coordinate system. So frames correspond at best to classes of coordinate systems."

Let us look at the steps involved in the solution of a sample problem -

The global Symbol for time and the global frame are defined. For a given RefFrame, the sense of 'time' is SAME for all coordinate systems defined in it. Hence, all coordinate systems defined within a global RefFrame use the same symbol for time.

>>> t = Symbol('t')
>>> frame = RefFrame('frame', dim = 3, time_var = t)

Comment from Gilbert: I think it would be easier to have the time symbol defined under the ReferenceFrame class (not instances); e.g. ReferenceFrame.t == symbols('t').

Define a rectangular coordinate system as an instance of CoordSysRect, specifying variables you define. If the coordinate variables are not defined by the user, they default to symbols of type 'CoordinateSymbol(systemname_e1)'.

>>> A = Point('A')
>>> a_x, a_y, a_z, b_x, b_y, b_z = CoordinateSymbols('a_x, a_y, a_z, b_x, b_y, b_z')
>>> system_a = frame.attach_system(name = 'a', 'rectangular', origin = A, vars = [a_x, a_y, a_z])

The variables a_x, a_y and a_z are automatically linked with the system they are used in. Hence a_x.system() returns system_a. Plus, system_a.varlist() returns a tuple of all the variables defined in it Because a coordinate system needs a frame of reference, only a frame can initialize a coordinate system. The type of the system is defined while attaching via a string. By default, its rectangular. The attach function automatically initialises an object of the corresponding class.

Comment from Gilbert: Is the introduction of Point here really necessary, vs just introducing a translation vector on coordinate system creation?

[I would prefer a translation vector. I just introduced Point because Stefan used points to denote origins of the systems. However, using just translation vectors is a better option. Stefan agrees? But again, if Prasoon is thinking of allotting coordinates to every point in a system wrt an origin, its existence may be necessary. He should comment on this.]

To get the ith basis vector(starting from 1), use systemname[i]

>>> system_a[1] == Vector(system_a, [1,0,0])
True

Define system_b, which has an origin at B (which has a certain position vector with respect to Point A, the origin of system_a). The definition of origins helps to determine translation of one system with respect to another. (As of now, the API of Points is same as that of mechanics module)

>>> B = Point('B')
>>> B.set_pos(A, system_a[1])
#Define axis_AB. Taken simple for now, could be any vector defined in A
>>> axis_AB = system_a.basis(2)
>>> theta_AB = Symbol('theta_AB')

Orient system_b in a certain way with respect to system_A. Like the current mechanics framework, the orientation of one system with respect to another can be changed. Rotational orientation can be changed with 'orient' function, for translation, just redefine the position of origin of one with respect to the other.

'orient' function will change the behaviour of functions like vector.express(system). In a way, the, the systems are mutable, but only as far as orientations and location of origin is concerned. All other factors remain same.

>>> system_b = frame.attach_system(name = 'b', frame, origin = B, vars = [b_x, b_y, b_z])
>>> system_b.orient(system_a, translate = [], rotate = ['Axis', [theta_AB, axis_AB]])

Define system C as required

>>> C = Point('C')
>>> C.set_pos(A, 2 * system_a[2] + system_b.basis[1])
>>> system_c = frame.attach_system(name = 'c', 'spherical', frame', vars = [r, phi, theta]))

CoordSysRect, CoordSysSph are classes of coordinate systems. However, for reasons already explained, a user will not initialise a coordinate system, he will do it via a frame. Hence,

>>> system_c == frame['c']
True
>>> type(system_c)
...CoordSysSph

Comment from Gilbert: Does a user do anything with a coordinate system besides access its basis vectors?

[CoordSys will have functions for accessing basis vectors, the coordinate variables, and orienting it with respect to some other system attached to the SAME frame]

TODO: Rotate C

Define the scalar fields. CoordinateSymbols of more than one coordinate systems can be used in definition of one scalar field.

>>> sfield1 = ScalarField(a_x + 2 * a_y)
>>> sfield2 = ScalarField(b_x**2)

Comment from Gilbert: Can you list all the possible operations someone would do with a scalar field? e.g. re-expression, gradient, ...?

Operations on scalar fields-

  1. Re-expression in other systems

  2. Gradient

  3. Laplacian

  4. Testing equality (supposing the same scalar field has been expressed in different ways)

  5. Substitution (finding value at a particular point in a particular system)

A vector field will be an instance of (a subclass of?) Vector

>>> vfield = 3 * r * system_c[1] + r * system_c[2]
>>> sfield3 = sfield1 + sfield2

ScalarFields and Vectors can be expressed in any system

>>> sfield3.express(system_a)
1 + a_x + 2*a_y + ((a_x + 1)*cos(theta) + a_y*sin(theta))**2
>>> gradient(sfield3, system_b)
(2*b_x + sin(theta) + cos(theta)) * system_b[1] + (-sin(theta) + cos(theta)) * system_b[2]

Like the sympy.physics.mechanics module, sfield3.express(system_a) creates another instance. Doing sfield3.express(system_a)==sfield3 returns True. This will be done similar to the way '==' has been overloaded in the current mechanics module for Vector. Instead of a structural comparison, it will focus on the VALUEs of the fields in a common system (that's how equality of fields should be tested). Functions like 'subs', 'atoms' will work on the class.

Comment from Gilbert: equality comparisons between scalar fields seem like they would be impossible to perform without their own class (i.e., just representing scalar fields as a SymPy Expr wouldn't work). Comment from Gilbert: Actually, thinking about this, if we did CoordSys.express(some_scalar_field), that would work if some_scalar_field was an Expr, as the CoordSys would know what its coordinate variables are, and those of the other coordinate systems that have been defined relative to it.

[Yup, system.express(sfield) would work. Its that or the other way round (I prefered that way to be consistent with the current vector.express() method) The only problem remaining is the equality of fields. A user could check by -

  1. Find the system to which any one of the CoordinateSymbols in the expression belongs

  2. Express both fields in same system and then equate in the usual way

We could have a function like 'equate_fields' if we are intent on avoiding a new class.]

Comment from Gilbert: Right now, in sympy.physics.mechanics, ReferenceFrames allow the user to set custom names for basis vectors, and custom LaTeX strings. which should be kept. See: http://docs.sympy.org/0.7.2/modules/physics/mechanics/advanced.html#referenceframe

TODO

Just add an 'x' to whatever you consider finished (like this [x]).

  • Try to format the wikipage somewhat better. Instead of comments (hard to read) just use text. Think of this as if it is writing a tutorial for you future module meant for people that have never seen it. Do not answer questions below, just make the answer obvious in the tutorial above.
  • The origins of system_a, system_b and system_c are at different points (A, B and C respectively). In the example above all origins seem to be the same point.
  • The python variables a_x, a_y, a_z and b_x, b_y, b_z are never created. CoordSysRect('system_a', frame, vars = [a_x, a_y, a_z]) will just raise NameError: name 'a_x' is not defined
  • You use system_b.orient(...) does this mean that coordinate systems are mutable objects?
  • What is sfield3.express(system_a) doing? Is it modifying the sfield3 object?
  • Rotate system_c as Stefan suggests
  • sfield3.express(system_a) == sfield3 being True means that == will not do structural comparison. This is not the way sympy works, mainly because it will break all tree-manipulation routines.
  • Will subs, atoms, free_symbols, has, etc work on these?
  • Why is the ScalarField class necessary? I do not see any use for it. All the information is already present in the Symbol instance that you use.

Prasoon

Alright, so here's my explanation of things. Note that I didn't read Sachin's explanation prior to writing this so that my thought process didn't get influenced.

First, I shall describe a few things about the Vector module in brief. This really I considered necessary so that everyone knows clearly exactly why is a part needed and how is it to be used. Then, I have tried to provide the mock session.

The Vector calculus module in brief:

So first, I'll explain the basics once. Then using that, we will try to give a mock session of the stress test.

The basics:

First, I will explain using vectors that do not depend on time and are used in just one 'global' frame (in the sense that all calculations will be done wrt to the same origin of coordinates).

Now, to create a vector, we require two things:

  • A coordinate system

A coordinate system will be represented with a CoordSys object. Please note that a CoordSys object is not the same as a reference frame. (I shall explain this in a bit). For the time being, the CoordSys object provides a container to hold the basis elements (whether constant or space dependent), the h parameters (that are used to provide general formula for grad, div, curl etc), several methods that help in conversion between the different coordinate systems and methods that provide functionality that is a there due to usage of a particular coordinate system. So as you can see, this object is really concerned with providing functionality for a particular type of a coordinate system. For example, in the rectangular system, the basis is constant everywhere - that is a feature of that particular coordinate system. The CoordSys object represents this feature. So in the end, the CoordSys object is really used to just represent the features of a certain coordinate system and to provide a mathematical representation of that system. Nothing more. This means that this object won't be translated, rotated etc. This object just provides a way to represent points in space. That is it.

  • The Vector itself

A Vector object holds the components of the vector we want to represent. So, in effect, using these two classes (Vector and CoordSys), we have separated the components of a vector and the basis elements of it. Now, the vector class will contain methods that operate on the vector. Note that these methods will work in conjunction with the methods of the CoordSys class.


Now, let me talk about the RefFrame class a bit since it seems to have caused quite a bit of confusion. This class is supposed to be a sort of a wrapper over the CoordSys object. What is a reference frame? There can be many definitions. According to me, a reference frame is just different set of coordinate axes to measure quantitites in. These coordinate axes may have been translated, rotated and/or given kinematic variables wrt to a frame of our choosing. Now, if we assess from this definition, then we could just as well have added some more functionality in the CoordSys class and we could have gotten the same results without needing a different RefFrame class. But, I decided to make different classes. One of the attributes of the RefFrame class will be a CoordSys object. Why I decided to make a separate implementation is more of matter of design. Traditionally, a vector module should just have to deal with vectors that change with coordinates, not time. Therefore, again going with the modular nature of SymPy, I decided that it was best to implement a separate RefFrame class that takes care of time related kinematic variables, translation, rotation etc.

Now, let's talk about functionality of the RefFrame class and how things will work.

First thing is about having different reference frames in the same coordinate system. For example, let us say we want two frames of reference for measuring a certain vector quantity. We also want both of those frames to measure this vector in Cartesian coordinates. So, we will define two RefFrame objects, provide these objects with details about their translational/rotational parameters, pass the coordinate system we are to use (here Cartesian) and then measure the required quantity. Using this example, it is easy to see the distinction between a coordinate system and a reference frame in that coordinate system. The coordinate system provides the methods to represent and measure quantities in space whereas the reference frame provides a frame of reference in any given coordinate system and allows us to associate kinematic variable to that frame.

Now let us see how things will work internally. When we define a vector in a reference frame, the vector will be stored just as if it was defined in the global frame. The reason for this is quite obvious - since the vector is defined in a certain frame, obviously the components of the vector will be in that frame too. The usage of the frame computationally comes in when we try to express this vector in another frame. Let me take an example for this:

Let's have two frames A and B. A measures quantities in Cartesian coordinates while B does it in spherical coordinates. A has its origin at Point_A and B has its origin at Point_B, each wrt to the global coordinates. Orientations of all frames are the same. Now, we define a vector in A. The vector has components (a, b, c) in A. So the Vector object will just store the components as (a, b, c). The ref_frame attribute of Vector will contain the RefFrame for A. So far so good. Note that we haven't done any conversions (that is computing how the vector will appear in another frame) yet.

Now, if we want to express the vector in the reference frame B, then we will need to perform some computations. Now, if this vector is a free vector, then we will just convert it to spherical coordinates and return the new vector as being in frame B. If it is a vector field, then we will need to first compute the vector field in reference frame B and then we can convert it to spherical coordinates. This will be the general procedure of things internally.

I hope now the distiction between a coordinate system and a reference frame is clear. Now, we can have a mock session of the problem

>>> from sympy import *
>>> from sympy.vector import *
>>> x, y, z = symbols('x y z')
>>> r, phi, theta = symbols('r phi theta')

>>> c_rect = CoordSysRect(name='c_rect', dim=3)
>>> c_sph = CoordSysSph(name='c_sph')

>>> point_A = Point(name='point_A', coordinates=(1, 1, 1), coord_system='rect')
>>> point_B = Point(name='point_B', coordinates=(2, 2, 2), coord_system='rect')
>>> point_C = Point(name='point_C', coordinates=(-1, -1, -1), coord_system='rect')

# Let us keep RefFrame as mutable
>>> frame_A = RefFrame(name='frame_A', coord_sys=c_rect, vars=(x, y, z))
>>> frame_A.set_position(point_A)

>>> frame_B = RefFrame(name='frame_B', coord_sys=c_rect, vars=(x, y, z))
>>> frame_B.set_position(point_B)

>>> frame_C = RefFrame(name='frame_C', coord_sys=c_sph, vars=(r, phi, theta))
>>> frame_C.set_position(point_C)

Now we rotate the frames A and B. An axis can always be defined by using two points. So first we need two axis points.

>>> a_1 = Point(coordinates=(0, 0, 0), coord_system='rect')
>>> a_2 = Point(coordinates=(3, 4, 5), coord_system='rect')

# Define clock wise
>>> theta_AB = Symbol('theta_AB')
>>> frame_A.rotate(a_1, a_2, theta_AB, '+')
>>> frame_B.rotate(a_1, a_2, theta_AB, '-')

>>> alpha, beta, gamma = symbols('alpha beta gamma')

# I do not get the point mentioning Euler angles. I'll add it once clear on that.

>>> sc_A = ScalarField(frame=frame_A, x+y+z)
>>> sc_B = ScalarField(frame=frame_B, x)
>>> vc_C = Vector(coords=(r**2, r*sin(phi), cos(theta)), frame=frame_C)

We need a reference frame to show the sum of scalar fields in.

Since no frame has been supplied, hence, return results in the global frame. Also since both sc_A and sc_B have the same coordinate systems underneath, hence we do not need to supply a coordinate system either.

>>> G = grad(sc_A + sc_B)

# If we want it in any other frame, say in A, then
>>> G_A = grad(sc_A + sc_B, frame=frame_A)
>>> vc_new = custom_op(G, 'dot', 'nabla', oper_on=vc_C, frame=frame_A)

# To express to any other coordinate system or frame
# Express vc_C in in frame B
>>> vc_C_in_frameB = vc_C.convert_to(frame=frame_B)
# To express a scalar field in another frame (or coordinate system)
>>> sc_A_in_frameB = sc_A.convert_to(frame=frameB)

TODO

  • The way you use Coordinate Systems and References are very confusing. It is a nice idea to abstract away the difference between rectangular, spherical, etc systems, however if this is a hardcoded abstraction why at all should the user bother to create instances of CoordSys which are only created to be used in RefFrame? Just create appropriate subclasses (RectangularRefFrame, SpericalRefFrame, etc).

Comment from Sachin : I think this entire discussion on frames and systems is getting confusing. At this rate, reaching a consensus seems difficult. Prasoon, could you provide a link to the source which you are using as a reference for these definitions like I have? It will make things much clearer and get us all on the same page. I am saying this because my source has very different ideas from yours, so we need to take a decision pronto.

Comment from Prasoon : I didn't take the definition from a source exactly; I was just specifying the things that I think are intuitive and natural for defining a frame of reference. But let us not dwell on the definitions. And okay, I think we can proceed with subclassing the RefFrame class for each particular reference frames.

  • I see no reason to specify the same variables on the creation of your RefFrames. This leads to the extremely confusing definition of scalar fields: sc_A = ScalarField(frame=frame_A, x+y+z); sc_B = ScalarField(frame=frame_B, x) where the meaning of x, y, z is completely different in the two cases.

Comment from Sachin : We don't need to specify the frame/system during definition of ScalarField. The definition may contain coordinate variables from different systems. To express in one particular frame, the user may then use sfield.express().

Comment from Prasoon : Well, I suppose this can get confusing. So we can perhaps use something like we used in mechanics before : like R.x for x coordinate in R.

Comment from Sachin : Yup. But sfield = R_x + R1_x should be allowed. So no need to give the system. But ya, the variable will have information about the system it is attached to.

  • Why you require RefFrame to be mutable?

Comment from Sachin : Considering the RefFrame needs to have parameters of its motion with respect to some other frame, it will be mutable in the sense that we can change its velocity etc. if needed.

Comment from Prasoon : It is as Sachin says.

  • On numerous occasions you use strings like 'rect'? Why? You already have coordinate systems that are defined to be rectangular - if the whole point of CoordSys was to abstract away the type of coordinate system that is used, why are you then using some ad-hoc string interface?

Comment from Prasoon : The only place where I have used the string is while defining points - the coordinates of which can be provided in any coordinate system. That is why we had those strings there.

  • You require a global reference frame. Why? The whole point of the module is to express vector calculus naturally without depending on some global reference frame in which the answers will always be simpler or in which it will always be easier to work. For instance you specify positions of reference frames not wrt each other but wrt the global frame.

Comment from Sachin : I have stressed this before, we don't need a global frame. The user can define one on his own to be the 'global'. But mathematically speaking, there is no need. The only thing global, as Gilbert pointed out, must be the time variable.

Comment from Prasoon : Let's see what a global frame is. It is nothing special - it is just the lack of a proper defined frame. So from user's point of view, s/he doesn't need to do anything more to define a global frame. So, there is no global frame - it is just an abstraction.

Another Proposal

Because of the unnatural parametrization of SO(n) (the rotation group) in usual 3D mechanics the proposal below works only for 2D and 3D. If nD is necessary, the user will be expected to create his own class BlahReferenceFrame. I believe Gilbert will have a better solution (see the line defining cart_B).

Prepare symbols for coordinates:

x_a, y_a, z_a, x_b, y_b, z_b, rho, theta, phi = symbols(...)

Prepare symbols for base vectors:

e_x_a, e_y_a, e_z_a, e_x_b, e_y_b, e_z_b, e_rho, e_theta, e_phi = symbols(...)

Prepare symbols for relative positions and speeds:

delta, v, alpha, omega = symbols(...)

And a symbol for time dependence:

t = Symbol(...)

Define a Cartesian reference frame:

cart_A = CarthesianReferenceFrame(x_a, y_a, z_a, e_x_a, e_y_a, e_z_a, t) 
#                                 \___________/  \_________________/  |
#                                   coords          base vectors     time

Define another Carthesian system related to the first one:

cart_B = CarthesianReferenceFrame(x_a, y_a, z_a, e_x_a, e_y_a, e_z_a, t, cart_A, delta, 0, 0, alpha, 0, 0, v, 0, 0, omega, 0, 0)
#                                 \___________/  \_________________/     \____/  \_________/  \_________/ \___ ___/  \_________/
#                                   coords          base vectors          wrt     translation  rotation    velo-  angular velocity
#                                                                        which                              city
#                                                                      coordinate
#                                                                        system

Define a Spherical reference system with the same origin and orientation as cart_A:

sph = SphericalReferenceFrame(rho, theta, phi, e_rho, e_theta, e_phi, t, cart_A)

Maybe define a world:

world = World(sph, cart_A, cart_B)

Define a scalar field and take its gradient:

scalar = x_a+y_b+theta
grad_scalar = world.gradient(scalar)
print grad_scalar
# e_x_a+e_y_b+1/r*e_theta

A dot product (at some point we might want to have a container Dot, but it does not seem necessary for the moment):

print world.dot(e_x_b, e_rho)
#cos(alpha+omega*t-theta)

And nonsense detection:

print world.is_vector(x_a)
# False
world.gradient(e_x_a)
# TypeError
nonsense = x_a + e_x_a # No error when executing Add
world.<any method>(nonsense)
# TypeError
⚠️ **GitHub.com Fallback** ⚠️