source: CGBLisp/examples/henon-bifurcation.lisp@ 1

;;
;; A computation of the places at which flip bifurcation occurrs in the
;; quadratic family
;; Note: period 3 finished quite easily, period 4 runs forever
;; Two methods: grobner and resultant
;;


(setf order (elimination-order
             2
             :primary-order #'grevlex>
             :secondary-order #'grevlex>
             ))

(setf f (cdr (string-read-poly "[a-x^2+b*y,x,a,b]" '(x y a b) :order order)))
(setf id (cdr (string-read-poly "[x,y,a,b]" '(x y a b) :order order)))
(setf one (string-read-poly "1" '(x y a b) :order order))

(defun f-composition (n) (poly-dynamic-power f n order))

;; Flip bifurcations occur when derivative is -1 at some fixed point
;; g
(defun g (n) (subseq (mapcar #'(lambda (x y) (poly- x y order))
                             (f-composition n) id) 0 2))
(defun f-jacobian (n) (characteristic-polynomial (jacobi-matrix (f-composition n) 2 2)
                                                 order))
(defun flip-value (n) (poly-scalar-composition
                       (f-jacobian n)
                       (cdr (string-read-poly "[x,y,a,b,-1]" '(x y a b) :order order))
                       order))


(defun ideal (n) (cons (flip-value n) (g n)))

(defun print-ideal (n) (poly-print (cons '[ (ideal n)) '(x y a b)) (terpri))

(defun bifurcation (n)
  (mapcar #'(lambda (x) (poly-contract x 2))
          (elimination-ideal (ideal n) 2
                             :order order
                             )))

(defun print-bifurcation (n)
  (poly-print (cons '[ (bifurcation n)) '(a b))
  (terpri))