First time here? Check out the FAQ!

Ask Your Question
4

Integrating geodesics backwards in time

asked 3 years ago

Bubbz gravatar image

updated 3 years ago

slelievre gravatar image

Is it possible to integrate a geodesic equation backwards in time in SageMath? I'm using the following code in order to ray trace null geodesics:

v0 = initial_vector(r0, b, al, ph0=0, inward=False)
geod = M.integrated_geodesic(g, (s, 50, 0), v0, across_charts=True)
sol = geod.solve_across_charts(step=0.1, parameters_values={a: 0, q: 0},
                               verbose=True)

where 's' is my affine parameter.

From what i've understood, solve_across_charts does use scipy.odeint which supports backwards integration but when i try to run it i get the following error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-58-13e5246d5cee> in <module>
     14         geod = M.integrated_geodesic(g, (s, Integer(50), Integer(0)), v0, across_charts=True)
     15         sol = geod.solve_across_charts(step=RealNumber('0.1'),parameters_values={a:Integer(0),q:Integer(0)},
---> 16                                        verbose=True) 
     17 
     18 

/opt/sagemath-9.2/local/lib/python3.7/site-packages/sage/manifolds/differentiable/integrated_curve.py in solve_across_charts(self, charts, step, solution_key, parameters_values, verbose, **control_param)
   1719         ics = initial_pt_coords + initial_tgt_vec_comps
   1720         times = np.linspace(t_min, t_max, int((t_max - t_min) / step) + 1,
-> 1721                             endpoint=True)
   1722         nt = len(times)
   1723 

<__array_function__ internals> in linspace(*args, **kwargs)

/opt/sagemath-9.2/local/lib/python3.7/site-packages/numpy/core/function_base.py in linspace(start, stop, num, endpoint, retstep, dtype, axis)
    113     num = operator.index(num)
    114     if num < 0:
--> 115         raise ValueError("Number of samples, %s, must be non-negative." % num)
    116     div = (num - 1) if endpoint else num
    117 

ValueError: Number of samples, -499, must be non-negative.

Which does make sense. Any workaround? I've been also trying to set my affine parameter as a negative but this also raises an error.

Thanks in advance.

Preview: (hide)

1 Answer

Sort by » oldest newest most voted
4

answered 3 years ago

Florentin Jaffredo gravatar image

Hi!

Yes, it should be possible. You can find an example in this notebook.

It's done by defining the initial 4-velocity with a negative time component (cells 26-27 in the example).

You didn't post your full code I can't be completely sure, but my guess is that the integration parameter inside integrated_geodesic is expected to be increasing. ​Could you try to replace the second line in your code with

geod = M.integrated_geodesic(g, (s, 0, 50), -v0, across_charts=True)

It's a bit of a workaround, but it should work. I didn't expect people to call solve_across_charts in this way, but it's in fact logical. I think we should support it (or at least throw an error when declaring the geodesic). It's probably just a matter of adding some absolute values in some places, and maybe fixing the chart limit detection.

Preview: (hide)
link

Comments

It does seem to do the trick for the time being. Simulations still running but so far so good. Thanks a lot!

Bubbz gravatar imageBubbz ( 3 years ago )

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: 3 years ago

Seen: 379 times

Last updated: Mar 30 '21