Derivation of Euler-Lagrange equations for Lagrangian with dependence on second order derivatives Suppose we have a Lagrangian that depends on second-order derivatives:
$$L = L(q, \dot{q}, \ddot{q},t).\tag{1}$$
If we're working on the variational problem for this Lagrangian, then I know that we'll wind up with the following Euler-Lagrange equation:
$$\frac{\partial L}{\partial q} - \frac{d}{dt} \frac{\partial L}{\partial \dot{q}} + \frac{d^2}{dt^2} \frac{\partial L}{\partial \ddot{q}} = 0.\tag{2}$$
However, I can't see how to derive this equation. Obviously, the final term is supposed to come from integrating by parts the contribution of the $\ddot{q}$ dependence to the variation in the Lagrangian; doing that yields (writing $T$ for the period of time over which we're extremising the action):
$$\int_T \frac{\partial L}{\partial \ddot{q}} \delta \ddot{q} = \left[ \frac{\partial L}{\partial \ddot{q}} \delta \dot{q} - \frac{\partial L}{\partial \ddot{q}} \delta q \right]_{\partial T} + \int_T \frac{d^2}{dt^2} \frac{\partial L}{\partial \ddot{q}}.\tag{3}$$
Now, the term in square brackets presumably has to vanish. The right-hand term will do so, since $\delta q$ vanishes on the boundary of $T$; but why should the left-hand term vanish? Is it just a condition of dealing with the variations for such a Lagrangian that we consider only variations for which $\delta \dot{q}$ vanishes on the boundary of integration as well? Or is there something I'm missing?
 A: *

*A Lagrangian $L(q,\dot{q},\ddot{q}, \ldots,q^{(N)},t)$ [that depends on up to $N$th-order time derivatives] has an infinitesimal variation of the form
$$ \begin{align}\delta L&~~~=~\sum_{k=0}^N\sum_{j=1}^n
\frac{\partial L}{\partial q^{j(k)}}\delta q^{j(k)}, \qquad q^{j(k)}~:=~\frac{d^kq^j}{dt^k},\cr
&~~~\stackrel{(D)}{=}~
\sum_{k=0}^N\sum_{j=1}^n\left(\frac{d P_{(k)j}}{dt}+P_{(k-1)j}\right)\delta q^{j(k)}\cr 
&\stackrel{\text{Leibniz}}{=}~P_{(-1)j}~ \delta q^j+  \frac{d}{dt}\sum_{k=0}^{N-1}\sum_{j=1}^nP_{(k)j}~\delta q^{j(k)}.\end{align}\tag{A} $$
In eq. (A) we have used that infinitesimal variations $\delta$ and time-derivatives $\frac{d}{dt}$ commute, and we have defined a sequence of quantities
$$ P_{(k)j}~:=~\sum_{m=k+1}^N\left(-\frac{d}{dt} \right)^{m-(k+1)}\frac{\partial L}{\partial q^{j(m)}}, $$
$$\qquad k~\in~\{-1,0,1,2,\ldots\} , \qquad j~\in~\{1,\ldots,n\} .\tag{B}$$
[The tail $P_{(N)j}$, $P_{(N+1)j}$, $\ldots$ of the sequence (B) vanish identically.] In the sequence (B), the first element
$$P_{(-1)j}~\stackrel{(B)}{:=}~\sum_{m=0}^N\left(-\frac{d}{dt} \right)^{m}\frac{\partial L}{\partial q^{j(m)}}, \qquad j~\in~\{1,\ldots,n\}, \tag{C}$$
is the Euler-Lagrange (EL) expression itself; and the subsequent elements $P_{(0)j}$, $P_{(1)j}$, $P_{(2)j}$, $\ldots$ are the higher Ostrogradsky momenta$^1$; which satisfy a recurrence relation
$$ \frac{d P_{(k)j}}{dt}+P_{(k-1)j}~\stackrel{(B)}{=}~\frac{\partial L}{\partial q^{j(k)}}, $$
$$\qquad k~\in~\{0,1,2,\ldots\}, \qquad j~\in~\{1,\ldots,n\} . \tag{D}$$

