Ask Your Question
4

Integrating geodesics backwards in time

asked 2021-03-30 10:51:29 +0100

Bubbz gravatar image

updated 2021-03-30 11:44:30 +0100

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.

edit retag flag offensive close merge delete

1 Answer

Sort by » oldest newest most voted
4

answered 2021-03-30 16:32:33 +0100

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.

edit flag offensive delete link more

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 ( 2021-03-31 11:55:48 +0100 )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-03-30 10:51:29 +0100

Seen: 342 times

Last updated: Mar 30 '21