I am trying to find hamiltonian for system described by EOM $$ \ddot{x}(t) + \beta \dot{x}^2(t)\sqrt{\dot{x}^2(t)+\dot{y}^2(t)} = 0, \\ \ddot{y}(t) + \beta \dot{y}^2(t)\sqrt{\dot{x}^2(t)+\dot{y}^2(t)} = 0. $$
I wanted to find Lagrangian first and use Legendre transformation but I have no idea how to even start. My first thought was to solve this analytially. Can somebody give me some hints how to start and how to do this?
The Lagrangian is:
$$\mathcal L=\frac 12 (\dot x^2+\dot y^2)-x\,f_x-y\,f_y$$
$\Rightarrow$
the Hamiltonian
$$H=\frac 12\left(\,{p_{{x}}}^{2}\,+{p_{{y}}}^{2}\right)+x{\it f_x}+y{\it f_y}$$
where $$ f_x= \beta \dot{x}^2(t)\sqrt{\dot{x}^2+\dot{y}^2} , \\ f_y=\beta \dot{y}^2(t)\sqrt{\dot{x}^2+\dot{y}^2} $$
I don't think you got the right equations. In one of your comments above you said that basically you want the two dimensional version of the equation $$\ddot{x} \, + \, \beta\, \dot{x}^2 = 0$$ which is the equation for which the magnitude of the air drag force is proportional to the square magnitude of the velocity Well, the higher dimesional version, in vector form, should be $$\ddot{\vec{r}} \, + \, \beta \, \big|\, \dot{\vec{r}} \,\big| \, \dot{\vec{r}} \, = \, \vec{0} $$ In dimension two, the component-wise equations become \begin{align} &\ddot{x} \, + \, \beta \, \dot{x}\, \sqrt{\dot{x}^2 + \dot{y}^2} \, = \, 0\\ &\ddot{y} \, + \, \beta \, \dot{y}\, \sqrt{\dot{x}^2 + \dot{y}^2} \, = \, 0 \end{align} In vector form, your system is a different system: $$\ddot{\vec{r}} \, + \, \beta \, \big|\, \dot{\vec{r}} \,\big|^2 \, \dot{\vec{r}} \, = \, \vec{0} $$
Either way, you can solve analytically any system of the type $$\ddot{\vec{r}} \, + \, \beta \, \big|\, \dot{\vec{r}} \,\big|^m \, \dot{\vec{r}} \, = \, \vec{0} $$ where $m$ is a positive integer. First, since position is not explicitly present in the system, we can substitute $\vec{v} = \dot{\vec{r}}$. Then the vector equation becomes $$\dot{\vec{v}} \, + \, \beta \, \big|\, \vec{v} \,\big|^m \, \vec{v} \, = \, \vec{0} $$ or in another notation: $$\frac{d\vec{v}}{dt} \, = \, - \, \beta \, \big|\, \vec{v} \,\big|^m \, \vec{v}$$
One way to go on from here is to change the time-parameter $t$. Assume that $t$ itself can be written as a strictly increasing function $t = t(s)$ of a scalar parameter $s$, so that the function $t(s)$ satisfies the differential equation $$\frac{dt}{ds} \, = \, \frac{1}{\big|\vec{v}\big(t(s)\big)\big|^m}$$ where one can impose the initial condition $t(0) = 0$. Then $$\frac{d\vec{v}}{ds} \, = \, \frac{dt}{ds}\,\frac{d\vec{v}}{dt} \, = \frac{dt}{ds} \left(- \, \beta \, \big|\, \vec{v} \,\big|^m \, \vec{v}\right) \, = \, \frac{1}{\,\,\big|\vec{v}\big|^m} \, \left(- \, \beta \, \big|\, \vec{v} \,\big|^m \, \vec{v}\right) \, = \, - \, \beta\, \vec{v} $$ If you put the last two equations together you get the extended system \begin{align} &\frac{d\vec{v}}{ds} \, = \, - \, \beta\, \vec{v} \\ &\frac{dt}{ds} \, = \, \frac{1}{\,\,|\vec{v}|^m} \end{align} The first vector equation can be solved immediately: $$\vec{v} = \vec{v}(s) = e^{-\beta s}\,\vec{v}_0$$ where $\vec{v}_0$ is a constant vector (initial condition). Knowing $\vec{v}(s)$ as a function of $s$ allows us to substitute it in the second, scalar, equation: $$\frac{dt}{ds} = \frac{1}{\,\,\big|\, e^{-\beta s} \,\vec{v}_0 \,\big|^m} = \frac{e^{m \beta s}}{\,\,\big| \,\vec{v}_0 \,\big|^m}$$ The solution to the latter equation, with initial condition $t(0) = 0$, is $$t = t(s) = \frac{e^{m \beta s} - 1}{\,\,m \beta \,\big| \,\vec{v}_0 \,\big|^m}$$ When you solve for $s$, i.e. invert the function, writing $s = s(t)$ you get $$s = \frac{1}{m \beta}\,\ln\big( 1 + m\beta \, \big| \,\vec{v}_0\big|^m t \big)$$ Substitute the latter in $\vec{v} = e^{-\beta s}\,\vec{v}_0$ and you have the solution of the original vector equation $$\vec{v} = \left(\, e^{-\beta \frac{1}{m \beta}\,\ln\big( 1 + m\beta \, | \,\vec{v}_0|^m t \big)}\,\right)\vec{v}_0 = \big(1\, + \, m\beta \, \big| \,\vec{v}_0\big|^m t \big)^{- \, \frac{1}{m}}\, \vec{v}_0$$ Finally, you can rewrite it as $$\vec{v} \, = \, \left( \frac{1}{\sqrt[m]{1 \, + \, m\beta \big| \,\vec{v}_0\big|^m t \,\, }} \right) \vec{v}_0 $$ To find the position as a function of time $t$, integrate with respect to $t$ the last equation:
$$\vec{r} \, = \, \vec{r}_0 \, + \, \left(\int\, \frac{dt}{\sqrt[m]{1 \, + \, m\beta \big| \,\vec{v}_0\big|^m t \,\, }} \right) \vec{v}_0 $$ If $m > 2$, this integral does not look explicitly solvable (but I cannot bet on it). However, when $m=1$ or $m=2$ you get explicit solutions $$m=1: \,\,\,\,\, \vec{r} \, = \, \vec{r}_0 \, + \, \left(\frac{1}{\beta}\,{\ln\big(1 \, + \, \beta \big| \,\vec{v}_0\big|\, t \big)} \right) \frac{\vec{v}_0}{|\vec{v}_0|} $$ $$m=2: \,\,\,\,\, \vec{r} \, = \, \vec{r}_0 \, + \, \left(\frac{1}{2\beta}\,\sqrt{1 \, + \, 2\beta \big| \,\vec{v}_0\big|^2 t} \right) \frac{\vec{v}_0}{|\vec{v}_0|^2} $$
Working in terms of velocity variables: $v_x=\dot{x}, v_y=\dot{y}$. The second suggestion redily reduces the order of the system: $$ \dot{v}_x + \beta v_x^2\sqrt{v_x^2+v_y^2}=0,\\ \dot{v}_y + \beta v_y^2\sqrt{v_x^2+v_y^2}=0 $$ One now has now multiple options available:
