GSoC 2017 Application Sagar Bharadwaj : Implementation of multiple types of coordinate systems for vectors - sympy/sympy GitHub Wiki

Personal Information

Heading Details
Name Sagar Bharadwaj KS
University National Institute of Technology Karnataka, Surathkal
Major Computer Science
E-Mail [email protected]
GitHub SagarB-97
Github Page https://sagarb-97.github.io/
Timezone IST (UTC + 5:30)

Background

I am a second year undergraduate studying at National Institute of Technology Karnataka, Surathkal. I have always had a keen interest in Mathematics and Algebra since higher primary school. I have completed three undergraduate Mathematics courses including Engineering mathematics 1, Engineering Mathematics 2 and Concrete Mathematics. Other courses of relevance I have taken include Data Structures and Algorithms course and Design and Analysis of Algorithms.

Programming Skills

  • I primarily work on an Ubuntu 16.04 (LTS) system. I prefer Atom over other editors because it's open source and has a plethora of plugins and packages I can choose from.

  • I have been programming for about 5 to 6 years now. Other languages I frequently use for my projects include Java, C and C++. I have involved myself in a lot of programming projects. My Github Page has details about all the projects I've taken up.

  • I am pretty comfortable with Python as I've been using it for over a year now. My favourite feature of python is List comprehension. It allows us to create lists with a syntax similar to what a Mathematician would use to create Sets.

  • My favourite features of SymPy are the series and rewrite functions.

    • The series function defined in sympy.core.expr.Expr helps in expanding an expression f(x) around x = x0.
      Example :

      >>> e = sin(x)*cos(y)     
      >>> e.series(x)     
      x*cos(y) - x**3*cos(y)/6 + x**5*cos(y)/120 + O(x**6)     
      >>> e.series(x,y)     
      sin(y)*cos(y) + (x - y)*cos(y)**2 - (x - y)**2*sin(y)*cos(y)/2 - (x - y)**3*cos(y)**2/6 +
      (x - y)**4*sin(y)*cos(y)/24 + (x - y)**5*cos(y)**2/120 + O((x - y)**6, (x, y))     
      >>> e.series(y,x)     
      sin(x)*cos(x) - (-x + y)*sin(x)**2 - (-x + y)**2*sin(x)*cos(x)/2 + (-x + y)**3*sin(x)**2/6 + 
      (-x + y)**4*sin(x)*cos(x)/24 - (-x + y)**5*sin(x)**2/120 + O((-x + y)**6, (y, x))     
    • The rewrite function in sympy.core.basic.Basic which is implemented using _eval_rewrite_as_ functions spread across many classes is one of the features unique to SymPy. Example :

      >>> e = fibonacci(n)     
      >>> e.rewrite(GoldenRatio)     
      (GoldenRatio**n - (-GoldenRatio)**(-n))/(-1 + 2*GoldenRatio)     
      >>> e = sin(x)*cos(x)     
      >>> e.rewrite(tan)     
      2*(-tan(x/2)**2 + 1)*tan(x/2)/(tan(x/2)**2 + 1)**2
  • I am comfortable with git and GitHub as I use them to coordinate with my project mates in most of my projects.

The Project

Implementation of multiple coordinate systems in SymPy has been in discussion for a long time. I went through the discussions in Pull Requests #9937, #10074, #10109 and also some discussions that happened on the mailing list. I referred to a paper by Victor I. Piercey and Alan Bromborsky's document on the topic. Also, the existing sympy.diffgeom module provided me some implementation ideas. After some extensive reading on the topic, here is my proposal of the idea.

Motivation

Cartesian coordinate system is the most widely used coordinate system because of its simplicity and ease of use. However, many seemingly complex problems with some identifiable obvious symmetry are reduced to trivial cases by merely changing the coordinate system. Therefore, implementing multiple coordinate systems in a Computer Algebra System is inevitable. The importance of this project for the evolution of SymPy into a fully featured CAS is what motivates me to take it up.

Theory

This proposal mainly focuses on the implementation of Orthogonal Curvilinear Coordinate systems in SymPy. Since the Cartesian coordinate system is in place and almost the complete implementation is available in sympy.vector module, it would be best to base the implementation of other coordinate systems on Cartesian system.

