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

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

• I wish I could thank both(@LutzL and @JJacquelin) of you more for your efforts. Thanks. Kudos! – mm-crj Dec 17 '18 at 19:06

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

• Shouldn't the 2nd boundary condition be $u(1/2)-u(0)=p$? – mm-crj Dec 16 '18 at 19:27
• No, that would leave the integration parameter arbitrary. Fixing $u(0)=0$ so that $u(t)=\int_0^tv(s)^2\,ds$ removes that ambiguity. – Lutz Lehmann Dec 16 '18 at 19:46