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.

| cite | improve this answer | |
  • $\begingroup$ I wish I could thank both(@LutzL and @JJacquelin) of you more for your efforts. Thanks. Kudos! $\endgroup$ – 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.

enter image description here enter image description here

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.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.plot(sol[0], sol[1]); plt.xlabel("p"); plt.ylabel("$\lambda$"); plt.grid(); 
| cite | improve this answer | |
  • $\begingroup$ Shouldn't the 2nd boundary condition be $u(1/2)-u(0)=p$? $\endgroup$ – mm-crj Dec 16 '18 at 19:27
  • 1
    $\begingroup$ 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. $\endgroup$ – Lutz Lehmann Dec 16 '18 at 19:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.