ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Wed, 31 Mar 2021 11:55:48 +0200Integrating geodesics backwards in timehttps://ask.sagemath.org/question/56427/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.
Thanks in advance.Tue, 30 Mar 2021 10:51:29 +0200https://ask.sagemath.org/question/56427/integrating-geodesics-backwards-in-time/Answer by Florentin Jaffredo for <p>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:</p>
<pre><code>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)
</code></pre>
<p>where 's' is my affine parameter.</p>
<p>From what i've understood, <code>solve_across_charts</code> does use <code>scipy.odeint</code>
which supports backwards integration but when i try to run it
i get the following error:</p>
<pre><code>---------------------------------------------------------------------------
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.
</code></pre>
<p>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.</p>
<p>Thanks in advance.</p>
https://ask.sagemath.org/question/56427/integrating-geodesics-backwards-in-time/?answer=56431#post-id-56431Hi!
Yes, it should be possible. You can find an example in [this notebook](https://nbviewer.jupyter.org/github/sagemanifolds/SageManifolds/blob/master/Notebooks/SM_black_hole_rendering.ipynb).
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.Tue, 30 Mar 2021 16:32:33 +0200https://ask.sagemath.org/question/56427/integrating-geodesics-backwards-in-time/?answer=56431#post-id-56431Comment by Bubbz for <p>Hi!</p>
<p>Yes, it should be possible. You can find an example in <a href="https://nbviewer.jupyter.org/github/sagemanifolds/SageManifolds/blob/master/Notebooks/SM_black_hole_rendering.ipynb">this notebook</a>.</p>
<p>It's done by defining the initial 4-velocity with a negative time component (cells 26-27 in the example).</p>
<p>You didn't post your full code I can't be completely sure, but my guess is that the integration parameter inside <code>integrated_geodesic</code> is expected to be increasing. Could you try to replace the second line in your code with</p>
<pre><code>geod = M.integrated_geodesic(g, (s, 0, 50), -v0, across_charts=True)
</code></pre>
<p>It's a bit of a workaround, but it should work. I didn't expect people to call <code>solve_across_charts</code> 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.</p>
https://ask.sagemath.org/question/56427/integrating-geodesics-backwards-in-time/?comment=56445#post-id-56445It does seem to do the trick for the time being. Simulations still running but so far so good. Thanks a lot!Wed, 31 Mar 2021 11:55:48 +0200https://ask.sagemath.org/question/56427/integrating-geodesics-backwards-in-time/?comment=56445#post-id-56445