New to sympy: Implement Spherical Coordinates in Vector module
1,032 views
Skip to first unread message
Lukas Zorich
Nov 19, 2014, 5:05:35 PM11/19/14
to sy...@googlegroups.com
Hi everyone,
I am Lukas from Chile, pursuing a degree in Computer Science and with interests in mathematics. I would like to start contributing to sympy. Reading through the mailing list, I saw this post that says that the Spherical Coordinates System isn't implemented yet, so I would like to start from there. I wanted to know who can I ask to give me further information about this specific functionality, or things I need to know before implementing this (not about the theory, but about the concrete implementation of the spherical system in the vector module).
Thanks :)!
Lukas
Jason Moore
Nov 19, 2014, 5:17:24 PM11/19/14
to sy...@googlegroups.com
Lukas,
Sachin Joglekar developed the vector package this past summer. I was his GSoC mentor on the project. You can check out this implemented of the Cartesian coordinate system:https://github.com/sympy/sympy/blob/master/sympy/vector/coordsysrect.py
--
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.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/d6a7e710-7d96-452d-8a7b-6518c0226225%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
Lukas Zorich
Nov 19, 2014, 10:34:23 PM11/19/14
to sy...@googlegroups.com
Thanks Jason!
Alan Bromborsky
Nov 20, 2014, 8:45:50 AM11/20/14
to sy...@googlegroups.com
You might wish to be more general and
implement vector analysis as a function of a metric tensor. See
the following link -
http://en.wikipedia.org/wiki/Curvilinear_coordinates#Vector_operations
http://en.wikipedia.org/wiki/Curvilinear_coordinates#Vector_operations
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/86d19d32-44b7-4530-b24a-3742c3276099%40googlegroups.com.
Lukas Zorich
Nov 21, 2014, 11:53:13 AM11/21/14
to sy...@googlegroups.com
El jueves, 20 de noviembre de 2014 10:45:50 UTC-3, brombo escribió:You might wish to be more general and implement vector analysis as a function of a metric tensor. See the following link -
http://en.wikipedia.org/wiki/Curvilinear_coordinates#Vector_operations
Can you explain me a bit more? you mean that instead of calculating dot product or cross product in a specific coordinate system, I calculate that using the metric tensor as it appeared in the wikipedia article?
Alan Bromborsky
Nov 21, 2014, 12:31:45 PM11/21/14
to sy...@googlegroups.com
Yes. It would be a good idea to
restrict yourself to diagonal metric tensors (orthogonal
coordinate system) then inverting the metric tensor is simple and
the Jacobian matrix is also diagonal. If you have a metric tensor
you would use the Christoffel symbols to calculate the derivatives
of the basis vectors. For a spherical coordinate system
(r,theta,phi) the diagonal of metric tensor would be
(1,r**2,r**2*sin(theta)**2).
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/7105867e-e256-44d3-a2a8-8fed73f11798%40googlegroups.com.
Aaron Meurer
Nov 21, 2014, 2:58:14 PM11/21/14
to sy...@googlegroups.com
Note that a challenge with spherical coordinates is picking a
convention (see
https://en.wikipedia.org/wiki/Spherical_coordinate_system#Conventions).
We either need to pick one and document it clearly, and always use it,
or find a clean way to support multiple conventions.
Aaron Meurer
> https://groups.google.com/d/msgid/sympy/CAP7f1AjM26X9RRu%2B5g0_rndJWyYZ%2B6t7%3D6zvkg6XEQ24jjbv4Q%40mail.gmail.com.
convention (see
https://en.wikipedia.org/wiki/Spherical_coordinate_system#Conventions).
We either need to pick one and document it clearly, and always use it,
or find a clean way to support multiple conventions.
Aaron Meurer
Francesco Bonazzi
Nov 25, 2014, 3:58:49 AM11/25/14
to sy...@googlegroups.com
On Wednesday, November 19, 2014 11:17:24 PM UTC+1, Jason Moore wrote:
Sachin Joglekar developed the vector package this past summer. I was his GSoC mentor on the project. You can check out this implemented of the Cartesian coordinate system:
https://github.com/sympy/sympy/blob/master/sympy/vector/coordsysrect.py
In sympy.diffgeom coordinate systems are defined as instances of the CoordSystem class, once more coordinate systems are defined, the transformation rules are set using the connect_to method:
In [1]: from sympy.diffgeom import *
In [2]: r, theta, phi = symbols('r theta phi', positive=True)
In [4]: M = Manifold('M', 3)
In [5]: P = Patch('P', M)
In [7]: rect = CoordSystem('rect', P, ['x', 'y', 'z'])
In [8]: spherical = CoordSystem('spherical', P, ['r', 'theta', 'phi'])
In [9]: spherical.connect_to(rect, [r, theta, phi], [r*sin(phi)*cos(theta), r*sin(phi)*sin(theta), r*cos(phi)], inverse=False)
I put inverse=False in the last line, because SymPy's solvers are unable to reverse the spherical-to-rect transformation.
Manifold and Patch are some boilerplate, Manifold specifies the space you're considering, while Patch is a connected subspace of the manifold. Patch is used because on some manifolds there does not exist a homeomorphism to a coordinate system (e.g. on the sphere).
I think that some algorithmic code could be shared between the vector and the diffgeom modules.
Aaron Meurer
Nov 27, 2014, 1:17:49 AM11/27/14
to sy...@googlegroups.com
solve() compute things (similar to how sympy.stats works).
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.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit
Francesco Bonazzi
Dec 6, 2014, 9:58:52 AM12/6/14
to sy...@googlegroups.com
I published a gist on how to get the gradient in any coordinate system, given the backward transformation to rectangular coordinates.
I was also trying to write the code to get the divergence, but I ran into problems with sympy's simplify( ) function.
https://gist.github.com/Upabjojr/34f75ce115f1cf5a6afc
This gist defines a function get_grad( ... ) which returns the gradient, it requires as an input a lambda function (backwards transformation).
Example:
Notice that the function from_spherical is simply the easiest way to define the spherical coordinates.
I get the following output:
SymPy doesn't know that in theta's range of values its sin(theta) is always positive, this is the reason why it does not get simplified.
I have only tried the spherical coordinates, but this should work with any coordinate transformation in three dimensions. It should not be too complicated to generalize this to any dimension.
That's why I would like coordinate systems to be instances of a class, and not classes themselves. It is really uncomfortable to define a new class for every new coordinate system.
I was also trying to write the code to get the divergence, but I ran into problems with sympy's simplify( ) function.
https://gist.github.com/Upabjojr/34f75ce115f1cf5a6afc
This gist defines a function get_grad( ... ) which returns the gradient, it requires as an input a lambda function (backwards transformation).
Example:
from_spherical = lambda r, theta, phi: (r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta))
grad = get_grad(from_spherical)
Notice that the function from_spherical is simply the easiest way to define the spherical coordinates.
I get the following output:
⎡ ∂_θ │sin(θ)│⋅∂ᵩ⎤
⎢∂ᵣ, ───, ───────────⎥
⎢ r 2 ⎥
⎣ r⋅sin (θ) ⎦
SymPy doesn't know that in theta's range of values its sin(theta) is always positive, this is the reason why it does not get simplified.
I have only tried the spherical coordinates, but this should work with any coordinate transformation in three dimensions. It should not be too complicated to generalize this to any dimension.
That's why I would like coordinate systems to be instances of a class, and not classes themselves. It is really uncomfortable to define a new class for every new coordinate system.
Alan Bromborsky
Dec 6, 2014, 3:15:42 PM12/6/14
to sy...@googlegroups.com
This is what I do for sqrt of real
expression x (I want sqrt(abs(x)) without the abs() function in
the result unless it is unavoidable)
def square_root_of_expr(expr):
if expr.is_number:
if expr > 0:
return(sqrt(expr))
else:
return(sqrt(-expr))
else:
expr = trigsimp(expr)
(coef, pow_lst) = sqf_list(expr)
if coef != S(1):
if coef.is_number:
coef = square_root_of_expr(coef)
else:
coef = sqrt(abs(coef))
for p in pow_lst:
(f, n) = p
if n % 2 != 0:
return(sqrt(abs(expr)))
else:
coef *= f ** (n / 2)
return coef
This was done specifically for dealing with cylindrical and spherical coordinate systems where you have to normalize the basis vectors
and express the derivatives of the basis vectors in terms of the normalized basis vectors.
If you have a metric for spherical coordinates -
diagonal matrix(1,r**2,r**2*sin(th)**2)
and do not normalize the basis vectors you never have to take the square root of a function.
def square_root_of_expr(expr):
if expr.is_number:
if expr > 0:
return(sqrt(expr))
else:
return(sqrt(-expr))
else:
expr = trigsimp(expr)
(coef, pow_lst) = sqf_list(expr)
if coef != S(1):
if coef.is_number:
coef = square_root_of_expr(coef)
else:
coef = sqrt(abs(coef))
for p in pow_lst:
(f, n) = p
if n % 2 != 0:
return(sqrt(abs(expr)))
else:
coef *= f ** (n / 2)
return coef
This was done specifically for dealing with cylindrical and spherical coordinate systems where you have to normalize the basis vectors
and express the derivatives of the basis vectors in terms of the normalized basis vectors.
If you have a metric for spherical coordinates -
diagonal matrix(1,r**2,r**2*sin(th)**2)
and do not normalize the basis vectors you never have to take the square root of a function.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/14ce54ae-897d-41bd-9a35-36e9a0b11822%40googlegroups.com.
Alan Bromborsky
Dec 7, 2014, 12:13:58 PM12/7/14
to sy...@googlegroups.com
Here is commented version of my
"square_root_of_expr" function. The functions determines if the
input expression is a perfect square and if it is returns the
positive square root of the perfect square. The other cases for
the input expression are now documented in the code. I assume
the expression is real and if
the expression is a number and negative you want the square root of the absolute value of the expression. This is done for the case of normalizing the basis
vectors in a Minkowski space. That is if the square of a basis vector is negative then the square of the normalized basis vector should be -1. This only works for
basis vectors whose squares are numbers and not functions of the coordinates since in general sympy cannot determine if a function is positive or negative over
a given range.
def square_root_of_expr(expr):
"""
If expression is product of even powers then every power is divided
by two and the product is returned. If some terms in product are
not even powers the sqrt of the absolute value of the expression is
returned. If the expression is a number the sqrt of the absolute
value of the number is returned.
the expression is a number and negative you want the square root of the absolute value of the expression. This is done for the case of normalizing the basis
vectors in a Minkowski space. That is if the square of a basis vector is negative then the square of the normalized basis vector should be -1. This only works for
basis vectors whose squares are numbers and not functions of the coordinates since in general sympy cannot determine if a function is positive or negative over
a given range.
def square_root_of_expr(expr):
"""
If expression is product of even powers then every power is divided
by two and the product is returned. If some terms in product are
not even powers the sqrt of the absolute value of the expression is
returned. If the expression is a number the sqrt of the absolute
value of the number is returned.
"""
if expr.is_number:
if expr > 0:
return(sqrt(expr))
else:
return(sqrt(-expr))
else:
expr = trigsimp(expr)
(coef, pow_lst) = sqf_list(expr)
if coef != S(1):
if coef.is_number:
coef = square_root_of_expr(coef)
else:
coef = sqrt(abs(coef)) # Product coefficient not
a number
for p in pow_lst:
(f, n) = p
if n % 2 != 0:
return(sqrt(abs(expr))) # Product not all even
powers
else:
coef *= f ** (n / 2) # Positive sqrt of the square of an expression
return coef
else:
coef *= f ** (n / 2) # Positive sqrt of the square of an expression
return coef
On 12/06/2014 09:58 AM, Francesco Bonazzi wrote:
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/14ce54ae-897d-41bd-9a35-36e9a0b11822%40googlegroups.com.
Francesco Bonazzi
Dec 8, 2014, 7:43:29 AM12/8/14
to sy...@googlegroups.com
Interesting function, can be really a useful workaround.
Anyways, I think that the most polite way is to use assumptions, such as:
Anyways, there's a bug while simplifying expressions made up of diffgeom's coord functions, example:
The variables x, y should behave exactly as SymPy's symbols in the expressions they appear in, but their type is BaseScalarField, they contain information in the arguments of which coordinate system they belong, on which patch, and on which manifold, and especially which coordinate they are. Is there any way to instruct simplify( ) that it should not enter the arguments of BaseScalarField in its algorithm?
Anyways, I think that the most polite way is to use assumptions, such as:
>>> refine(sqrt(abs(sin(t)**2)), Q.positive(sin(t)))
sin(t)
Anyways, there's a bug while simplifying expressions made up of diffgeom's coord functions, example:
In [1]: from sympy import *
In [2]: from sympy.diffgeom import *
In [3]: M = Manifold('M', 2)
In [4]: P = Patch('P', M)
In [5]: rect = CoordSystem('rect', P, ['x', 'y'])
In [6]: x, y = rect.coord_functions()
In [7]: x*y
Out[7]: x*y
In [8]: simplify(x*y)
/usr/local/lib/python2.7/dist-packages/sympy/core/function.pyc in count_ops(expr, visual)
2250 while args:
2251 a = args.pop()
-> 2252 if a.is_Rational:
2253 #-1/3 = NEG + DIV
2254 if a is not S.One:
AttributeError: 'int' object has no attribute 'is_Rational'
The variables x, y should behave exactly as SymPy's symbols in the expressions they appear in, but their type is BaseScalarField, they contain information in the arguments of which coordinate system they belong, on which patch, and on which manifold, and especially which coordinate they are. Is there any way to instruct simplify( ) that it should not enter the arguments of BaseScalarField in its algorithm?
Reply all
Reply to author
Forward
0 new messages