Feeds:
Posts
Comments

## SAGE tip: More about Differential Forms

In a previous post we discuss the definition of the coordinated patch on a manifold, how to define differential forms, wedge them or calculate their exterior derivative… even simplify’em.

This time a zero form will be defined and a list of forms will be created… So, let’s begin!

# Define a 0-form

Once created the coordinated patch and the differential forms algebra

sage: reset()
sage: var('t,x,y,z')
sage: U = CoordinatePatch((t,x,y,z))
sage: Omega = DifferentialForms(U)

A 0-form is defined as an element of $\Omega^0(U)$, but the value of the 0-form is given inside the declaration command,

sage: A = DifferentialForm(Omega, 0, exp(x*y))

I tried addition, multiplication, wedge product and exterior differentiation on 0-forms and they worked!

Of course you can combine them with forms of different degrees.

# New method of defining a form

I wrote to Joris this morning… but before he was able to answer, from the documentation of the differential form package.

When one calls the generators of the differential form,

sage: Omega.gen(1)
dx

and the result is a differential form… Thus, one can assign a form as follow,

sage: A = sin(x)* Omega.gen(2)
sage: B = cos(y) * Omega.gen(0)
sage: C = sin(z) *Omega.gen(1)
sage: D = cos(y) * Omega.gen(2)

And this forms can be wedged, differentiated, et cetera. ðŸ™‚

# List of forms

Finally, after discovering the above behavior I tried the following, a list of differential forms ;-), for example,

sage: pro = matrix([[A, B], [C, D]])
sage: for i in range(2):
...       for j in range(2):
...           show(pro[i,j].diff())

returns
$cos(x)dx\wedge dz$
$sin(y)dx\wedge dy$
$-cos(z)dy\wedge dz$
$-sin(y)dy\wedge dz$

This implies that somehow one can manage a series of forms by using list properties… I expect to go deeper on this subject in the future! ðŸ˜€

Enjoy people!!!

Dox

Advertisements

Read Full Post »

## SAGE tip: GRmodule. Day 07.

This code is supposed to be (if some one does the work in the future) located in sage.tensor.differential_form_element.

The code presented below is a slight modification of Joris code for differential forms manipulation on SAGE.

## Needed modules

from sage.symbolic.ring import SymbolicRing, SR
from sage.rings.ring_element import RingElement
from sage.algebras.algebra_element import AlgebraElement
from sage.rings.integer import Integer
from sage.combinat.permutation import Permutation

The advantage of using this is that, tensors defined here are an Algebra element, not just a python object as in the previous code.

The sage.combinat.permutation won’t be used (yet), but could be useful if tensor symmetries are defined.

#### TensorFormatter

class TensorFormatter:
r"""
This class contains all the functionality to print a tensor in a
graphically pleasing way.  This class is called by the _latex_ and
_repr_ methods of the Tensor class.
"""
def __init__(self, space):
r"""
Construct a tensor formatter.  See
TensorFormatter for more information.
"""
self._space = space

def repr(self, comp, fun):
r"""
String representation of a primitive tensor, i.e. a function
times a tensor product of d's of the coordinate functions.

INPUT:

- comp -- a subscript of a differential form.

- fun -- the component function of this form.

EXAMPLES::

sage: from sage.tensor.tensor_element import TensorFormatter
sage: x, y, z = var('x, y, z')
sage: U = CoordinatePatch((x, y, z))
sage: D = TensorFormatter(U)
sage: D.repr((0, 1), z^3)
'z^3*dx@dy'

"""

str = "@".join( \
[('d%s' % self._space.coordinate(c).__repr__()) for c in comp])

if fun == 1 and len(comp) > 0:
# We have a non-trivial form whose component function is 1,
# so we just return the formatted form part and ignore the 1.
return str
else:
funstr = fun._repr_()

if not self._is_atomic(funstr):
funstr = '(' + funstr + ')'

if len(str) > 0:
return funstr + "*" + str
else:
return funstr

def latex(self, comp, fun):
r"""
Latex representation of a primitive differential form, i.e. a function
times a tensor product of d's of the coordinate functions.

INPUT:

- comp -- a subscript of a differential form.

- fun -- the component function of this form.

EXAMPLES::

sage: from sage.tensor.tensor_element import TensorFormatter
sage: x, y, z = var('x, y, z')
sage: U = CoordinatePatch((x, y, z))
sage: D = TensorFormatter(U)
sage: D.latex((0, 1), z^3)
'z^{3} d x \otimes d y'

"""

from sage.misc.latex import latex

str = " \otimes ".join( \
[('d %s' % latex(self._space.coordinate(c))) for c in comp])