Consider a coordinate system W with (w1, w2, w3) being the orthogonal curvilinear coordinates or Base scalars.

(x, y, z) refers to Cartesian coordinates. Let the coordinates in the two systems be related by the following 3 functions :

  • x = f (w1, w2, w3)
  • y = g (w1, w2, w3)
  • z = h (w1, w2, w3)

A good implementation always strives to generalise expressions and encompass a wide range of situations. In order to generalise expressions for different coordinate systems, we need a finite number of coefficients that can wholly represent a coordinate system. After having done this, obtaining specific expressions from generalised expressions would merely involve substituting the respective values for these coefficients.

Lamé Parameters are one set of such coefficients that can completely represent the differential behaviour of a coordinate system.

The Jacobian Matrix J can be defined as :

Jacobian Matrix

Now the Metric Tensor M for the coordinate system is given by :
Metric Tensor

If the given system is an Orthogonal Curvilinear coordinate system, the metric tensor is a Diagonal Matrix. The Lamé Parameters for such a system can be easily derived from the Metric Tensor by taking the square root of the diagonal elements.

That is, Lamé Parameters h1, h2, h3 are defined as :
h1
h2
h3

Also, the base vectors (which are mutually orthogonal and unit vectors) for this new coordinate system can be defined as :
Base Vectors

Where,
R

Now, position vectors and vector fields can be defined in terms of the new base vectors in the new coordinate system.

Expressions for many differential operations in Coordinate systems can be generalised using Lamé Parameters.

  • Del Operator
    Del Operator

  • Gradient of a Scalar function
    Gradient

  • Divergence of a vector field
    Let
    Vector
    Then,
    Divergence

  • Curl of Vector field
    Curl

  • Laplacian of Scalar Function
    Laplacian

Other essential functions such as conversion of base scalars from one coordinate system to Cartesian and vice versa can be achieved by simply substituting values in f, g, h functions (described above) or their counterparts relating w1, w2, w3 to x, y, z.
Similarly vector transformations can be achieved by using the mutually orthogonal base vector relationships discussed above.

The best way to achieve inter-conversion between different coordinate systems is to internally convert them to Cartesian.

Implementation

