# Integration of differential equation with orthogonality constraint

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.

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