if fun == 1 and len(comp) > 0:
return str
else:
funstr = latex(fun)

if not self._is_atomic(funstr):
funstr = '(' + funstr + ')'

return funstr + " " + str

def _is_atomic(self, str):
r"""
Helper function to check whether a given string expression
is atomic.

EXAMPLES::

sage: x, y, z = var('x, y, z')
sage: U = CoordinatePatch((x, y, z))
sage: from sage.tensor.tensor_element import TensorFormatter
sage: D = TensorFormatter(U)
sage: D._is_atomic('a + b')
False
sage: D._is_atomic('(a + b)')
True
"""
level = 0
for n, c in enumerate(str):
if c == '(':
level += 1
elif c == ')':
level -= 1

if c == '+' or c == '-':
if level == 0 and n > 0:
return False
return True

The only I’ve changed here is “DifferentialForm” by “Tensor” and “\wedge” by “\otimes”

The above code allows to write the tensor product in a basis, $dx^1\otimes\cdots\otimes dx^n$. The chosen symbol for denoting the tensor product was @.

## Tensor Class

This code is incomplete due to:

• I’ve not defined the TensorsAlgebra, which should be done in parallel.
• There are a lot of attributes not presented in this class.
• class Tensor(AlgebraElement):
r"""
Tensor class.
"""

def __init__(self, parent, degree, fun = None):
r"""
Construct a tensor.

INPUT:

- parent -- Parent algebra of tensors.

- degree -- Degree of the tensor.

- fun (default: None) -- Initialize this differential form with the given function.  If the degree is not zero, this argument is silently ignored.

EXAMPLES::

sage: x, y, z = var('x, y, z')
sage: F = Tensors(); F
Algebra of tensors in the variables x, y, z
sage: f = Tensor(F, 0, sin(z)); f
sin(z)

"""

from sage.tensor.tensorss import Tensors
if not isinstance(parent, Tensors):
raise TypeError, "Parent not an algebra of tensors."

RingElement.__init__(self, parent)

self._degree = degree
self._components = {}

if degree == 0 and fun is not None:
self.__setitem__([], fun)

def __getitem__(self, subscript):
r"""
Return a given component of the tensor.

INPUT:

- subscript: subscript of the component.  Must be an integer
or a list of integers.

EXAMPLES::

sage: x, y, z = var('x, y, z')
sage: F = Tensors(); F
Algebra of tensors in the variables x, y, z
sage: f = Tensor(F, 0, sin(x*y)); f
sin(x*y)
sage: f[()]
sin(x*y)
"""

if isinstance(subscript, (Integer, int)):
subscript = (subscript, )
else:
subscript = tuple(subscript)

dim = self.parent().base_space().dim()
if any([s >= dim for s in subscript]):
raise ValueError, "Index out of bounds."

if len(subscript) != self._degree:
raise TypeError, "%s is not a subscript of degree %s" %\
(subscript, self._degree)

"""sign, subscript = sort_subscript(subscript)"""

if subscript in self._components:
return sign*self._components[subscript]
else:
return 0

def __setitem__(self, subscript, fun):
r"""
Modify a given component of the tensor.

INPUT:

- subscript: subscript of the component.  Must be an integer or a list of integers.

EXAMPLES::

sage: F = Tensors(); F
Algebra of tensors in the variables x, y, z
sage: f = Tensor(F, 2)
sage: f[1, 2] = x; f
x*dy@dz
"""

if isinstance(subscript, (Integer, int)):
subscript = (subscript, )
else:
subscript = tuple(subscript)

dim = self.parent().base_space().dim()
if any([s >= dim for s in subscript]):
raise ValueError, "Index out of bounds."

if len(subscript) != self._degree:
raise TypeError, "%s is not a subscript of degree %s" %\
(subscript, self._degree)

"""sign, subscript = sort_subscript(subscript)"""
self._components[subscript] = SR(fun)