*The corresponding action $$S[q]~=~\int_{t_i}^{t_f} \! dt ~L(q,\dot{q},\ddot{q}, \ldots,q^{(N)},t)\tag{E}$$ 
has an infinitesimal variation of the form
$$ \delta S~=~\int_{t_i}^{t_f} \! dt~\delta L
~\stackrel{(A)}{=}~\text{bulk-terms}+\text{boundary-terms},\tag{F}$$
where
$$\text{bulk-terms}~=~\int_{t_i}^{t_f} \! dt \sum_{j=1}^n P_{(-1)j}~ \delta q^j,\tag{G}$$
and 
$$\text{boundary-terms}
~=~\sum_{k=0}^{N-1}\sum_{j=1}^n\left[P_{(k)j}~\delta q^{j(k)} 
 \right]_{t=t_i}^{t=t_f}.\tag{H}$$

*Now, in order to deduce the EL equations$^2$
$$ P_{(-1)j}~\stackrel{(G)}{\approx}~0,\qquad j~\in~\{1,\ldots,n\}, \tag{I}$$
which are $n$ ODEs of $2N$th order, we need the boundary-terms (H) to vanish by specifying $2Nn$ boundary conditions (BCs), i.e. $Nn$ initial conditions and $Nn$ final conditions. There are various possibilities:


*

*Essential BCs: $q^{j(k)}$ are fixed at the boundary, where $j\in\{1,\ldots,n\}$ and $k\in\{0,\ldots,N\!-\!1\}$.

*Natural BCs: $P_{(k)j}$ vanish at the boundary, where $j\in\{1,\ldots,n\}$ and $k\in\{0,\ldots,N\!-\!1\}$.

*Some combination of essential and natural BCs.


