source: CGBLisp/examples/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
;;
;;Define the map f: (z,c)->(z^2+c,c)
(setf f (cdr (string-read-poly "[z^2+c,c]" '(z c))))
;;
;;Define the identity map as a polynomial
(setf id (cdr (string-read-poly "[z,c]" '(z c))))
;;
;;Define a constant polynomial 1 in variables z and c
(setf one (string-read-poly "1" '(z c)))
;;
;;(f-composition n) returns f o g o ... o f (n-times) as polynomial
(defun f-composition (n) (poly-dynamic-power f n))
;;
;; g = f^n-id
(defun g (n) (car (mapcar #'poly- (f-composition n) id)))
;;
;; (f^n)' (derivative over z = 0-th variable)
(defun df (n) (car (partial (f-composition n) 0)))
;;
;; Flip bifurcations occur when derivative is -1 at some fixed point
;; (ideal n) is the ideal spanned by f^n-id and f'+1 and its zeros
;; are clearly the locus of flip bifurcations
(defun ideal (n) (list (g n) (poly+ (df n) one)))
;;
;; Printer of the n-th ideal
(defun print-ideal (n) (poly-print (cons '[ (ideal n)) '(z c)) (terpri))
;;
;; Eliminate z from the equations, because we are just after the values of
;; the parameter c
(defun bifurcation (n)
  (mapcar #'poly-contract (elimination-ideal (ideal n) 1)))
;;
;; The same, but elimination done using resultant
(defun resultant-bifurcation (n)
  (poly-contract (poly-resultant (g n) (poly+ (df n) one))))
;;
;; Print the polynomial whose zeros are the values of c
;; for which flip bifurcation occurrs
(defun print-bifurcation (n)
  (poly-print (cons '[ (bifurcation n)) '(c))
  (terpri))
;;
;; Same for resultant version
(defun print-resultant-bifurcation (n)
  (poly-print (resultant-bifurcation n) '(c))
  (terpri))