System of ODEs with integral constrains

Question

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.

Answer 1

$$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.

Answer 2

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()
 .