# Geodesics/tangent vector evaluation interpolation error

Hi there. I am trying to evaluate the tangent vector at a specific value of the affine parameter. I have tried everything from https://doc.sagemath.org/html/en/refe... but it raises an error i can't find a solution to.

edit: Just to be precise and avoid taking up anyone's time. I am actually trying to create a poincare section of my orbit so i think i could use the numerical solution by calculating the velocity where the sol meets the criteria. Maybe there's an easier way (?). Couldnt think of anything else that can help me find the poincare section after and not during the integration.

M = Manifold(4, 'M', latex_name=r'\mathcal{M}', structure='Lorentzian')
X.<t,r,th,ph> = M.chart(r't r:(2.0001,+inf) th:(0,pi):\theta ph:(0,2*pi):periodic\phi')
X.coord_range()

var('m,b,a', domain='real')
g = M.metric()
m=1

g[0,0] = -(1-2*m*r/(r^2 + (a*cos(th))^2))
g[0,3] = -2*a*m*r*sin(th)^2/(r^2 + (a*cos(th))^2)
g[1,1] = (r^2 + (a*cos(th))^2)/(r^2 -2*m*r + a^2)
g[2,2] = (r^2 + (a*cos(th))^2)
g[3,3] = (r^2+a^2+2*m*r*(a*sin(th))^2/(r^2 + (a*cos(th))^2))*sin(th)^2
g.display()

def g00(m,r,th,a):return g[0,0](m,r,th,a)
def g03(m,r,th,a):return g[0,3](m,r,th,a)
def g11(m,r,th,a):return g[1,1](m,r,th,a)
def g22(m,r,th,a):return g[2,2](m,r,th,a)
def g33(m,r,th,a):return g[3,3](m,r,th,a)
def D(m,r,th,a): return (g03(m,r,th,a))^2-g00(m,r,th,a)*g33(m,r,th,a)

E.<x,y,z> = EuclideanSpace()
phi = M.diff_map(E, [r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)])
phi.display()

def initial_vector(r0, b,al, ph0=0, E=1, inward=False):
t0,th0=0,pi/2
L = -b*E
vt0=(E*g33(m,r0,th0,a)+L*g03(m,r0,th0,a))/(D(m,r0,th0,a))
vr0=sqrt(((1/D(m,r0,th0,a))*((L^2)*g00(m,r0,th0,a)+2*E*L*g03(m,r0,th0,a)+(E^2)*g33(m,r0,th0,a)))/g11(m,r0,th0,a))
if inward:
vr0 = - vr0
vth0 = al/(100*r0)
vph0 = -(1/D(m,r0,th0,a))*(E*g03(m,r0,th0,a)+L*(g00(m,r0,th0,a)))
p0 = M((t0, r0, th0, ph0), name='p_0')
return M.tangent_space(p0)((vt0, vr0, vth0, vph0), name='v_0')

r0 = 100
th0=pi/2
s = var('s')

v0 = initial_vector(r0, al=1,b=1, ph0=0, inward=True)
geod = M.integrated_geodesic(g, (s, 0, 500), v0, across_charts=True)
sol = geod.solve_across_charts(step=0.01,method='odeint',parameters_values={a:1},
solution_key='sol',verbose=True)

interp = geod.interpolate(solution_key='sol',
interpolation_key='interp 1', verbose=True)


where using:

v = geod.tangent_vector_eval_at(200, verbose=True)


raises :

Evaluating tangent vector components from the interpolation associated with the key 'interp 1' by default...

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-24-6fbce2379e9a> in <module>
----> 1 v = geod.tangent_vector_eval_at(Integer(200), verbose=True)

/opt/sagemath-9.2/local/lib/python3.7/site-packages/sage/manifolds/differentiable/integrated_curve.py in tangent_vector_eval_at(self, t, interpolation_key, verbose)
2275             # partial test, in case future interpolation objects do not
2276             # contain lists of instances of the Spline class
-> 2277             raise TypeError("unexpected type of interpolation object")
2278
2279         interpolated_coordinates=[coordinate_curve_spline(t)

TypeError: unexpected type of interpolation object

edit retag close merge delete

Can you share the whole code that you are using, so that one can reproduce your issue? In particular, what are M and initial_vector?

Of course. I've edited both into the original post. Haven't really used anything that different than the SageManifolds' examples.

I am actually trying to create a poincare section of my orbit so i think i could use the numerical solution by calculating the velocity where the sol meets the criteria. Maybe there's an easier way (?). Couldnt think of anything else that takes advantage of the fact that i already have the geodesic solution after the integration.

The code is still incomplete; m is not defined and the metric is not initialized.

i've included everything on the original post! thanks for your time! really appreciate it.

Sort by » oldest newest most voted

This is a bug, now tracked here. Thank you for reporting it.

It seems there was an oversight when adding the solve_across_chart option, the method tangent_vector_eval_at seems to have been forgotten. (I am the one to blame).

You have 2 possible workarounds here:

1/ Either set the option across_charts to False in the manifold declaration, and use geod.solve instead of geod.solve_across_charts. Doing this you will lose the option to detect the edge of the chart during integration (that you are using to avoid the singularity I assume), and the option to switch charts (that you don't seem to use in this example).

 geod._interpolations['interp 1'] = geod._interpolations['interp 1']


it will only work if the integration is done in single chart (which is the case here), and if the point you are looking for is inside this chart (not the case for your geod.tangent_vector_eval_at(200), but up to geod.tangent_vector_eval_at(97, verbose=True) works fine.

more