Lagrange equation of motions for a particle moving in a surface in the presence of gravity I have to model the dynamic behaviour of a particle solid in a gravitational field.It is for a control theory course. And my background in dynamics is not the greatest. The particle can move left and right (so there are 2 generalized coordinates x,y or $x_1,x_2$ and i have the height in each point in a function $h(x,y) = h(x_1,x_2)$. To help you understand, i have a plot of what i want to achieve at the end (The result is simulated with the final equation of motion).

I have done the analysis in one dimension, where it can move along the x-axis and i have a h(x) function for the height, and the result seemed orderly. There, one can see the inertial, gravitational and centrifugal force. I write the final equation of motion for the 1d case:
$$ m\ddot{x_1} +m\frac{ h'(x_1)h''(x_1) } {1+h'(x_1)^2 } \dot{x}_1^2 + mg  \frac{h'(x_1)}{1+h'(x_1)^2} = 0 $$
The first term is the inertial force, the second the centrifugal force and the third is the gravitational force.
The problem i have is that i cannot make sense of the final equation of motion. I expected more simplified (one lest term).
So i am writing the equation of motions using lagrange formulation:
$$ T = \frac{1}{2} m ( \dot{ x}_{r1}^2 +\dot{ x}_{r2}^2 ) = \frac{1}{2} \dot{x}_r^T [m] \dot{x}_r $$
Here, i want to change variables.$x_r$ is the velocity of the surface and it is not the same as the projected velocity in the x,y plane.
It can be proven that:
$$ ||x|| = cos(\theta ) ||x_r|| = \frac{1}{\sqrt{ 1+||\nabla h||^2} } ||x_r|| \rightarrow ||x_r|| = \sqrt{ 1 + ||\nabla h||^2} ||x|| $$
I leave an image so you can see that (DS = $\nabla h$):
$$ cos( \theta )  = \frac{1}{ \sqrt{ 1+||\nabla h||^2}} $$

So the kinetic energy becomes (From here on, x is a vector $x = \{x_1,x_2\}^T$):
$$ T = \frac{1}{2}(1+ ||\nabla h||^2)m ( \dot{ x_1}^2 +\dot{ x_2}^2 ) = \frac{1+||\nabla h||^2}{2} \dot{x}^T [m] \dot{x}  = T(x(t),\dot{x}(t)) $$
The potential energy of the system is:
$$ V = \int _{x_o} ^{x} mg \nabla h (s) \cdot ds = mg( h(x) - h(x_0) ) = V(x)$$
Let's start differantiating with respect to $x_1$:
$$\frac{ \partial T}{\partial \dot{x_1}} =  (1+ ||\nabla h||^2)m  \dot{ x}_1$$
$$ \frac{ \partial} {\partial t} (\frac{ \partial T}{\partial \dot{x_1}}) =  (1+ ||\nabla h||^2)m  \ddot{ x_1} + m\dot{x_1} \frac{ \partial  ( ||\nabla h||^2)}{\partial t} $$
Here, we should observe that the following is true before we continue:
$$  \frac{ \partial  ( ||\nabla h||^2)}{\partial t} = \frac{ \partial}{\partial t}   \left( \frac{\partial h}{\partial x_1} ^2(x_1(t),x_2(t))+\frac{\partial h}{\partial x_2} ^2(x_1(t),x_2(t)) \right)$$
So:
$$  \frac{ \partial  ( ||\nabla h||^2)}{\partial t} = 2\frac{\partial h}{\partial x_1} \left (\frac{\partial^2 h}{\partial x_1^2} \dot{x_1}+\frac{\partial^2 h}{\partial x_2x_1}\dot{x_2} \right) + 2\frac{\partial h}{\partial x_2} \left (\frac{\partial^2 h}{\partial x_1x_2} \dot{x_1}+\frac{\partial^2 h}{\partial x_2^2}\dot{x_2} \right) = 2 \nabla h ^T H_h \dot{x} $$
So in total the first term of the lagrange equation of motions is:
$$\frac{ \partial} {\partial t} (\frac{ \partial T}{\partial \dot{x_1}}) =  (1+ ||\nabla h||^2)m  \ddot{ x_1} + 2m \left( \nabla h ^T H_h \dot{x} \right) \dot{x_1}$$
The second term :
$$ \frac { \partial (T-V) } {\partial x_1} = \frac{1}{2}m(\dot{x}_1^2+\dot{x}_2^2)  \frac{ \partial }{\partial x_1} \left(1+ ||\nabla h|| ^2  \right)  - mg \frac{\partial h}{ \partial x_1}$$
$$  \frac{ \partial  ( ||\nabla h||^2)}{\partial x_1} = 2\frac{\partial h}{\partial x_1} \frac{\partial^2 h}{\partial x_1^2} +  2\frac{\partial h}{\partial x_2} \frac{\partial^2 h}{\partial x_1x_2} $$
So in total, in one dimension:
$$ (1+ ||\nabla h||^2)m  \ddot{ x_1} + 2m \left( \nabla h ^T H_h \dot{x} \right) \dot{x_1} + mg \frac{ \partial h}{\partial x_1} - m(\dot{x_1}^2+x_2^2)\left(\frac{\partial h}{\partial x_1} \frac{\partial^2 h}{\partial x_1^2} +  \frac{\partial h}{\partial x_2} \frac{\partial^2 h}{\partial x_1x_2}  \right )  $$
In matrix form that becomes :
$$ \ddot{ x} +  2\frac{ \left( \nabla h ^T H_h \dot{x} \right) }{1+ ||\nabla h||^2} \dot{x} + g\frac{1}{1+ ||\nabla h||^2 }\nabla h - ||\dot{x}||^2\frac{1}{1+ ||\nabla h||^2 }   \nabla h ^T H_h = 0 $$
The first term is the inertial force, the second is the centrifugal, however it has a factor of 2, the third is the gravitational force and the fourth, well i do not know. It seems like the centrifugal force too.
Edit Generally, i would like to tread the problem using 2 variables in the whole process. Like, when you study the movement of an object in an inclined plane:

Here, i write the equation of motion in the relative frame (which is constant) and then i change the variables. I could do the same with the kinetic energy.
If the inclined plane has variable angle, so it is like having an $h(x_1)$ i could do the following:
$$ dx_r = \sqrt{ 1 + h'(x_1)^2 } dx_1$$
So the kinetic energy becomes:
$$ T = \frac{1}{2} m \dot{x}_r^2 = \frac{1}{2} m (1 + h'(x_1)^2) \dot{x}_1^2$$

 A: Start by considering the position vector $~\vec R~$ to the mass
$$\vec R=\left[ \begin {array}{c} x_{{1}}\\ x_{{2}}
\\ h \left( x_{{1}},x_{{2}} \right) \end {array}
 \right] \;.
$$
From here, the kinetic energy is:
$$T=\frac{m}{2}\vec v\cdot\vec v\;,
$$
where
$$
\vec v=\vec{\dot{R}}=\frac{\partial \vec R}{\partial \vec q}\,\vec{\dot{q}}\quad,\quad
\vec q= \left[ \begin {array}{c} x_{{1}}\\ x_{{2}}
\end {array} \right]\;,
$$
The kinetic energy can be calculated explicitly in terms of $x_1$ and $x_2$. For example, explicitly, the velocity is:
$$\vec v=\left[ \begin {array}{c} {\dot{x}}_{{1}}\\ {\dot{x}}_
{{2}}\\   \left( {\frac {\partial }{\partial x_{{1}}}
}h \left( x_{{1}},x_{{2}} \right)  \right) {\dot{x}}_{{1}}+ \left( {
\frac {\partial }{\partial x_{{2}}}}h \left( x_{{1}},x_{{2}} \right)
 \right) {\dot{x}}_{{2}}\end {array} \right]$$
The potential energy is:
$$U=m\,g\,h(x_1~,x_2)$$
Then, using $L=T-U$, you obtain the equations of motion in the usual way from the Euler-Lagrange procedure as:
$$\left[ \begin {array}{c} {\ddot x}_{{1}}\\ {\ddot x
}_{{2}}\end {array} \right]= -
 \left[ \begin {array}{c} {\frac { \left( {\frac {\partial }{\partial
x_{{1}}}}h \left( x_{{1}},x_{{2}} \right)  \right)  \left( {{\dot{x}}_{
{1}}}^{2}{\frac {\partial ^{2}}{\partial {x_{{1}}}^{2}}}h \left( x_{{1
}},x_{{2}} \right) +2\,{\dot{x}}_{{1}} \left( {\frac {\partial ^{2}}{
\partial x_{{2}}\partial x_{{1}}}}h \left( x_{{1}},x_{{2}} \right)
 \right) {\dot{x}}_{{2}}+ \left( {\frac {\partial ^{2}}{\partial {x_{{2
}}}^{2}}}h \left( x_{{1}},x_{{2}} \right)  \right) {{\dot{x}}_{{2}}}^{2
}+g \right) }{1+ \left( {\frac {\partial }{\partial x_{{2}}}}h \left(
x_{{1}},x_{{2}} \right)  \right) ^{2}+ \left( {\frac {\partial }{
\partial x_{{1}}}}h \left( x_{{1}},x_{{2}} \right)  \right) ^{2}}}
\\ {\frac { \left( {\frac {\partial }{\partial x_{{2
}}}}h \left( x_{{1}},x_{{2}} \right)  \right)  \left( {{\dot{x}}_{{1}}}
^{2}{\frac {\partial ^{2}}{\partial {x_{{1}}}^{2}}}h \left( x_{{1}},x_
{{2}} \right) +2\,{\dot{x}}_{{1}} \left( {\frac {\partial ^{2}}{
\partial x_{{2}}\partial x_{{1}}}}h \left( x_{{1}},x_{{2}} \right)
 \right) {\dot{x}}_{{2}}+ \left( {\frac {\partial ^{2}}{\partial {x_{{2
}}}^{2}}}h \left( x_{{1}},x_{{2}} \right)  \right) {{\dot{x}}_{{2}}}^{2
}+g \right) }{1+ \left( {\frac {\partial }{\partial x_{{2}}}}h \left(
x_{{1}},x_{{2}} \right)  \right) ^{2}+ \left( {\frac {\partial }{
\partial x_{{1}}}}h \left( x_{{1}},x_{{2}} \right)  \right) ^{2}}}
\end {array} \right]
$$

the EOM's with Newton approach
$$m\,\mathbf J^T\,\mathbf J\,\mathbf{\ddot{q}}=\mathbf J^T\left(
\begin{bmatrix}
   0 \\
   0 \\
   -m\,g \\
 \end{bmatrix}-m\,\mathbf F_z\right)$$
where
$$\mathbf J=\frac{\partial \mathbf R}{\partial \mathbf q}\\
\mathbf F_z=\frac{\partial \mathbf v}{\partial \mathbf q}\,\mathbf{\dot{q}}\\
\mathbf v=\mathbf J\,\mathbf{\dot{q}}$$

\begin{align*}
  &\text{Newton Equation}\\
  &m\,\mathbf{\ddot{R}}=\mathbf{F}+\mathbf{F}_c\tag 1\\
\end{align*}
with
\begin{align*}
  &\mathbf{R}=\mathbf{R}(~\mathbf{q})\quad\Rightarrow\\
  &\mathbf{v}=\mathbf{\dot{R}}=\underbrace{\frac{\partial\mathbf{R}}{\partial\mathbf{q}}}_{\mathbf{J}}\,\mathbf{\dot{q}}
  \quad,
  \mathbf{\ddot{R}}=\mathbf{J}\,\mathbf{\ddot{q}}+
  \underbrace{\frac{\partial\mathbf{v}}{\partial\mathbf{q}}\,\mathbf{\dot{q}}}
  _{\mathbf{a}_q}
\end{align*}
thus equation (1)
\begin{align*}
  & m\,\mathbf{\ddot{R}}=\mathbf{F}+\mathbf{F}_c=
  m\left[~ \mathbf{\mathbf{J}\,\mathbf{\ddot{q}}+\mathbf{a}_q}~\right]=
  \mathbf{F}+\mathbf{F}_c
\end{align*}
rearrange
\begin{align*}
  &m\,\mathbf{\mathbf{J}}\,\mathbf{\ddot{q}}=-m\,\mathbf{a}_q+\mathbf{F}+\mathbf{F}_c
\end{align*}
to eliminate the constraint forces $~\mathbf{F}_c~$ from the above equation, you multiply the  equation from the left with $~\mathbf{J}^T~$
\begin{align*}
  &\boxed{~m\,\mathbf{J}^T\,\mathbf{\mathbf{J}}\,\mathbf{\ddot{q}}=
  -m\,\mathbf{J}^T\mathbf{a}_q+\mathbf{J}^T\,\mathbf{F} ~}
\end{align*}
notice
\begin{align*}
  &\frac{\partial\mathbf{a}}{\partial\mathbf{b}}=\underbrace{\frac{\partial a_i}{\partial b_j}}
  _ {\mathbf A_{~ij}}\quad,\text{Matrix}~n\times m
\end{align*}

*

*$~\mathbf R~$ Position vector

*$~\mathbf q~$ generalized coordinates vector

*$~\mathbf F~$ External forces

*$~\mathbf F~$ Constraint  forces

A: First, I agree that starting with the Lagrangian :
$$L =  \frac 12 (1+\|\nabla h\|^2)\|\dot x\|^2 - gh \tag 1$$
one ends up with your equation of motion (one detail :  the last term should probably be written $H_h \nabla h$ to be a column vector). We can rearrange the different terms to write (setting $m = 1$ for simplicity) :
$$(1+\|\nabla h\|^2)\ddot x + g\nabla h + 2 \dot x (\dot x \cdot H_h \cdot \nabla h) - \|\dot x\|^2 H_h\nabla h=0 \tag{2}$$
The two last terms cannot be interpreted as centrifugal forces, as those appear in rotating frames. Here, those forces appear because the kinetic term in the Lagrangian has this $1+\|\nabla h\|^2$ factor. In the language of differential geometry (or general relativity), this is a metric different from the Euclidean one, which can induce curvature. If we set $g=0$, the equation of motion is the geodesic equation (with the term due to curvature bilinear in $\dot x$).
However, this is not the correct Lagrangian to describe a particle constrained on $\Sigma$. More precisely, the kinetic term is wrong.
Notations : Let us write $2$d vectors like $x = (x_1,x_2)^t$ and $3$d vectors like $\vec x = (x_1,x_2,x_3)^t = (x,x_3)^t$.
In $3$d space, the velocity is :
$$\vec v = \begin{pmatrix} \dot x  \\ \nabla h\cdot \dot x\end{pmatrix}$$
Therefor the kinetic energy is :
$$T = \frac12 \|\vec v\|^2 = \frac 12 (\|\dot x\|^2 + (\nabla h \cdot \dot x)^2)$$
and the full Lagrangian is :
$$L =\frac 12 (\|\dot x\|^2 + (\nabla h \cdot \dot x)^2) - mh \tag 3$$
Then, we can compute :
\begin{align}
\frac{\partial L}{\partial x} &= (\dot x\cdot \nabla h) H\dot x - g\nabla h \\
\frac{\partial L}{\partial \dot x} &= \dot x + \nabla h(\dot x \cdot \nabla h)
\end{align}
so the Euler-Lagrange equation is :
$$\ddot x  +(\ddot x \cdot \nabla h + \dot x \cdot H \cdot \dot x + g)\nabla h = 0 \tag 4$$
To double check that $(4)$ is the correct equation of motion and $(3)$ the correct Lagrangian, we can derive the equation of motion of the particle constrained on $\Sigma$ from the fundamental principle of dynamics.
Fundamental principle of dynamics and normal forces
Since the constraint is not exerting work on the particle, it is enforced by a force normal to the surface $\Sigma$. Some quick calculations tells us that the unit normal vector to $\Sigma$ at a point $\vec x = (x,h(x))$ is :
$$\vec n = \frac{1}{\sqrt{1+\|\nabla h\|^2}} \begin{pmatrix} -\nabla h\\1\end{pmatrix}$$
The acceleration is :
$$\vec a = \begin{pmatrix} \ddot x \\ \nabla h \cdot \ddot x  + \dot x \cdot H_h \cdot \dot x \end{pmatrix}$$
The force of gravity is :
$$\vec g  =\begin{pmatrix} 0 \\ -g\end{pmatrix}$$
Therefore, the fundamental principle of dynamics is $\vec a  = \vec g  +R \vec n$ (where $R$ is the intensity of the normal force), or in other words :
$$\ddot x = -\frac{R}{\sqrt{1+\|\nabla h\|^2}}\nabla h$$
and :
$$\nabla h \cdot \ddot x  + \dot x \cdot H_h \cdot \dot x = -g + \frac{R}{\sqrt{1+\|\nabla h\|^2}}$$
This last equations set $R$ equal to the value which will make the particle stay on $\Sigma$ (in an engineering setting, this would allow you to see at what point the required force exceeds what your materials can take, or maybe you particle is just sitting on the surface so you can only have $R\geq 0$ and when it vanishes the particle jumps off the surface). Its value is :
$$R = \sqrt{1+\|\nabla h\|^2}\left( \nabla h \cdot \ddot x  + \dot x \cdot H_h \cdot \dot x + g \right)$$
Plugging this in the other equation of motion, we get :
$$\ddot x+ \left( \nabla h \cdot \ddot x  + \dot x \cdot H_h \cdot \dot x + g \right)\nabla h = 0$$
which is the same as $(4)$.
The second and third terms in this last equation are due to the fact that the particle has some kinetic energy due to its vertical motion.
Edit : From the discussion with @Eli in the comments below his answer, it appears that we can rewrite this equation in a slightly different way. Equation $(4)$ is rewritten :
$$(I_2 + \nabla h \otimes \text dh) \ddot x + (\dot x H\dot x + g)\nabla h = 0 \tag 5$$
where $\nabla h\otimes \text dh$ is the $2\times 2$ matrix such that  $(\nabla h\otimes \text dh)v = \nabla h (\nabla h \cdot v)$.
Then, we can compute the inverse of the matrix $I_2 + \nabla h \otimes \text dh$, we comes out to :
$$(I_2 + \nabla h \otimes \text dh)^{-1} = I_2 - \frac{\nabla h\otimes \text dh}{1+\|\nabla h\|^2}$$
Then, multiplying $(5)$ by this matrix and using :
$$(I_2 + \nabla h \otimes \text dh)^{-1} \nabla h = \frac{1}{1+\|\nabla h\|^2} \nabla h \tag{6}$$
we get (and this equations is equivalent to $(4)$) :
$$\ddot x +\frac{1}{1+\|\nabla h\|^2}\left( \dot x \cdot H_h \cdot \dot x + g \right)\nabla h = 0$$
This results in a nice equation because of $(6)$ which tells us that $\nabla h$ is a eigenvector of $(I_2 + \nabla h \otimes \text dh)^{-1}$ and the fact that equations $(4)$ is of the form $ (I_2 + \nabla h \otimes \text dh)\ddot x\propto\nabla h$. This would be true if we added an arbitrary potential $V(x)$ to the Lagrangian.
