# Integrating geodesics backwards in time

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.

edit retag close merge delete

Sort by » oldest newest most voted

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.

more

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