Implementation of Newton Raphson algorithm in clojure

newton-raphson.clj

(ns newton-raphson.core)

;;utility function for list manipulation
(defn keep-list [expr]
  (if (and (coll? expr) (not (= 1 (count expr)))) (list expr) expr))


(defn to-list [expr]
  (cond (list? expr) expr
(coll? expr) (seq expr)
:default (list expr)))

(defn replace-in-list [coll old-val new-val]
  (map #(cond 
         (coll? %) (replace-in-list % old-val new-val)
         (= % old-val) new-val  
         :default %) 
       coll))
  

(defn derivative-prod [expr variable]
  (conj 
   (let [elts (rest expr)]
     (loop [members elts
            res ()
            tmp (first elts)
            previous ()
            next (rest elts)]
       (cond (= () members) res
             :default
             (recur
              (rest members)
              (conj res (concat '(*) previous (keep-list (to-list (derivative tmp variable))) next)) ;;inelegant
              (second members)
              (conj previous tmp)
              (rest next)))))
   '+))

;;messy implementation of mathematical derivative
(defn derivative [expr variable]
  ;;not simplified derivative calculations
  (cond (number? expr) 0;;(list 0)
        (symbol? expr) (if (= variable expr) 1 0);;(list 1) (list 0))
        :default
        (let [op (first expr)
              args (rest expr)
              der #(derivative % variable)]
          (cond (= '+ op) (conj (map #(derivative % variable) (rest expr)) '+)
                (= '* op) (derivative-prod expr variable)
                (= '- op) (conj (map #(derivative % variable) (rest expr)) '-)
                ;;(= '- op) (list '- (derivative (second expr) variable) (map #(derivative % variable) (rest (rest expr))))
                (= '/ op) (if (= 2 (count args)) (let [u (first args)
                                                       v (second args)]
                                                   (list '/ (list '- (list '* (der u) v) (list '* u (der v))) (list '* v v)))
                              (der (list '/ (first args) (concat '(*) (rest args)))))
                (= 'sin (first expr)) (list '* (derivative (second expr) variable) (list 'cos (second expr)))
                (= 'cos op) (list '* -1 (der (second expr)) (list 'sin (second expr)))
                (= 'exp op) (list '* (der (first args)) (list 'exp (first args)))
                :default (do (println (first expr) "not handled yet...") nil)))))



(defn my-eval [expr variable value]
  (let [new-expr (replace-in-list expr variable value)]
    (eval new-expr)))

(defn abs [x] 
  (if (< x 0) (- x) x))


(defn to-bigdec [x] 
  (if (ratio? x) (/ (bigdec (numerator x)) (bigdec (denominator x))) (bigdec x)))


;;implementation of newton-raphson algorithm
(defn nr-solve [expr variable start tol]
  "solve (expr = 0) as an equation in *variable*"
  (let [der (derivative expr variable)
        run-algo 
        (fn [s] 
          (let [slope (my-eval der variable s)
                k (my-eval expr variable s)
                candidate (- s (/ k slope))]
            (if (< (abs (my-eval expr variable candidate)) tol) 
              candidate
              (recur  candidate ))))
        sol (run-algo start)]
    sol))