# How to solve ODEs with constraints using BVP4C?

I am using BVP4C to solve a system of ODEs which is as follows. \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.

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?

• 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? Mar 23 '15 at 17:18
• @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$? Mar 24 '15 at 7:56
• @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. Mar 24 '15 at 8:42
• So what is $f$? Is it a given function? Please describe all objects that you use.
– cfh
Mar 24 '15 at 9:14

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.

• @yagoo, if this solve your problem why don't you accept the answer? Mar 30 '15 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:

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

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

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

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: \left\{ \begin{aligned} \dot{y}u &= \dot{x}v\\ 1 &= \dot{x}^2 + \dot{y}^2\\ \end{aligned} \right.

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

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

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

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

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

• @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. Mar 30 '15 at 7:45