Ok, so again I’ve changed “DifferentialForm(s)” by “Tensor(s)”, drop the permutation of indices (’cause tensors do not need to be neither symmetric nor anti-symmetric.

I’ll keep working with this code… It’s all by now. Oh! I’ll post next week the rest of the code based in Sergey’s GRPy.

Enjoy.

Dox

Read Full Post »

From time to time when one works with operators, such as in Quantum mechanics, something like an exponential of the operator appears (this is also the case in many areas of Mathematics like group theory of differential geometry). This exponentiation of a matrix should be understood as the series expansion of the exponential.

SAGE knows how to do this exponentiation,

sage: reset()
sage: var('a,b,c', domain=RR)
sage: A = a*I*matrix([[0,1],[1,0]])
sage: B = b*I*matrix([[0,-I],[I,0]])
sage: C = c*I*matrix([[1,0],[0,-1]])
sage: A.exp()
[1/2*(e^(2*I*a) + 1)*e^(-I*a) 1/2*(e^(2*I*a) - 1)*e^(-I*a)]
[1/2*(e^(2*I*a) - 1)*e^(-I*a) 1/2*(e^(2*I*a) + 1)*e^(-I*a)]
sage: B.exp()
[   1/2*(e^(2*I*b) + 1)*e^(-I*b) -1/2*(I*e^(2*I*b) - I)*e^(-I*b)]
[ 1/2*(I*e^(2*I*b) - I)*e^(-I*b)    1/2*(e^(2*I*b) + 1)*e^(-I*b)]
sage: C.exp()
[ e^(I*c)        0]
[       0 e^(-I*c)] 

The only problem here is that, even when the relations are simple in this example, I’ve not found a trivial’ way of simplifying the matrix elements of the exponentiation, not even with the procedure post in here. I didn’t try with the rewrite package

Comments are welcome!

Enjoy.

Dox

Read Full Post »

## SAGE tip: GRmodule. Day 06.

I’ve just updated the SAGE worksheet which uses the definitions described in the previous posts.

• There are some explanations in text format
• The code has been hidden… because is long.
• Moreover… I’ve discover something really amazing! Joris Vankerschaver‘s code of the differential form package in SAGE. Thus, I could use some ideas from Joris’ code to improve GRmodule. Nice, Isn’t it?

Let’s hope I could so something nice this weekend!

Don’t forget check the worksheet, and post some comment for feedback! ðŸ˜‰

Enjoy!

Dox

Read Full Post »

Two days ago I was trying to expand in series a lot of functions… so I ask myself, Could it be done in SAGE? It should be possible… but, How? ðŸ˜›

So I ask on google group.

#### Solution by Andrzej Chrzeszczyk

sage: var('r');
sage: f=2*r/sinh(2*r)
sage: f.taylor(r,0,5)
14/45*r^4 - 2/3*r^2 + 1
sage: maxima(f).powerseries('r',0)
-4*r*'sum((2^(2*i2-1)-1)*2^(2*i2-1)*bern(2*i2)*r^(2*i2-1)/
(2*i2)!,i2,0,inf)

This solution uses a power series expansion from maxima… really nice feature! Isn’t it?
Ah… and this expansion is around $r=0$.

If one would like the asymptotic expansion $r\to\infty$,

sage: maxima(f).powerseries('r',infinity)
-4*r*'sum((2^(2*i3-1)-1)*2^(2*i3-1)*bern(2*i3)*r^(2*i3-1)/
(2*i3)!,i3,0,inf)

However, note that this expansion coincides with the previous one, i.e., it’s the function itself. It couldn’t be that perfect. ðŸ˜‰

#### Solution by Francois Maltey

Use the Taylor command of SAGE,

• Around zero
sage: taylor (2*x/sinh(2*x), x, 0, 10)
-292/13365*x^10 + 254/4725*x^8 - 124/945*x^6 + 14/45*x^4 - 2/3*x^2 + 1
• Around infinity… a trick! change $x\mapsto 1/x$ and expand around zero ðŸ™‚
sage: taylor (2*1/x/((exp(2/x)-exp(-2/x))/2), x, 0, 12)
4*e^(-10/x)/x + 4*e^(-6/x)/x + 4*e^(-2/x)/x
• Thank you guys!

Enjoy!

Dox

Read Full Post »

## SAGE tip: GRmodule. Day 05

Hi everyone!

This time the Christoffel connection will be defined.

# The code

As usual, here is the code:

class Christoffel(Tensor):
'''The class to represent Christoffel Symbols of the second kind. Please
note that while it inherits from Tensor, Christoffel symbols are
NOT tensors'''

def __init__(self,metr,symbol='C',rank=(1,2),sh=(1,-1,-1)):

# The metric
self.g_down = metr

# Since we have a metric we do indeed have a coordinate system
self.rep  = self.g_down.rep

self.g_up = metr.inverse

# Please note that this call will trigger a call to allocate in
# the Tensor class, but the allocate will actually be the allocate
# defined below
super(Christoffel,self).__init__(symbol,rank,sh,coords=metr.coord)

def allocate(self,rank):
Tensor.allocate(self,rank)
# Now that we have allocated things, time to actually compute things
for i in xrange(self.dim):
for k in xrange(self.dim):
for l in xrange(self.dim):
sum = 0
for m in xrange(self.dim):
term1 = diff(self.g_down[m,k],self.g_down.coords[l])
term2 = diff(self.g_down[m,l],self.g_down.coords[k])
term3 = diff(self.g_down[k,l],self.g_down.coords[m])

