Implement bifurcation tracking
Tracking of bifurcations can be achieved by building augmented/extended systems around the existing equation (see e.g. here).
With the generalized Equation approach, it should now be simple to implement such augmented systems.
Basically, two constraints are added to the system:
- The bifurcation's null-eigenvector y is orthogonal to the linearized system: Jacobian * y = 0 (N equations, y are unknowns)
- The bifurcation's null-eigenvector is normalized y*y_original - 1 = 0 (1 equation, the first continuation parameter p1 is the corresponding unknown, y_original is given)
This leads to a system of (2N+1) equations with the unknowns (u, y, p1) = original unknowns, null-eigenvector and first continuation parameter.
Upon initialization of the system, y_original is either passed as the eigenvector of the eigenvalue that crossed zero or generated by solving Jac * y_original = d(rhs)/d(p2), where p2 is the second continuation parameter.