Lets say I have a system of differential equations which has the form $$\dot{C}_{\alpha,\beta,m} = f_{\alpha,\beta,m}(C_{\alpha,\beta,1},\ldots,C_{\alpha,\beta,N};t).$$ The $f$s are some functions of time and of the $C$s. The exact form is not important.

In principle, this can be solved with a simple RK4 algorithm. However, there is now the additional (orthogonality) constraint $$\sum_{\alpha,\beta}C_{\alpha,\beta,m}^* C_{\alpha,\beta,n} \stackrel{!}{=} \delta_{m,n},$$ which should hold at any time $t$. The star denotes the complex conjugate.

My question is now, how this constraint can be implemented in my integration scheme/RK 4 algorithm.


1 Answer 1


In general, one cannot expect rk4 to maintain quadratic invariants of the system, it simply doesn't do that. Methods that do maintain specific invariants have to be specially devised — this usually goes by the name of geometric integration. Symplectic integrators are the most common kind of these methods, although not for your problem.

See, for example, the book Geometric Numerical Integration by Hairer, Wanner, Lubich, or the notes Solving Differential Equations on Manifolds by Hairer (free online).

Probably the simplest way to incorporate the constraint is to project the approximate solution onto the constraint manifold ("Projection methods", Example 1.3, p. 21 in the Hairer notes above). Given a system of ODEs $$ \dot X = f(t, X), \qquad X^*X=I, $$ define $\pi(X)$ to be the closest unitary matrix to $X$. If $X = UH$ is the polar decomposition of $X$, then $\pi(X) = U$.

The simplest is to project onto the constraint manifold after each RK step, $\tilde X_{n+1} = \pi(X_{n+1})$. The second simplest is to rewrite the ODE as $$ \dot X = f(t, \pi(X)), $$ which will incorporate the projection into the RK substeps. (This doesn't actually modify the ODE because $\pi(X)=X$ for any exact solution.)

Beyond this, the constraint manifold is called a Stiefel manifold, and there already seem to be some number of articles on google scholar on numerical integration on Stiefel manifolds.

  • $\begingroup$ Thanks for your answer. Is there also a way using Lagrange multipliers? I think this is used in the Car-Parrinello method: en.wikipedia.org/wiki/Car%E2%80%93Parrinello_method Do you now something about this way of treating the problem? $\endgroup$
    – Merlin1896
    Feb 3, 2016 at 8:19
  • $\begingroup$ @Merlin1896 I don't know. $\endgroup$
    – Kirill
    Feb 3, 2016 at 13:06

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.