Welcome to Cadabra Q&A, where you can ask questions and receive answers from other members of the community.
0 votes

At the end of the page Programming in Cadabra, there is a function expand_nabla, which expands the covariant derivative in terms of partial derivative and contractions with the connection.

In that implementation the dummy index contracting the connection with the tensor is called p, i.e. is fixed. Since the dummy index is fixed, consecutive application of the covariant derivative cannot be expanded. For example, if we consider the code below,

foo := \nabla_{a}{ \nabla_{b}{ h^{c}_{d} } };

an error is raised since the expansion of the inner covariant derivative introduces the index p, and then the expansion of the second covariant derivative attempts to introduce another p index (already in use).


How could the definition of the expand_nabla function be generalised to modify automatically the name of the dummy index?

Moreover, assuming that we have a set of symbols with the property Indices assigned, Is it possible to preserve the type of index? i.e. to introduce a new dummy index using a symbol of the list.

in General questions by (12.5k points)


Was this ever resolved? Two questions please:

  1. Is there an agreed general expand_nabla() or expand_covariant_derivative() function that does not have the issue with expanding

    \nabla_{a}{ \nabla_{b}{ h^{c}_{d} } }
  2. Is there such a function that also includes covariant derivatives for the Tetrad/Vielbein formalism or only for the usual GR formalism?

Thank you! GPN

4 Answers

+1 vote
Best answer

The answer I'm posting (after two years) is based on an idea by @dominicprice (former student of Kasper, and developer of many aspects of cadabra).

The code