Implementing a single Python class that can represent all the coordinate systems is the best way forward.
The class would take in the type of the coordinate system (str) or a Lambda object representing functions f, g, h described above in case a custom coordinate system needs to be designed. (It would also take in the other optional arguments including location, parent, vector_names, variable_names)

  • Example
    coordSystem = CoordinateSystem('W', type = 'cylindrical') or
    coordSystem = CoordinateSystem('W', type = Lambda((r, theta, z), (r*cos(theta), r*sin(theta), z))

The type of the coordinate system can be a string for built in coordinate systems. They would again be internally mapped onto a Lambda object representing f, g and h functions. For such built in coordinate systems, fastest performance will be ensured if Jacobian Matrix, Lamé Parameters and base vector relations are also hard coded into the class for immediate retrieval.

In case of custom coordinate system, the Lamé Parameters, Jacobian Matrix and base vector relations can be derived using the relations described above.
The advantage of generalising expressions (described above) and using a single Python class to represent all coordinate systems is evident here. Now, there is no difference with respect to implementation between custom coordinate systems and inbuilt coordinate systems.

Expressing coordinates in one system in terms of another

Achieving scalar transformation(passive) would involve a simple substitution as implemented in coord_tuple_transform_to and connect_to in sympy.diffgeom.diffgeom.CoordSystem.
Similarly, vector transformation(passive) would require using the base vector relations mentioned previously.
As mentioned in theory section, converting to and fro from Cartesian would be the best method for inter-conversion.

Differential Geometry Functions

The Del operator's definition needs to be changed to a generic one as mentioned in the theory section. The implementation of gradient, curl and divergence in sympy.vector.functions conveniently use delop. So, they might not need major changes as dot product and cross product definitions will be accordingly changed to suit all coordinate systems. (These definitions also are similar to Cartesian ones already implemented in sympy.vector module).

The Del class must be rewritten to return the Del operator itself (The current implementation only has differential geometry functions and does not return the actual del operator). This would enable one to overload the * and ^ operators and calculate the differential function results without using the gradient, curl, divergence or dot and cross functions. For example :

  • divergence = (del operator) * (Vector Filed)
  • curl = (del operator) ^ (Vector Field)

Scalar Potential

The vector fields would now be described by using the new generalised base vectors. The infinitesimals and the vector field must be represented in the same coordinate system after which similar integration would follow.

Orienters

Orienters would require some reimplementation as they have been primarily written for Cartesian coordinate system.
I would prefer a simple redefinition of Rotation Matrix that would be general for all coordinate systems. Such a Rotation Matrix can be calculated from Rodrigues' rotation formula. The Rodrigues formula takes the axis of rotation and angle of rotation as input and gives a general rotation matrix irrespective of coordinate system. The formula is :
Rodrigues
where,
R -> Rotation Matrix
K -> Cross product vector for a unit vector denoting rotation of axis

Further details of Implementation in different modules is presented in detail in the Timeline.

Clarifications

  • The coordinate systems that can be 'built in' include the orthogonal curvilinear coordinate systems mentioned here:

    • Spherical polar coordinates
    • Cylindrical polar coordinates
    • Parabolic cylindrical coordinates
    • Paraboloidal coordinates
    • Elliptic cylindrical coordinates
    • Prolate spheroidal coordinates
    • Oblate spheroidal coordinates
    • Ellipsoidal coordinates
    • Bipolar coordinates
    • Toroidal coordinates
    • Conical coordinates
  • The given base vectors are Covariant basis. Contravariant Basis can also be represented by dividing covariant basis vectors by Lamé Parameters as they both are in same direction in case of orthogonal curvilinear coordinate systems.

Timeline

I would be able to set aside around 40 hours for SymPy per week on an average.

  • I'll devote around 4-5 hours on weekdays and 7-8 hours on weekends.
  • My college reopens at the end of July. But that would not limit the number of hours I can work as the beginning of the semester would not involve a hectic schedule.
  • The timeline might seem a little non uniform. This is because I have allotted myself sparse work for the period I might have other commitments.

May 4 - May 30 (Community Bonding Period)

I would spend this time studying the already existing sympy.vector module thoroughly and comprehending the changes required to be implemented. I would also spend some time on the sympy.diffgeom module which has implementation of some differential geometry functions. I would interact with mentors and get concrete details of expected work to be completed.

In short : Comprehend thoroughly the existing code and interact with mentors.

May 31 - June 6

Rewrite the class CoordSysCartesian in sympy.vector.coordsysrect. The name of the class would be changed. The class would be rewritten such that it forms the main Python class that represents all Orthogonal Curvilinear coordinate systems (described above). This would involve considerable alterations.
For Example : The user will no longer prefer to refer to the base scalars as coordsystem.x, coordsystem.y and coordsystem.z as these are cartesian conventions. Rather, provision will have to be made to refer to them using the vector names passed. (Currently they are used only in str printing).
eg. coordsystem.r, coordsystem.theta, coordsystem.phi.
Same argument would hold true for base vectors.

The class will also contain hard coded dictionaries with built in coordinate system names (string) mapped to f, g and h functions described above. There will also be dictionaries mapping them to Lame parameters. Other orient_new_* functions use sympy.vector.Orienter class and can remain unchanged.

In Short : Change the basic definition of CoordSysCartesian so that it can represent any coordinate system.

June 7 - June 11

Write functions to calculate Lame Parameters, Jacobian Matrix and Metric Tensor in the new rewritten class. These functions would be used to set these parameters for a custom defined orthogonal curvilinear coordinate system.

In Short : Take care of custom defined coordinate system.

June 12 - June 14

Write tests and documentation for the altered class. Some basic tests would also be written while coding the class.

June 15 - June 22

Work on sympy.vector.deloperator.Del class. This class would need almost a complete redefinition because all differential operations associated with Coordinate systems is done here. This would involve rewriting the dot, cross and gradient functions to suit multiple coordinate systems. A new function to calculate laplacian of a scalar function should also be written.

Alternative : The Del class can be rewritten so that it only returns the del operator. The calculation of differential geometry functions like gradient, curl, divergence and laplacian can then be done only in sympy.vector.functions. This would also allow calculation of these functions using * and ^ operators and del objects directly by the user.
The alternative approach seems better when compared to just rewriting the existing class architecture.

In Short : Take care of differential geometry functions in multiple coordinate systems.

June 23 - June 25

Write tests and documentation for the differential geometry functions coded earlier. Some basic tests would be written while coding also.

June 26 - June 30

Phase 1 Evaluation time.
I would also spend some time writing more tests for differential functions because of the enormous number of cases that should be considered. (As there are many possible coordinate systems and various differential functions for each)

July 1 - July 9

Work on sympy.vector.orienter. This file would require a complete redefinition to accommodate multiple coordinate systems. The classes use rot_axis1, rot_axis2 and rot_axis3 functions in sympy.matrices.dense which were written with Cartesian system in mind. Rodrigues method to calculate rotation matrix must be implemented and all the Orienter class's sub classes must be updated.

In Short : Work on Orienters for different coordinate systems.

July 10 - July 12

Write tests and documentation for the new orienter functionality.

July 13 - July 20

Work on sympy.vector.point.Point class. Functions such as position_wrt and express_coordiantes need to be redefined.
New parameters indicating the output coordinate system needs to be included. That is, if the user wants the relative position vector in terms of a coordinate system other than Cartesian.
The algorithm where position is calculated based on the Point tree will be retained. Additional code should be added that converts the calculated Position vector (which is in Cartesian) to the desired system.

In Short : Redefine the Point class

July 21 - July 24

Write tests and documentation for the redefined Point class. Also write any tests and documentation that were not written earlier due to lack of time and prepare for Evaluation 2.

July 25 - July 28

Second Evaluation period.
I might be busy in some courses selection and registration for new College semester. So I'll spend this time correcting and completing any previously incomplete work.

July 29 - August 5

Work on sympy.vector.basisdependent. This file requires many changes.
For example : Integration of a vector expressed in Spherical coordinate system can be done using Cartesian coordinates as Integration variables. Such cases might require the use of Jacobian Matrix defined earlier. The same holds true for differentiation as well. Details of all changes required in basisdependent.py will be collected during community bonding period.

In Short : Work on basisindependent.py

August 6 - August 10

Some functionality defined earlier might have given wrong results as they depend on functionality defined in basisdependent.py. All such wrong results are noted and retested in this period. Any changes that should be reflected back to the previously completed work is also done.

More tests are written to exclusively test redefined functionalities in basisdependent.py.

August 11 - August 17

Work on sympy.vector.vector and sympy.vector.scalar. All changes required in these files will be determined in the course of writing other modules. Such changes must be implemented in this period.

August 18 - August 21

Complete any unfinished work.
Write well formatted documentation and tests for the Multiple coordinate system for vectors.

August 22 - August 29
Final evaluation Period

Pull Request Submissions

In the period mentioned under each division, I'll be sending a new commit to my pull request (or creating a new PR, whatever is appropriate). Some commits might fail many tests and might not be merged because of inter dependencies. But all such anomalies will be cleared eventually.

  • After redefinition of CoordSysCartesian class before writing tests and docs for it.

    • Date : June 12 - June 13
  • After redefinition of Del class and taking care of differential geometry functions.

    • Date : June 23 - June 24
  • After defining new Orienter classes.

    • Date : July 10 - July 11
  • After redefining Point class

    • Date : July 21 - July 23
  • After redefining basisdependent and taking care of previous fails that depended on functions in basis dependent.

    • Date : August 7 - August 10
  • Final submission after correcting all previous fails

    • Date : August 21 - August 22

Contributions to SymPy

Pull Requests

  • (Merged) str printing of logic expressions changed #12108
  • (Merged) fibonacci rewritten in terms of GoldenRatio #12167
  • (Merged) Evaluation of a fibonacci limit #12167
  • (Merged) Corrections to str set printing #12112

Concerned Issues

  • (Fixed) Needed to change str printing of logic expressions #11435
  • (Fixed) fibonacci limit evaluation #10382
  • (Fixed str part) str printing of sets #10672

References

⚠️ **GitHub.com Fallback** ⚠️