I am using BVP4C to solve a system of ODEs which is as follows. \begin{equation} \left\{ \begin{aligned} \frac{\partial f(x,y)}{\partial x} &- \frac{d}{ds}\big(\dot{x} f(x,y)\big) = \lambda \ddot{x}(s)\\ \frac{\partial f(x,y)}{\partial y} &- \frac{d}{ds}\big(\dot{y} f(x,y)\big) = \lambda \ddot{y}(s)\\ \end{aligned} \right. \end{equation}

There is a constraint which is of the form $$\dot{x}^2 + \dot{y}^2 = 1$$.

The boundary conditions are $$x(0) = x_A, \, y(0) = y_A, \, x(l) = x_B, y(l) = y_B$$.

What should I do to deal with above constraint?

  • 2
    $\begingroup$ Welcome to SciComp Exchange. A better description of your problem is needed to provide you suggestions that are related to what you want. What is your system of equations? $\endgroup$
    – nicoguaro
    Mar 23, 2015 at 17:18
  • $\begingroup$ @yagoo instead of commenting, please edit your question. Also, please indicate what the boundary conditions are. $\endgroup$ Mar 24, 2015 at 5:28
  • $\begingroup$ @yagoo: This system of equations doesn't look like a system of ODEs; rather, it looks like a system of PDEs. I suppose $s$ is supposed to be the "time" variable? Are you discretizing the PDEs in terms of $x$ and $y$? $\endgroup$ Mar 24, 2015 at 7:56
  • $\begingroup$ @Geoff Oxberry: It looks like a system of PDEs, however, it is a system of ODEs. As you said, $s$ is supposed to be the "time" variable. $\endgroup$
    – yagoo
    Mar 24, 2015 at 8:42
  • $\begingroup$ So what is $f$? Is it a given function? Please describe all objects that you use. $\endgroup$
    – cfh
    Mar 24, 2015 at 9:14

2 Answers 2


Introduce a new unknown $\phi(s)$ such that $$ (\dot x(s), \dot y(s)) = (\cos\phi(s), \sin\phi(s)), $$ then rewrite the problem as a system of coupled first-order ODEs. You now got rid of the constraint and can use any standard ODE integrator.

  • 1
    $\begingroup$ @yagoo, if this solve your problem why don't you accept the answer? $\endgroup$
    – nicoguaro
    Mar 30, 2015 at 18:24

You call this system an ODE (ordinary differential equation), but this sort of system is actually called an DAE (differential algebraic equation).

What should I do to deal with above constraint?

The constraint $\dot{x}^2 + \dot{y}^2 = 1$ is a smaller problem than the fact that you need to determine the time evolution of $\lambda$, or more precisely analytically compute the constraint satisfied by $\lambda$. One way to do this is to first write the system in integral form:

\begin{equation} \left\{ \begin{aligned} \frac{d}{ds}\big(\dot{x} (f(x,y)-\lambda)\big) &= \frac{\partial f(x,y)}{\partial x}\\ \frac{d}{ds}\big(\dot{y} (f(x,y)-\lambda)\big) &= \frac{\partial f(x,y)}{\partial y}\\ \end{aligned} \right. \end{equation}

and then introduce dummy variables for converting it into semi-explicit form:

\begin{equation} \left\{ \begin{aligned} \frac{dx}{ds} &= \dot{x}\\ \frac{dy}{ds} &= \dot{y}\\ \frac{du}{ds} &= \frac{\partial f(x,y)}{\partial x}\\ \frac{dv}{ds} &= \frac{\partial f(x,y)}{\partial y}\\ \end{aligned} \right. \end{equation} \begin{equation} \left\{ \begin{aligned} u &= \dot{x} (f(x,y)-\lambda)\\ v &= \dot{y} (f(x,y)-\lambda)\\ 1 &= \dot{x}^2 + \dot{y}^2\\ \end{aligned} \right. \end{equation}

You are lucky that the algebraic equation system can be solved for $\dot{x}$, $\dot{y}$, and $\lambda$. In general, you must use Pantelides algorithm (for example) to generate more constraints (and dummy variables) until your original dynamic variables are uniquely determined by the constrains.

You don't even need $\lambda$, so let's eliminate it: \begin{equation} \left\{ \begin{aligned} \dot{y}u &= \dot{x}v\\ 1 &= \dot{x}^2 + \dot{y}^2\\ \end{aligned} \right. \end{equation}

This algebraic equation system now allows you to compute $\dot{x}$ and $\dot{y}$ from $u$ and $v$. So BVP4C will only see $x$, $y$, $u$, and $v$, and you will solve for $\dot{x}$ and $\dot{y}$ yourself and use it where it is needed.

Edit: Warning, the solution above is wrong! Writing the system in integral form is not as straightforward as suggested, because we actually have

\begin{equation} \left\{ \begin{aligned} \frac{d}{ds}\big(\dot{x} f(x,y)\big)-\lambda\frac{d}{ds}\dot{x} &= \frac{\partial f(x,y)}{\partial x}\\ \frac{d}{ds}\big(\dot{y} f(x,y)\big)-\lambda\frac{d}{ds}\dot{y} &= \frac{\partial f(x,y)}{\partial y}\\ \end{aligned} \right. \end{equation}

We could introduce dummy variables $\ddot{x}$ and $\ddot{y}$ to get an integral form

\begin{equation} \left\{ \begin{aligned} \frac{d}{ds}\big(\dot{x} f(x,y)\big) &= \frac{\partial f(x,y)}{\partial x}+\lambda\ddot{x}\\ \frac{d}{ds}\big(\dot{y} f(x,y)\big) &= \frac{\partial f(x,y)}{\partial y}+\lambda\ddot{y}\\ \frac{d}{ds}\dot{x} &= \ddot{x}\\ \frac{d}{ds}\dot{y} &= \ddot{y}\\ \end{aligned} \right. \end{equation}

And now we really need to apply Pantelides algorithm...

  • $\begingroup$ @yagoo I just noticed that my answer is wrong. The transformation to integral form omitted a term involving a derivative of $\lambda$. This probably means that the index of the system is bigger than one, which will require some mastery of the techniques to solve. $\endgroup$ Mar 30, 2015 at 7:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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