I'm trying to solve an ODE problem. Let's say $\mathbf{x}(t)$ represents the position of a particle at time $t$, and $\mathbf{u}(\mathbf{x},t)$ is a velocity field defined in Cartesian coordinates on some manifold $\mathbb{M} \subset \mathbb{R}^3$ represented as a level set of the form $F(\mathbf{x}) = 0$. Then, the path taken by the particle can be given by: \begin{align} \frac{d \mathbf{x}}{d t} &= \mathbf{u}(\mathbf{x},t), \\ F(\mathbf{x}(t)) &= 0,\\ \mathbf{u}(\mathbf{x},0) &= u_0. \end{align} How do I go about solving this constrained IVP? Assume I can evaluate $F$ and $\mathbf{u}$ analytically. Are there Runge-Kutta or linear multistep methods designed for this sort of ODE with a nonlinear solution constraint?

  • 1
    $\begingroup$ I'm not an expert on the field, but does that even have a solution? Without the constraint with $F$ that should already have a unique solution. If that solution does not satisfy $F(x) = 0$... $\endgroup$ Mar 22, 2020 at 20:35
  • $\begingroup$ I'm pretty sure that given the right constraint on $\mathbf{u}$, this should have a solution. After all, it's simply an ODE saying a particle is moving on a surface in some given velocity field. $\endgroup$ Mar 22, 2020 at 20:45
  • $\begingroup$ It's not an ODE I think, It's a PDE in fact, why? See here: $$\frac{d \mathbf{x}}{dt} = \mathbf{u} (\mathbf{x},t)$$ so: $$\frac{d \mathbf{u}}{dt} = \mathcal{F}(\mathbf{x},t)$$ where $\mathcal{F}(\mathbf{x},t)$ is the force field acting on your particle. You must know this force field, otherwise it's not possible to solve this problem. But you have: $$\frac{d \mathbf{u}}{d t} = \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u}$$ , so finally: $$\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = \mathcal{F}(\mathbf{x},t)$$ So it's a PDE... $\endgroup$ Mar 22, 2020 at 21:14
  • $\begingroup$ I already said I know $\mathbf{u}$ and $F$ analytically. However, the ODE system is in the Lagrangian frame, and there are numerical benefits to working in that frame. $\endgroup$ Mar 22, 2020 at 22:33
  • 3
    $\begingroup$ What you have is what's called a Differential-Algebraic equation (DAE). The form of $\mathbf u$ is then typically $\mathbf u(t)=\mathbf w(t) + \mathbf \lambda(t)$ where $\mathbf w$ is known, and $\lambda$ is a Lagrange multiplier corresponding to the constraint. The (time-dependent) Stokes equations are an example, as are many equations in mechanics. There are specialized solvers for DAEs -- which you should be able to find now that you know the correct term :-) $\endgroup$ Mar 23, 2020 at 0:36

1 Answer 1


As @WolfgangBangerth already commented, this is usually referred to as a differential algebraic equation (DAE). These have their own challenges and there are special numerical methods for them.

In general, for a system of DAEs you might have fewer differential equations than unknowns, so the algebraic constraint(s) serve to close the system and identify a unique solution. It seems that in your case you have an ODE for each unknown, so presumably the exact solution of your system of ODEs satisfies the constraint exactly. Your numerical solution will not, and it sounds like it is important that the constraint be satisfied exactly (or to machine precision). In this setting it is not always best to use the formalism of DAEs (though you certainly can). If your constraint is low-dimensional -- especially if you have just a single constraint -- then there are perhaps more useful methods coming from the field of geometric numerical integration. This deals for instance with systems where it is important to conserve energy or momentum. A common technique is to use orthogonal projection onto the constraint manifold after each step, but another approach I've developed recently is that of relaxation, which requires a very small and simple modification to whatever numerical method you are already using.

  • $\begingroup$ Brilliant! I've studied DAEs (years and years ago), but I seem to have forgotten about them. This is great. The constraint manifold projection trick is something I've used for linear constraints, and I think I can figure out how to do it for nonlinear ones too. I'll take a look at the relaxation paper. Thank you! $\endgroup$ Mar 24, 2020 at 1:36

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.