Ask Your Question
2

Geodesics/tangent vector evaluation interpolation error

asked 2021-05-08 12:50:02 +0200

Bubbz gravatar image

updated 2021-05-09 18:53:28 +0200

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 flag offensive close merge delete

Comments

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?

eric_g gravatar imageeric_g ( 2021-05-09 12:16:01 +0200 )edit

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

Bubbz gravatar imageBubbz ( 2021-05-09 12:40:11 +0200 )edit

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.

Bubbz gravatar imageBubbz ( 2021-05-09 13:00:10 +0200 )edit

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

eric_g gravatar imageeric_g ( 2021-05-09 18:07:40 +0200 )edit

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

Bubbz gravatar imageBubbz ( 2021-05-09 18:54:35 +0200 )edit

1 Answer

Sort by ยป oldest newest most voted
2

answered 2021-05-11 16:06:36 +0200

Florentin Jaffredo gravatar image

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).

2/ More convenient in your case: between your last 2 lines, add:

 geod._interpolations['interp 1'] = geod._interpolations['interp 1'][0][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.

edit flag offensive delete link more

Comments

Workaround works just fine ! thank you very much!

Bubbz gravatar imageBubbz ( 2021-05-14 14:30:03 +0200 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2021-05-08 12:50:02 +0200

Seen: 85 times

Last updated: May 11