def expand_nabla(ex, *, index_name=None, index_set=None):
    if index_name is None:
        index_name = "expandnablaidx"
        Indices($expandnablaidx, expandnablaidx#$, $expandnablaindices$)
    index_counter = 1
    for nabla in ex[r'\nabla']:
        dindex = nabla.indices().__next__() 
        for arg in nabla.args():
            for index in arg.free_indices():
                t2:= @(arg);
                newidx = Ex(index_name + str(index_counter))
                index_counter += 1
                if index.parent_rel==sub:
                    t1:= -\Gamma^{@(newidx)}_{@(dindex) @(index)};
                    t2[index]:= _{@(newidx)};
                    t1:=  \Gamma^{@(index)}_{@(dindex) @(newidx)};
                    t2[index]:= ^{@(newidx)};
                ret += Ex(str(nabla.multiplier)) * t1 * t2
            nabla += ret
    distribute(ex, repeat=True)
    if index_set is not None:
        rename_dummies(ex, "expandnablaindices", index_set)
    return ex

How does it work?

The problem I found with the algorithm expand_nabla reported in the cadabra book https://cadabra.science/notebooks/ref_programming.html, was that when considering higher order derivatives, the name of the dummy index wouldn't change, resulting in a "repeated indices" error.

It is desirable that the algorithm could use one of the declared indices!

Therefore, we have to give a name to the set of indices:

{a, b, c, d, e, f, g, h, i, j, k ,l}::Indices(space, position=fixed).
{\mu,\nu,\lambda,\rho}::Indices(spacetime, position=fixed).

NOTE: there are two types of indices: space and spacetime.

Now, we declare an expression and expand the nabla operator (covariant derivative).

foo := \nabla_{c}{ h^{a}_{b} };
expand_nabla(_, index_set="space");

Since we have asked to expand with dummy indices of type space, the result is

$$\partial{c}{h^{a},{b}}+\Gamma^{a},{c d} h^{d},{b}-\Gamma^{d},{c b} h^{a},{d}$$

However, if we expand indices of type spacetime, we get

$$\partial{c}{h^{a}\,{b}}+\Gamma^{a}\,{c \mu} h^{\mu}\,{b}-\Gamma^{\mu}\,{c b} h^{a}\,{\mu}$$

Higher orders

Higher orders (I was able to calculate up to third order) can be computed with ease, for example:

foo := \nabla_{e}{\nabla_{d}{\nabla_{c}{ h^{a}_{b} }}};
#expand_nabla(_, index_set="vector");
expand_nabla(_, index_name=r'\lambda');

The index_name=r'\lambda' ask the algorithm to use the name \lambda to expand the dummy indices.

by (12.5k points)
selected by

Thank you. It works.

How to implement a zoom-sensitive algorithm
+1 vote

I faced the same problem in my time, so I ended up writing a function like this.

def select_index(used_indices):
    indeces = r'w v u s r q o n m l k j i h g f e d c b a'.split() 
    for uind in indeces:
        found = False
        for qind in used_indices:
            if qind == uind:
                found = True
        if not found:
            index = uind
    return Ex(index), used_indices

def one_nabla(ex, used_indices): 
    t3, used_indices = select_index(used_indices)
    free = dict()
    free['sub'] = set()
    free['up'] = set()
    for nabla in ex[r'\nabla']:
        dindex = nabla.indices().__next__() 
        for arg in nabla.args():             
            for index in arg.free_indices(): 
                if index.parent_rel==sub:
            for key in free.keys(): 
                for index in free[key]:
                    ind = Ex(index)
                    if key == 'sub': 
                        t1:= -\Gamma^{@[t3]}_{@(dindex) @[ind]};
                        t1:=  \Gamma^{@[ind]}_{@(dindex) @[t3]};
                    t2:= @[arg];
                    for term_index in arg.free_indices(): 
                        if str(term_index.ex()) == index:
                            if term_index.parent_rel==sub:
                                t2[term_index]:= _{@[t3]};
                                t2[term_index]:= ^{@[t3]};
                    ret += Ex(str(nabla.multiplier)) * t1 * t2 
            nabla += ret #
    return ex, used_indices

def nabla_calculation(ex, used_indices, count): 
    ex = ex.ex()
    for element in ex.top().terms(): 
        local_count = 0
        for nabla in element[r'\nabla']: 
            local_count += 1
        if local_count == 0: 
        elif local_count == 1: 
            new, used_indices = (one_nabla(element, used_indices))
            count -= 1
            if element.ex().top().name == r'\prod': 
                i = 0
                while i < 2:
                    for mult in element.ex().top().children():
                        local_count2 = 0
                        for nabla in mult[r'\nabla']:
                            local_count2 += 1
                        if local_count2 == 0:
                            new *= mult.ex()
                        elif local_count2 == 1: 
                            for nabla in mult[r'\nabla']:
                                nabla1, used_indices = one_nabla(nabla, used_indices)
                                new *= nabla.ex()
                            mult1, used_indices = nabla_calculation(mult, used_indices, local_count)
                            new *= mult1
                new *= Ex(str(element.multiplier))
                for nabla in element[r'\nabla']:
                    for arg1 in nabla.args():
                        arg2, used_indices = nabla_calculation(arg1, used_indices, count - 1)
                        index = nabla.indices().__next__() 
                        t := \nabla_{@(index)}{@[arg2]};
                    new = t
                    nabla1, used_indices = one_nabla(new, used_indices)
                    new = Ex(str(nabla.multiplier)) * nabla1
    return ex, used_indices

def expand_nabla(ex): 
    if ex.top().name == '\equals':
        for child in ex.top().children():
            for element in child.ex().top().terms():
                count = 0
                used_indices = set()
                for nabla in element.ex()[r'\nabla']:
                    count += 1
                if count == 0:
                    ret += element.ex()
        #new_ex += element.ex()
                    for n in element.ex():
                        for index in n.indices():
        #for nabla in element[r'\nabla']:
                    element1, used_indices = nabla_calculation(element, used_indices, count)
                    ret +=element1
        for element in ex.top().terms():
                count = 0
                used_indices = set()
                for nabla in element.ex()[r'\nabla']:
                    count += 1
                if count == 0:
                    for n in element.ex():
                        for index in n.indices():
                element1, used_indices = nabla_calculation(element, used_indices, count)
    return ex

To use it, you need to call the main function expand_nabla(ex). It is suitable for large expressions, does not duplicate already occupied indexes, and also knows how to work with equalities.

by (970 points)

This routine throws a RuntimeError: Free indices in different terms in a sum do not match.

Thank you for the comment @Arina. Sorry for my (super) later reply.

+1 vote

There's a very simple workaround for this.

{\mu, \nu, \lambda, \kappa, \rho, \sigma#}::Indices(space, position=independent).


A{#}::Depends(\partial{#}, \nabla{#}).

ex1:= A_\mu.
ex2:= B_{\mu \nu} = \nabla_{\mu}{ A_\nu }.


ex3:= C_{\mu \nu \rho} = \nabla_{\mu}{B_{\nu \rho }}.
substitute(ex3, ex2)


You essentially have to rename dummy indices after each use of expand_nabla.

by (740 points)

Thank you for the comment @jsem.

This only works if you have one nabla in the summand. If there is more than one, they all get the same dummy index.

What do you mean? If you need more nablas, you can act with one, rename dummies, act again with the next nabla, etc. It's not an elegant solution, hence only a workaround.

I mean that if you have something like \nabla_{\rho}(\nabla_{\nu}(\xi_{\mu})),

this function from Programming in Cadabra expands it with two identical dummies.

That is true. That was the original question. My code produces a workaround that seems to do the job.

Sorry, I was inattentive, that's what you offer the solution for. Unfortunately, if we are talking about expressions with ~100 summands, this approach is very problematic((

Indeed. It's not a real solution, just a workaround for simple cases.

Your code didn't seem to run at all for me, so I don't know of a better solution unfortunately.

Sorry for that, you can email me, maybe we can figure out what the problem is. arshtm@gmail.com

The example in the original question no longer throws an error, in Cadabra2 version 4+. But it still has a problem with getting indexes correctly. I think I solved it though - see answer below.

0 votes

OK, I went and fixed the example from the programming page and here is my suggestion. It actually takes the idea from @jsem to add a rename_dummies() after each operator expansion, and tries to do so systematically.

Please do review!

def expand_covariant_derivative(ex, operator=r'\nabla'):
    Expand the covariant derivative operator into partial
    coordinate derivatives and Christoffel symbols.

    :operator: which operator is used for the covariant derivative, default `\nabla`

    **Warning** whatever you choose for the `operator`, you must also add a definition
    of that operator as a derivative, otherwise indexes will not be handled correctly.
    For example:

    #                \nabla{#}::Derivative.

    This function deals with the covariant derivative with respect to generalized
    coordinates, not with respect to flat coordinates as used in the Vierbein (Tetrad)
    def expand_one_appearance(ex):
        Finds the first appearance of `operator` and expands it with a single
        partial derivative and Christoffel symbols
        for oper in ex[operator]:
            dindex = oper.indices().__next__() 
            for arg in oper.args():             
                for index in arg.free_indices():
                    t2:= @(arg).
                    if index.parent_rel==sub:
                        t1:= -\Gamma^{\gamma}_{@(dindex) @(index)}.
                        t2[index]:= _{\gamma}.
                        t1:=  \Gamma^{@(index)}_{@(dindex) \gamma}.
                        t2[index]:= ^{\gamma}.
                    ret += Ex(str(oper.multiplier)) * t1 * t2
                oper += ret
        return ex

    # loop that deals with the covariant-derivative operators, one at a time
    # First we expand one of them
    # Then we rename dummies, to avoid name conflicts
        while len(list(ex[operator])) > 0:
            ex_temp := @(ex).
            ex = expand_one_appearance(ex_temp)

    except RuntimeError as reason:

    return ex
by (1.2k points)