The calculation of reciprocal basis for non-orthogonal basis in GAlgebra deviate from standard definition and caused a bug when calculating inverse metric.
By definition,
$$e_{i} \cdot e^{j}=\delta_{i}^{j}$$where $\delta_{k}^{j}$ is the kronecker delta.
But GAlgebra, intentionally makes it:
$$e_{i} \cdot e^{j}={E_n}^{2}\delta_{i}^{j}$$where
$$ E_{n} = e_{1}\wedge \dots \wedge e_{n} $$In the docstring of Ga.build_reciprocal_basis(), it explicitly states:
For non-orthogonal basis $e^{j}$ is not normalized and must be divided by ${E_n}^{2}$ (
self.e_sq) in any relevant calculations.
The implementation of Ga.build_reciprocal_basis() claims to use formula 4.94 from GA4P by Doran and Lasenby:
But it looks like this:
# Take all (n-1)-blades
duals = list(self.blades_lst[-(self.n + 1):-1])
# After reverse, the j-th of them is exactly e_{1}^...e_{j-1}^e_{j+1}^...^e_{n}
duals.reverse()
sgn = 1
self.r_basis = []
for dual in duals:
dual_base_rep = self.blade_to_base_rep(dual)
# {E_n}^{-1} = \frac{E_n}{{E_n}^{2}}
# r_basis_j = sgn * duals[j] * E_n so it's not normalized, missing a factor of {E_n}^{-2}
r_basis_j = collect(expand(self.base_to_blade_rep(self.mul(sgn * dual_base_rep, self.e_obj))), self.blades_lst)
self.r_basis.append(r_basis_j)
# sgn = (-1)**{j-1}
sgn = -sgn
Let's verify this logic step by step using Kerr-Debney Metric:
In [1]:
from sympy import *
In [2]:
import sys
from galgebra.printer import Format, GaLatexPrinter
Format()
from galgebra.ga import Ga
from galgebra.mv import ONE, ZERO, HALF
W = symbols("W")
c = symbols("c")
g4coords = (u, x, y, z) = symbols("u x y z")
g = [
[0, 0, -exp(-z), 0],
[0, HALF * u ** 2 * exp(4 * z), 0, 0],
[-exp(-z), 0, 12 * exp(-2 * z), u * exp(-z)],
[0, 0, u * exp(-z), HALF * u ** 2],
]
g4 = Ga("e_u e_x e_y e_z", g=g, coords=g4coords, norm=False) # Create g4
(e_u, e_x, e_y, e_z) = g4.mv()
In [3]:
from IPython.display import display, Math
def show(x):
display(Math(GaLatexPrinter.latex(x)))
In [29]:
g4.g
Out[29]:
In [4]:
g4.basis
Out[4]:
In [5]:
for blade in g4.blades_lst:
show(blade)
Take all (n-1)-blades:
In [6]:
duals = list(g4.blades_lst[-(g4.n + 1):-1])
In [7]:
for dual in duals:
show(dual)
After reverse, the i-th of them is exactly $e_{1} \wedge \ldots \wedge \breve{e}_{i} \wedge \ldots \wedge e_{n}$
In [8]:
duals.reverse()
In [9]:
for dual in duals:
show(dual)
Turn them into base reprsentation:
In [10]:
for dual in duals:
show(g4.blade_to_base_rep(dual))
In [11]:
show(g4.e_obj)
In [12]:
for i,dual in enumerate(duals):
show(S(-1)**(i-1)*g4.base_to_blade_rep(g4.blade_to_base_rep(dual) * g4.e_obj))
In [13]:
for i,dual in enumerate(duals):
show(collect(expand(S(-1)**(i-1)*g4.base_to_blade_rep(g4.blade_to_base_rep(dual) * g4.e_obj)),g4.blades_lst))
So it turns out the formula used is actually
$$e^{i}=(-1)^{i-1}\left(e_{1} \wedge \ldots \wedge \breve{e}_{i} \wedge \ldots \wedge e_{n}\right) E_{n}$$it's different from
$$e^{i}=(-1)^{i-1}\left(e_{1} \wedge \ldots \wedge \breve{e}_{i} \wedge \ldots \wedge e_{n}\right) E_{n}^{-1}$$by a factor of $E_{n}^{-2}$ which is a scalar.
$E_n$ is
In [14]:
g4.e
Out[14]:
$E_n^2$ is
In [15]:
g4.e * g4.e
Out[15]:
In [16]:
g4.e_sq
Out[16]:
$E_n^{-2}$ is
In [17]:
1 / g4.e_sq
Out[17]:
It's clear that $$e_{i} \cdot e^{j}={E_n}^{2}\delta_{i}^{j}$$
In [18]:
for i,base in enumerate(g4.basis):
show(g4.dot(g4.basis[i], g4.r_basis[i]))
And g_inv is implemented like this:
# Calculate inverse of metric tensor, g^{ij}
for i in self.n_range:
rx_i = self.r_symbols[i]
for j in self.n_range:
rx_j = self.r_symbols[j]
if j >= i:
g_inv[i, j] = self.dot(self.r_basis_dict[rx_i], self.r_basis_dict[rx_j])
if not self.is_ortho:
g_inv[i, j] /= self.e_sq
else:
g_inv[i, j] = g_inv[j, i]
which divide the dot product by $E_n^2$ and it's not enough and caused:
$$g^{il} g_{lk} \neq \delta^i_k $$
In [21]:
g4.g*g4.g_inv
Out[21]:
The reason is that g_inv can be seen as:
In [24]:
g_inv = eye(g4.n)
g4.dot_mode = '|'
for i in g4.n_range:
rx_i = g4.r_symbols[i]
for j in g4.n_range:
rx_j = g4.r_symbols[j]
if j >= i:
if g4.is_ortho:
g_inv[i, j] = g4.dot(g4.r_basis_dict[rx_i], g4.r_basis_dict[rx_j])
else:
# NOTE: both reciprocal basis vectors should be devided by E_n^2
g_inv[i, j] = g4.dot(g4.r_basis_dict[rx_i] / g4.e_sq, g4.r_basis_dict[rx_j] / g4.e_sq)
else:
g_inv[i, j] = g_inv[j, i]
Now g_inv is correct and $$g^{il} g_{lk} = \delta^i_k $$
In [28]:
g4.g*g_inv
Out[28]:
In [ ]: