newton.scm

;; f(z) = z^3 + 2z^2 + 1
(define (f z) (+ (expt z 3) (* 2 (square z)) 1))

;; f'(z) = 3z^2 + 4z
(define (df z) (+ (* 3 (square z)) (* 4 z)))

;; ニュートン法。初期値zから出発。変化が充分少なくなったらその時のzがf(z)=0の解。
;; 微分が0に近くなった場合はエラーとする。
;; 循環してしまうケースは検出しない (止まらなくなる)
(define (newton z)
  (let ([fz  (f z)]
        [dfz (df z)])
    (when (< (magnitude dfz) 1.0e-15)
      (error "Unstable"))
    (if (< (magnitude (/. fz dfz)) 1.0e-15)
      z
      (newton (- z (/. fz dfz))))))

#|
;; 実行例。

gosh> (newton 1)
-2.2055694304005904
gosh> (newton 1+i)
0.10278471520029513+0.6654569511528134i
gosh> (newton -1+i)
0.10278471520029517+0.6654569511528134i

;; 解であることの確認 (誤差があるので完全に0にはならない)
gosh> (f -2.2055694304005904)
0.0
gosh> (f 0.10278471520029513+0.6654569511528134i)
1.1102230246251565e-16-1.1102230246251565e-16i
gosh> (f 0.10278471520029517+0.6654569511528134i)
1.1102230246251565e-16

|#