Let $M$ be a Riemmanian manifold, $p\in M$ and $V\in T\_p(M)$. 


Suppose $f^{-1}:U\_p \mapsto U$ is a diffeomorphism of a neighborhood of p to an open subset of $\mathbb{R}$ and define the sequence:
\begin{equation}
\{p\_i := G\_{p\_{i-1}}^{-1}(V)+ p\_{i-1} \}\_{n \in \mathbb{N}},
\end{equation}
where $p\_0=p$ and $G\_q^{-1}$ is the inverse of the pull-back of the metric tensor at $f(p)$ to $p$ in $U$. 


Then it seems *intuitive* to me that $\cup\_{n \in \mathbb{N}}\{f(p\_i)\}$ converges to the image of geodesic in $M$ with velocity $V$, but how can I prove this?


I was thinking of taking a discrete derivatives but still this is not obvious to me. Thank you. 



Convergence of Discrete Geodesic
I have some doubts about this discretization of the geodesic flow. I believe you should have another component of the map: one that transforms the velocity vector $V$. By the way, technically $V$ should be a momentum covector (belonging to the cotangent bundle). 


If you want a discrete version of the geodesic flow, I think you should start from a Lagrangian point of view. The geodesics on a Riemannian manifold with metric tensor $G(x)$ are the solutions of the Euler-Lagrange equations $\frac{d}{dt}\big(\frac{\partial L}{\partial \dot{x}}(x,\dot{x})\big) = \frac{\partial L}{\partial x}(x,\dot{x})$ with Lagrange function $L(x,\dot{x}) = \frac{1}{2}\dot{x}^T G(x) \dot{x}$. Now, a straight-forward discretization of this method is as follows: replace the derivative $\dot{x}$ by a difference $(x\_1 - x)/\varepsilon$ and consider the discrete Lagrangian $L\_{\varepsilon}(x,x\_1) := L(x, (x\_1-x)/\varepsilon)$. By analogy with the continuous case, where the action is $S[x] = \int\_{t\_0}^{t\_1} L(x,\dot{x})dt$ , in the discrete case the action is $S\_{\varepsilon}[x] = \sum\_{k=0}^{N} L\_{\varepsilon}(x\_k,x\_{k+1})$. Then, the critical discrete geodesic should simply be the solution to the algebraic equations $\nabla S[x] = 0$ which componentwise is basically the discrete Euler-Lagrange equations $\frac{\partial L\_{\varepsilon}}{\partial x\_k}(x\_{k-1}, x\_k) + \frac{\partial L\_{\varepsilon}}{\partial x\_k}(x\_{k}, x\_{k+1}) = 0$ for $k=1,...,N$. Now, we can look locally at these discrete equations and write them down using a simplified subscript notation $\frac{\partial L\_{\varepsilon}}{\partial x}(x\_{-1}, x) + \frac{\partial L\_{\varepsilon}}{\partial x}(x, x\_{1}) = 0$. Introduce the variable $p = \frac{\partial L\_{\varepsilon}}{\partial x}(x, x\_{1})$. Then the discrete Euler-Lagrange equations become $\frac{\partial L\_{\varepsilon}}{\partial x}(x\_{-1}, x) +p= 0$ which shifted by one subscript turn into $\frac{\partial L\_{\varepsilon}}{\partial x\_1}(x, x\_1) + p\_1= 0$. Thus we obtain the equations $$p = \frac{\partial L\_{\varepsilon}}{\partial x}(x, x\_{1})$$
$$p\_1 = -\frac{\partial L\_{\varepsilon}}{\partial x\_1}(x, x\_1).$$ If one can express $x\_1$ as a function of $(x,p)$ from the first equation, then the second equation also gives us $p\_1$ as a function of $(x,p)$. Thus we have obtained a map $\Phi : (x,p) \mapsto (x\_1,p\_1)$. Observe that this is a map $\Phi : T^\*M \to T^\*M$, which turns out to be symplectic, because the Lagrangian $L\_{\varepsilon}$ is in fact a generating function of the symplectomorphism $\Phi$. In the current situation $L\_{\varepsilon}(x,x\_1) = \frac{1}{2 \varepsilon^2} (x\_1 - x)^T G(x) (x\_1 - x)$, so the second equation above yields $p\_1 = -\frac{1}{\varepsilon^2} G(x)(x\_1-x)$, which is equivalent to $x\_1 = x - \varepsilon^2 G(x)^{-1}p\_1$, which resembles your equation. But as you can see, it does not involve a constant tangent vector $V$, but rather a variable cotangent vector $p\_1$, which can be found by taking into account the first equation as well. And in this construction you have convergency of $x\_n : n =0,1,2,...$ to a smooth geodesic. 


