I'd like tp quote the linked document "The majority of geodesics on the torus are not æsthetically pleasing. They are aperiodic and cover either the entire surface (if the geodesic is unbounded) or the outer region of the surface bounded by the barrier curves (if bounded). The rare exceptions are the geodesics which return to their starting point after just a few circuits around the z axis". The question arises; whay is found be you?

Also see that thread. There are other references (not only linked by you) there.

Can you kindly elaborate your claim "The crossing at g = 0 is spurious", giving us details?

One more question. D[Derivative[1][u][ v] == ((R + r Cos[u[v]]) Sqrt[-g^2 + (R + r Cos[u[v]])^2])/(g r), v] // FullSimplify produces Derivative[2][u][v] == ((g^2 - 2*R^2 - 2*r*Cos[u[v]]*(2*R + r*Cos[u[v]]))*Sin[u[v]]* Derivative[1][u][v])/(g*Sqrt[-g^2 + (R + r*Cos[u[v]])^2]), not u''[v] == ((2 + Cos[u[v]]) (-32 + g^2 - 32 Cos[u[v]] - 8 Cos[u[v]]^2) Sin[u[v]])/g^2. Are eq2 and eq1 equivalent?

BTW, I don't see your eq2 in the document linked by you.

Executing you code, after f = ParametricNDSolveValue[{eq2, u[0] == 0, eq1 /. v -> 0}, u[6 Pi/5] - 3 Pi/5, {v, 0, 6 Pi/5}, g]; Plot[f[g], {g, -6, 6}, AxesLabel -> {g, "f"}, LabelStyle -> {10, Bold, Black}] I obtain "ParametricNDSolveValue::ndinnt: Initial condition -((0.166673 (r+R) Sqrt[-35.9971+(r+R)^2])/r) is not a number or a rectangular array of numbers." I think you omitted the specification of R and r.

@user64494 Thanks for your comments. I added an explanation of why the g = 0 crossing is spurious and also included {r, R} = {2, 4};, which I had omitted. I am aware of the thread you mentioned in your second comment but only included one of its references to save space. My eq2 is not included in that reference. eq1 is equivalent to eq2 in the sense that it is a first integral of eq2. I do not understand the final sentence in your first comment.

Thank you. 1. Unfortunately, the result of Plot[f[g], {g, -6, 6}, AxesLabel -> {g, "f"}, LabelStyle -> {10, Bold, Black}, PlotRange -> All] looks fantastically (compare with the plots in thethred linked by me. ). 2. I am not a specialist in ODEs. Can you kindly elaborate your "eq1 is equivalent to eq2 in the sense that it is a first integral of eq2", giving us details? 3, h[u, v] is not specified by tou so I see empty plots.

@bbgodfrey thank you very much for this thoughtful analytical answer. It is much appreciated. It appears the geodesic approximation of NikiEster which I cite explicitly is reasonable.

Next, should not be ArcLength[%%, {v, 3*Pi/4, 6 Pi/5}] instead of ArcLength[%%, {v, 0, 6 Pi/5}] (see the question under consideration)? The former produces 4.97607.

@user64494 Multiplying both sides of eq2 by 2 u'[v] and integrating yields eq1. Your last comment is, I believe` incorrect.

I repeat, you don't define h[u, v], Is it h[u_, v_] = {(R + r Cos[u]) Cos[v], (R + r Cos[u]) Sin[v], r Sin[u]}; from the question?

Sorry for the typos. It should be "The question arises; what is found by you?" in my first comment to your answer.

I got the range {0, 6 Pi/5}, but then u[0] == 3*Pi/5 should be. This produces ArcLength equal to 12.8527.

@bbgodfrey I should stated have more correctly state one of the three geodesics was reasonably approximated. Again, thank you for the exposition. Although the question is a duplicate your answer provides the most comprehensive assessment.

@bbgodfrey: Sorry, somewhat more accurate execution gives ArcLength is equal to 9.94234. Details on demand.

@user64494 I can at least answer two parts of your question: 1. As to eq2, notice Integrate cannot handle equation (… == …) at least for now, so you need to write e.g. Integrate[ 2 u'[v] #, v] & /@ eq2. Also, pay attention to the line eq2 = ApplySides[Simplify[D[#^2, v]/(2 u'[v]), Trig -> False] &, eq1], that's where eq2 is deduced. (If you're not familiar with ApplySides, etc., please check the document by pressing F1, they're well-documented functions. ) 2. As to h, I think it's common practice in this site to omit definitions that are already in the question.

@user64494 "should not be ArcLength[%%, {v, 3*Pi/4, 6 Pi/5}] instead of ArcLength[%%, {v, 0, 6 Pi/5}]" Do you mean you think the right code should be ArcLength[%%, {v, 3*Pi/4, 6 Pi/5}]? If so, can you explain a bit about why you think the lower limit is 3 Pi/4? （0 and 6 Pi/5 is obviously from f[0] == h[0, 0], f[n] == h[3π/5, 6π/5] in the question. ）

Notice that, first argument of h is u, and second argument of h is v.

@bbgodfrey Feel free to point out if I've made any mistake in explaining :) .

For clarity: 1. My answer aimed to show a different algorithm to calculate a geodesic that was both relatively efficient and accurate. I posted it as a comment and am grateful to @bbgodfrey for his excellent answer

2. I was aware of the site discussing the complexity of geodesics on a torus and the lack of simple approach. I was intimidated by the DE’s but am very glad to learn from the answer provided and I think combined with my extended comment it is instructive.

3. I voted to close as a “duplicate” not as disrespectful to @bbgodfrey. I was concerned by the profusion of comments obfuscating matters. I, therefore, fully support converting to chat (as is advised and usual practice).

4. Ivoted to close as a “duplicate” not as disrespectful to @bbgodfrey. I was concerned by the profusion of comments obfuscating matters. I, therefore, fully support converting to chat (as is advised and usual practice).

@ubpdqn You gave a very good answer, I believe. Certainly, I took no offense at your decision to close the question.

@xzczd Thank you for clarifying several useful points.

@bbgodfrey thank you for clarifying. It is a very interesting (to me) challenge and i think the whole point of site is to stimulate and iteratively refine answers as well revelation of different approaches.

I too found this question an interesting challenge, particularly because I studied General Relativity in graduate school nearly six decades ago. One of the few things I recall is that geodesics can be described by differential equations and that every symmetry reduced the number of equations by one. Fortunately, I did not need to compute the Christoffel Symbols in this case, because others already had done so. Best wishes for the New Year.

@bbgodfrey if I ever get time I will try to deepen my superficial understanding. Currently, looking after 2 of the 4 grandchildren after coming home from work. Best wishes for 2025 to you too.

@bbgodfrey: Sorry, you still don't answer my comment " got the range {0, 6 Pi/5}, but then u[0] == 3*Pi/5 should be. This produces ArcLength equal to 12.8527".

@user64494 I've explained this above.

Since according to the question, f[0] == h[0, 0], f[n] == h[3π/5, 6π/5], and the argument of h is u and v, we have: at the beginning, u==0 and v==0; at the end, u==3Pi/5 and v==6Pi/5.

Notice that, for eq1 and eq2, the independent variable is v and dependent variable is u.

In other words, the boundary conditions for the boundary value problem (BVP) of eq2 are u[0] == 0, u[6 Pi / 5] == 3 Pi / 5, this has been included in the definition of f, please read it carefully.