tr = self.g_up[i,m] * (term1+term2-term3)

sum += tr
res = sum/2
self.components[i,k,l] = res
self.getNonZero()

This code is almost a copy of Sergey’s one, except for the use of xrange instead of np.arange, and the fact that I’ve dropped the minus signs denoting the shape of the tensors.

## Sage implementation

This time I won’t present a Python file, but a SAGE file, GRmodule.

Enjoy!

Dox

Read Full Post »

## SAGE tip: Working with Differential Forms.

A few weeks ago I realize there is a package for working with differential forms in SAGE. So, I’ve tried to take advantage of that!

This is what I’ve learnt so far.

# Define your manifold

If you work with differential forms you should know (probably better than me), that you define forms in a coordinate patch… rarely globally.

I’d like to start with something like this!

sage: reset()
sage: # The variables
sage: var('theta,phi,r')
sage: var("bps", latex_name=r"\bar{\psi}")
sage: var("bph", latex_name=r"\bar{\phi}")
sage: var("bth", latex_name=r"\bar{\theta}")
sage: # The coordinate system
sage: coords = [r,theta, phi, bps, bth, bph]
sage: U = CoordinatePatch((r,theta, phi, bps, bth, bph))
sage: Omega = DifferentialForms(U)

From a previous post a weird declaration of variables has been made.

But don’t lose yourself into the code!!!… the important here is:

• Seven variables have been declared,
• These are converted into a list,
• and also are the coordinated of the Coordinate Patch.
• Finally, $\Omega$ is the algebra of differential forms. Just like $\Omega^\bullet(M)$, where $M$ is the manifold (in this case the patch).
• ## Define the differential forms

Differential form have to be defined… like this,

sage: A = DifferentialForm(Omega, 1)
sage: A[1] = -1/4
sage: A[2] = 1/4*(-cos(theta) - sin(theta))
sage: A[3] = cos(bps)*sin(bth)*sin(bph)
sage: A[4] = sin(bps)*cos(bth)*sin(bph)
sage: A[5] = sin(bps)*sin(bth)*cos(bph)

The first line defines $A$ as an element of $\Omega^1(M)$.

The rest of the lines are for setting the components of $A$. Remember that the coordinates run from 0 to 6 in this example! ðŸ˜› And, of course, they are ordered strictly as we declared them, i.e.,

$A = -\frac{1}{4}d\theta -\frac{1}{4}(\cos(\theta) + \sin(\theta))d\phi + \cos(\bar{\psi})\sin(\bar{\theta})\sin(\bar{\phi})d\bar{\psi} +\cdots$

## Exterior Differentiation

Once the form has been declared, one might differentiate it by using the diff command,

sage: A.diff()

returns

$(\frac{1}{4} \, \sin\left(\theta\right) - \frac{1}{4} \, \cos\left(\theta\right)) d \theta \wedge d \phi$

## Showing the forms

It’s very useful to see the formulas in a nice written way… this is one of the features I love the more from SAGE, to see your form, use

 sage: show(A)

## Wedge Product

Of course one of the most important operations when working with forms is the wedge product. For using this, try

sage: C = A.wedge(A.diff())
sage: show(C)

it vanishes. One might try to define another form,

sage: B = DifferentialForm(Omega,2)
sage: B[0,1] = sin(bth)

Note that $B\in \Omega^2(M)$, and has only one non-vanishing component.

Thus, $A\wedge B$ is calculated,

sage: A.wedge(B)

One can also define multiple wedge product,

sage: A.wedge(B.diff()).wedge(B)

Needless to say one can multiply a form by a function, or number.

## Simplifying a form

Forms have not implemented the simplify_full attribute, but their components, which are functions, do. So, after a very complicated calculation one might try to implement a long routine of simplification, say,

sage: D = A.wedge(B.diff()).wedge(B)
for i in xrange(Omega.ngens()):
for j in xrange(i+1, Omega.ngens()):
for k in xrange(j+1, Omega.ngens()):
for l in xrange(k+1, Omega.ngens()):
for m in xrange(l+1, Omega.ngens()):
D[i,j,k,l,m] = D[i,j,k,l,m].simplify_full()


In the above, Omega.ngens()` returns the dimension of the patch. Additionally, each index runs from the value of the previous plus one, in order to avoid repetition (due to anti-symmetric property) or the zeros values.

if someone find a better way, doesn’t hesitate in post a comment!!! ðŸ™‚

Well, I think this is much of it! I’ll keep you posted, in case I learn more about ot.

Enjoy life! and Happy Valentine’s day! ðŸ˜‰

Dox

Read Full Post »

Older Posts »