I am trying to force a constraint on a ODE. I have the following equation:
$$ \begin{aligned} \dot x &= -0.04 x + 10 y z \\ \dot y &= 0.04 x - 10 y z - 3000 y^2 \\ \dot z &= 3000 y^2 \end{aligned} $$
This equation was derived assuming a continuity equation where $x+y+z=1$
Instead of solving for $\dot z$ I would like to solve the first two terms and apply the continuity constraint.
I think this should be done with a Lagrangian multiplier technique. However, I do not know how to derive the expression. Can anyone give me a hand with the derivation?
I know that for a function $f\left(x,y\right)$ subjected to a constraint $g\left(x,y \right)$, the lagrangian multiplier expression will be:
$$ \left\{\begin{matrix}\nabla f\left(x,y\right)=-\lambda\nabla g\left(x,y\right)\\g\left(x,y\right)=0\end{matrix}\right. $$
Will this still hold for constraints in ODE's?
The idea of first order optimality is that the directional derivative of the objective function should be zero in all directions in which you are permitted to move. So in the unconstrained setting, this means $\nabla f$ has to be perpendicular to everything, and thus must be zero. In the setting where all constraints are equality constraints, $\nabla f$ has to be perpendicular to any vector tangent to the constraint surface, and so must be a linear combination of the gradient(s) of the constraint function(s).
The same idea applies in the setting of constrained ODEs, but the result looks a bit different. The idea in the ODE setting is that the directional derivative of the constraint function should be zero in the direction of the flow. So for $g(x,y,z)=x+y+z$ you should have $(\dot{x},\dot{y},\dot{z})$ perpendicular to $\nabla g$. From there, given $\dot{x}$ and $\dot{y}$ you can solve for $\dot{z}$. If $g$ were more complicated, then this might be worth doing dynamically as you go through some kind of numerical simulation. But because $g$ is really simple in this problem, it is easier to just write $z=C-x-y$ where $C=x_0+y_0+z_0$ and substitute that into the equations for $\dot{x}$ and $\dot{y}$. That yields
$$\dot{x}=-0.04x+10y(C-x-y) \\ \dot{y}=0.04x-10y(C-x-y)-3000y^2 \\ z=C-x-y.$$