*So far we have only discussed the higher Lagrangian formulation with an $N$th-order Lagrangian. There is also a higher Hamiltonian formulation with independent phase space variables $q^{j(k)}$ and $P_{(k)j}$, where $j\in\{1,\ldots,n\}$ and $k\in\{0,\ldots,N\!-\!1\}$.
--
$^1$ NB: The labelling of higher Ostrogradski variables is shifted by one as compared to Wikipedia's notation.
$^2$ Here the $\approx$ symbol means equality modulo eqs. of motion.
A: Suppose now the Lagrangian $L$ is a function of $y(x), y'(x)$ and also $y''(x)$, i.e. it contains second derivatives w/r to the parameter $x$.
It is straightforward to adapt the 
usual procedure to this case: write
\begin{align}
Y(x,\epsilon)=y(x)+\epsilon\,\eta(x)
\end{align}
for an otherwise arbitrary function $\eta$.
We then have the parametrized integral
\begin{align}
I(\epsilon)=\displaystyle
\int_a^b\,dx\,L(x,Y(x,\epsilon),Y'(x,\epsilon),
Y''(x,\epsilon))
\end{align}
and we want to find $L$ at $\epsilon=0$
so that 
\begin{align}
0&=\frac{dI}{d\epsilon}\vert_{\epsilon=0}\, ,\\
&=\displaystyle\int_a^b\,dx\,
\left(
\frac{\partial L}{\partial Y}
\frac{\partial Y}{\partial \epsilon}+
\frac{\partial L}{\partial Y'}\frac{\partial Y'}{\partial \epsilon}
+\frac{\partial L}{\partial Y''}\frac{\partial Y''}{\partial \epsilon}
\right)\vert_{\epsilon=0}\, ,\\
&=\displaystyle\int_a^b\,dx\,
\left(
\frac{\partial L}{\partial y}\eta
+\frac{\partial L}{\partial y'}\eta'
+\frac{\partial L}{\partial y''}
\eta'' \right)\, .
\end{align}
We need to turn around the terms in 
$\eta'$ and $\eta''$.  A first integration
by parts will do this:
\begin{align}
\displaystyle\int_a^b\,dx\,
\frac{\partial F}{\partial y'}\frac{d\eta}{dx}
&=\frac{\partial F}{\partial y'}\eta(x)\Bigl\vert_a^b-
\displaystyle\int_a^b\,dx\,\eta
\frac{d}{dx}\left(\frac{\partial F}{\partial y'}\right)\, ,\\
\displaystyle\int_a^b\,dx\,
\frac{\partial F}{\partial y''}
\frac{d\eta'}{dx}&=
\frac{\partial F}{\partial y''}\eta'(x)
\Bigl\vert_a^b-
\displaystyle\int_a^b\,dx\,\eta'
\frac{d}{dx}\left(\frac{\partial F}{\partial y''}\right)\, .
\end{align}
We assume now that the function $\eta$ is
chosen so that $\eta(b)=\eta(a)=0$ as before.
In addition, we must also assume
that $\eta'(b)=\eta'(a)=0$, a new condition.
In this way we have
\begin{align}
\displaystyle\int_a^b\,dx\,
\frac{\partial L}{\partial y'}\eta'&=
-\displaystyle\int_a^b\,dx\,\eta\,
\frac{d}{dx}
\left(\frac{\partial L}{\partial y'}\right)\, ,\\
\displaystyle\int_a^b\,dx\,
\frac{\partial L}{\partial y''}\eta''&=
-\displaystyle\int_a^b\,dx\,\eta'\,
\frac{d}{dx}\left( \frac{\partial L}{\partial y''}\right)\, .
\end{align}
We still need to turn $\eta'$ around one last time. Using integration by parts again:
\begin{align}
-\displaystyle\int_a^b\,dx\,\eta'
\frac{d}{dx}\left(\frac{\partial F}{\partial y''}\right)= + \displaystyle\int_a^b\,dx\,
\eta\,\frac{d^2}{dx^2}\left(\frac{\partial F}{\partial y''}\right)
\end{align}
where the boundary condition $\eta(b)=\eta(a)$ has been used to eliminate the boundary term.
Thus, putting all this together, we get
\begin{align}
0=\frac{dI}{d\epsilon}\Bigl\vert_{\epsilon=0}
=\displaystyle\int_a^b \,dx\,\eta\,
\left(\frac{d^2}{dx^2}\left(\frac{\partial F}{\partial y''}\right)-\frac{d}{dx}
\left(\frac{\partial F}{\partial y'}\right)
+\frac{\partial F}{\partial y}\right)\, .
\end{align}
Since $\eta$ is arbitrary (up to the boundary conditions), we find therefore the function 
$L$ must satisfy the differential equation
\begin{align}
0=\frac{d^2}{dx^2}\left(
\frac{\partial L}{\partial y''}\right)-
\frac{d}{dx}\left(\frac{\partial L}{\partial y'}\right)+\frac{\partial L}{\partial y}\, .
\end{align} 
The generalization to $L$ containing yet
more derivatives is obvious: for
a the derivative of order $k$ we obtain a 
sign $(-1)^k$ as we need $k$ integrations
by parts.  Thus, we obtain the generalized
Euler-Lagrange equation in the form
\begin{align}
0&=\sum_{k}(-1)^k\frac{d^k}{dx^k}
\left(\frac{\partial L}{\partial y^k}\right)
\equiv E(L)\, ,\\
E&=\sum_k (-1)^k \frac{d^k}{dx^k}
\frac{\partial }{\partial y^k}\, .
\end{align}
There's some discussion of this (including how to properly define a conjugate momentum) in the following:

Riahi, F. "On Lagrangians with higher order derivatives." American Journal of Physics 40.3 (1972): 386-390.

and also in 

Borneas, M. "On a generalization of the Lagrange Function." American Journal of Physics 27.4 (1959): 265-267.

A: Since a user has not converted their comment which almost answers the question into an answer, I'm making this community-wiki post in case the comment is deleted:

For a Lagrangian that depends on first-order derivatives, we will find a second-order equation of motion. For such an equation we need two boundary conditions --- for instance, the position of the particle at an initial and final time. This condition 'fixes q at the end points'. For a Lagrangian that depends on second-order derivatives, we will find a fourth-order equation of motion, in general. So we will need four boundary conditions, and fixing the velocity (as well as the position) at the initial and final time will do the trick

