Can Lagrangian have a potential term proportional to the quadratic or higher of velocity? In general, classical Lagrangian $L(q,\dot q)=\frac{m}{2} \dot q^2-U(q,\dot q)$ has a $\dot q$-dependence. For example, potential term $U(q,\dot q)$ of the charged particle is given as follows:$U(q,\dot q)=e(\phi-\dot q_iA_i)$.
However, I don’t know Lagrangian which has velocity's quadratic and higher terms like:
$L(q,\dot q)=\frac{m}{2} \dot q^2 -b_{ijk}\dot q_i \dot q_j \dot q_k$.
My question is why such Lagrangians do not exist, or why should Lagrangians not contain terms which are proportional to the third or higher order of velocity?
Is there a good reason for this or is it simply because they don’t exist in nature?
 A: For what it's worth, higher powers of velocities (or momenta) do appear in many Lagrangians, e.g. as a power series in relativistic corrections, say, from square root terms.
A: In fact, there are a lot of systems, whose Lagrange contains higher (even infinite) order of velocity in non-relativistic mechanics, for example, if you have a simple equation of motion
$$
\ddot q+c_1 \dot q=0.
$$
Then the Lagrange can be written as
$$
L= \dot q f(\dot q) \exp(-c_2 q)
$$
where $f(\dot q)$ satisfies the equation
$$
\frac{d f}{d\dot q}=\frac{1}{\dot q^2}\exp\left(-\frac{c_2}{c_1} \dot q \right)
$$
where $c_1$ and $c_2$ are two constants.
A: The EL equations for your proposed potential would be
$$m\ddot q -\frac{d}{dt}\frac{\partial U}{\partial \dot q}= \frac{\partial U}{\partial q}$$
which can be re-expressed as
$$\left(m- \frac{\partial^2 U}{\partial \dot q^2} \right)\ddot q = -\frac{\partial U}{\partial q} +\frac{\partial^2 U}{\partial q \partial \dot q} \dot q$$
Non-linear dependance of $U$ on $\dot q$ would therefore result in a effective dynamical mass term, while cross-dependance on $q$ and $\dot q$ would yield a velocity-dependant force, which could describe some form of fluid friction or the Lorentz force, as you mention in the original question.
It should be noted, however, that the addition of such velocity dependance could make the Hessian matrix $\frac{\partial ^2 L}{\partial \dot q^i \partial \dot q^j}$ singular.  In this case, the transition to the Hamiltonian framework via Legendre transfom must be approached with much more caution - see e.g. this question.

Consider the relativistic Lagrangian
$$L = mc^2 \sqrt{1-\frac{v^2}{c^2}} - U(\mathbf x)$$
Expanding to second order in $v^2/c^2$, this becomes
$$L = \frac{1}{2} mv^2 +\frac{3}{8} m v^4/c^2 - U(\mathbf x)$$
which yields the EL equations
$$-\frac{\partial U}{\partial \mathbf x} = \frac{d}{dt}\left[m\left(1+ \frac{3v^2}{2c^2} \right)\mathbf v\right]$$
This is Newtonian mechanics, except that $m$ has been replaced by $m\left(1+\frac{3v^2}{2c^2}\right)$, giving the lowest-order relativistic correction.
