Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I want to solve (numerically) the geodesics of the upper half plane and plot. The results are quite known (i) straight lines parallel to $y$-axis and (ii) semicircles centered on the $x$-axis. Now the problem is how to solve the coupled ODEs. I tried with NDSolve but the results are not coming. Can you please write a code and help me.

Here are the coupled ODEs:

$$\begin{align*} \ddot{x} &= \frac{2\dot{x}\dot{y}}{y}\\ \ddot{y} &= \frac{-\dot{x}^2+\dot{y}^2}{y} \end{align*}$$

Thanking you in advance.

share|improve this question
dear chris can you please explain why the t in the NDSolve runs for 0 to 1 and then in parametric plot from 1 to 10. Also If I can put range of t 1 to 10 in NDSolve then two curves are coming ??? can you please explain – ricci1729 May 25 '13 at 10:01
thank you very much :) – ricci1729 May 25 '13 at 10:05

If you use NDSolve as follows

  NDSolve[{ x''[t] == 2 x'[t] y'[t]/y[t], 
          y''[t] == (y'[t]^2 - x'[t]^2)/y[t], 
          x[0] == 1, y[0] == 1, 
          x'[0] == 1, y'[0] == 1}, {x[t], y[t]}, 
           {t, 0, 8}]; 


 ParametricPlot[{x[t], y[t]} /. sol, {t, 0, 8}]


Note that if/when the solution is driven towards 'y=0' you will run into problems.

You might also want to have a look at EquationTrekker

<< EquationTrekker`
share|improve this answer
dear chris; that's very true . y=0 is the singularity that one can immediately see from the ODE's. But my question is, if you change the limit of t from 0 to 9 you see two curves. but why is that is not clear to me??? – ricci1729 May 25 '13 at 10:18
its a bit strange... – chris May 25 '13 at 10:28

I think you can solve the equations analytically. The change of variables y[t] -> Exp[z[t]] changes your system into

x''[t] == 2 x'[t] z'[t]
Exp[2 z[t]] z''[t] == - x'[t]^2

The first equation is easy and you get

z[t] == 1/2 Log[x'[t]] + a

with a a constant; the second equation, assuming x'[t]!=0 gives

1/2 Exp[2 a] x''[t]/x'[t] == -x[t] + b

with b a constant. This can be handled by Mathematica :

solx[t_] = x[t] /. First@
  DSolve[1/2 Exp[2 a] x''[t]/x'[t] == -x[t] + b, x[t], t] /. {C[1] -> c, C[2] -> d}

from which we can retrieve the solution for z[t] :

solz[t_] = 1/2 Log[D[solx[t], t]] + a

In order to check against @chris let's solve for the constants {a,b,c,d} by using the initial conditions :

const = FindRoot[
          {solx[0] == 1, Derivative[1][solx][0] == 1, 
           solz[0] == 0, Derivative[1][solz][0] == 1}, 
           {{a, 1}, {b, 2}, {c, 1}, {d, 1}}] //Chop ;

ParametricPlot[{solx[t] /. const, Exp[solz[t] /. const]}, {t, 0, 5}]


share|improve this answer
thank you very much but i know it is exactly solvable model. the question is how to generate this graph itself from the geodesic equation by NDSolve. – ricci1729 May 28 '13 at 8:00

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.