System of ODEs with integral constrains Can someone point me in a direction to solve this kind of integral constrained system of ODEs. As far as I know, there are no analytic methods that can solve this. So I will resort to numerical methods. But I can't break it down to a system of differential equations with algebraic constraints which can be solved numerically.
\begin{align}
&\int_0^{1/2}\dot{y}^2(t)=p\\
&2\lambda_1\ddot{y}(t)+\pi cos(\pi y(t))=0\\
&y(0)=0,y(1/2)=1/2
\end{align}
I have reduced it to 1st order:
\begin{align}
&\int_0^{1/2}x^2(t)=p\\
&\dot{y}=x \\
&2\lambda_1\dot{x}(t)+\pi cos(\pi y(t))=0\\
&y(0)=0,y(1/2)=1/2
\end{align}
but its still not suitable for a numerical solution. Any help will be appreciated.
Edit: P is a constant known, and $\lambda_1$ is a constant that has to be determined.
 A: $$2\lambda_1y''+\pi\cos(\pi y)=0$$
$$2\lambda_1y''y'+\pi\cos(\pi y)y'=0$$
$$\lambda_1(y')^2+\sin(\pi y)=c_1$$
$$y'=\frac{dy}{dt}=\sqrt{\frac{c_1-\sin(\pi y)}{\lambda}} \tag 1$$
Condition $$p=\int_{t=0}^{t=1/2}\left(\frac{dy}{dt}\right)^2dt=\int_{y(0)}^{y(1/2)}\frac{dy}{dt}dy=\int_0^{1/2}y'dy$$
$$p=\int_0^{1/2}\sqrt{\frac{c_1-\sin(\pi y)}{\lambda}}dy$$
$$p=-\frac{2}{\pi}\sqrt{\frac{c_1-1}{\lambda}}\text{E}\left(\frac{\pi}{4}\:\bigg|\;\frac{-2}{c_1-1}\right)$$
E$(\Phi\:|\:k)$ is the elliptic integral of the second kind with $\Phi=\frac{\pi}{4}$ and $k=\frac{-2}{c_1-1}$
http://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html
Solving $p=-\frac{2}{\pi}\sqrt{\frac{c_1-1}{\lambda}}\text{E}\left(\frac{\pi}{4}\:\bigg|\;\frac{-2}{c_1-1}\right)$ for $c_1$ leads to $c_1=c_1(p)$.
As far as I know, there is no standard closed form for the inverse function of 
$f\left(x\text{E}\left(\frac{\pi}{4}\:\big|\:\frac{1}{x}\right)\right)$. So, we cannot express $c_1$ as a function of $p$ on closed form. Numerical calculus is required. At this stage of the calculus we can consider that $c_1$ is know (as far as $p$ is a given value).
$$t=\pm\int \sqrt{\frac{\lambda}{c_1-\sin(\pi y) }}\:dy+\text{constant}$$
For $t\geq 0$ and $y\geq 0$  the condition $y(0)=0$ implies :
$$t=\int_0^y \sqrt{\frac{\lambda}{c_1-\sin(\pi \xi) }}\:d\xi$$
With $y(1/2)=1/2$ :
$$t=\frac12+2\sqrt{\frac{\lambda}{c_1-1}}\text{F}\left(\frac{\pi}{4}(1-2y)\:\bigg|\:\frac{-2}{c_1-1} \right) \tag 2$$
$\text{F}(\phi\:|\:k)$ is elliptic integral of the first kind with $\phi=\frac{\pi}{4}(1-2y)$ and $k=\frac{-2}{c_1-1}$ http://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html
The inverse function $y(t)$ involves the Amplitude Jacobi elliptic function. http://mathworld.wolfram.com/JacobiAmplitude.html
This is an arduous calculus. With the help of WolframAlpha :
$$y(t)=\frac12-\frac{2}{\pi}\text{am}\left(\frac{\pi}{2}\sqrt{\frac{c_1-1}{\lambda}}(t+c_2)\:\bigg|\:\frac{-2}{c_1-1}\right) \tag 3$$
The condition $y(1/2)=1/2$ implies $c_2=-\frac12$. The result is :
$$y(t)=\frac12-\frac{2}{\pi}\text{am}\left(\frac{\pi}{2}\sqrt{\frac{c_1-1}{\lambda}}(t-\frac12)\:\bigg|\:\frac{-2}{c_1-1}\right) \tag 4$$
NOTE : Eq.$(2)$ seems correct after checking. The Eqs.$(3-4)$ might be not correct. The analytical method is too ugly. The numerical method (such as LutzL did) definitively appears better in practice. 
A: You can reduce the problem to a boundary value problem for which there usually are solvers in a numerics library.
The first order system would be
\begin{align}
\dot y&=v\\
\dot v&=-πcos(πy)/(2λ_1)\\
\dot u &= v^2
\end{align}
with the boundary conditions
\begin{align}
y(0)&=0,& y(1/2)&=1/2\\
u(0)&=0,& u(1/2)&=p.
\end{align}
Tentatively, this can be implemented in python using scipy as
def ev_ode(t,w,param):
    y,v,u = w
    lam = param[0]
    return [ v, -pi*cos(pi*y)/(2*lam), v**2 ]

def ev_bc(w0, wh, param): return [w0[0], wh[0]-0.5, w0[2], wh[2]-p]

t_init = [0, 0.5]
w_init = [ [0,0.5], [1, 1], [0, 0.5] ]
lam_init = [0.3]
res = solve_bvp(ev_ode, ev_bc, t_init, w_init, p=lam_init)
print res.message
print "p =",p,", lambda =", res.p[0]

This problems seems to be very sensitive to initial data. Using $p=0.6$ ($p\ge0.5$ by Cauchy-Schwarz) gave once the successful result
The algorithm converged to the desired accuracy.
p = 0.6 , lambda = 0.26105387754

With this configuration also successful were
The algorithm converged to the desired accuracy.
p = 0.5 , lambda = 5135.44389598
p = 0.7 , lambda = 0.159001268888
p = 0.8 , lambda = 0.114078982598
p = 0.85 , lambda = 0.0994306061876

which seems also to cover the range of admissible parameters, or at least the local interval, as $p=0.9$ did not converge.

Per the computations of JJaqueline, a direct path to a solution is to chose a $c\ge 1$, compute
$$
\frac1{\sqrt{λ}}=\int_0^{1/2}\frac{2\,d\xi}{\sqrt{c-\sin(\pi\xi)}}
$$
and then use the solution of the BVP with these parameters or just $y'=\sqrt{(c-\sin(\pi y))/λ}$ to find the solution $y$ and the integral value.


def p_fun(c):
    res,err = quad(lambda x: 2*(c-sin(pi*x))**-0.5, 0, 0.5);
    lam = res**-2
    def p_ode(w,t): y,u=w; dudt = (c-sin(pi*y))/lam; return [dudt**0.5, dudt ]
    p = odeint(p_ode, [0,0], [0,0.5])[-1,1]
    return p, lam

arr_c = np.linspace(1.001,10,1000)
sol = np.array([ p_fun(c) for c in arr_c]).T

plt.figure(1)
plt.subplot(1,2,1); plt.plot(arr_c,sol[0]); plt.xlabel("c"); plt.ylabel("p"); plt.grid(); 
plt.subplot(1,2,2); plt.plot(arr_c,sol[1]); plt.xlabel("c"); plt.ylabel("$\lambda$"); plt.grid();
plt.figure(2)
plt.plot(sol[0], sol[1]); plt.xlabel("p"); plt.ylabel("$\lambda$"); plt.grid(); 
plt.show()
 .

