SICP 問題 1.36(不動点を求める方法を使ってx^x=1000の解を見つけ)

【問題】

問題1.22で示した基本の newline と display を使い、生成する近似値を順に印字するようfixed-pointを修正せよ。
次に、x├→ log(1000)/log(x) の不動点を探索することで、x^x = 1000 の解を見つけよ。
(自然対数を計算するschemeの基本log手続きを使う。)
平均緩和を使った時と使わない時のステップ数を比べよ。
(fixed-pointの予測値を 1 にして始めてはいけない。log(1) = 0 による除算を惹き起こすからだ。)

【解答】

求めた近似値を画面に印字するには、f を作用させた直後の値を表示すればいい。
なので、letで (f guess) を next 束縛した直後に表示すればいいと思う。

;精度?
(define tolerance 0.00001)

;不動点を求める手続き
(define (fixed-point f first-guess)
  ;十分近いかどうか判定する手続き
  (define (close-enough? v1 v2)
    (< (abs (- v1 v2) ) tolerance))
  ;不動点を求める核となるロジック。再帰します。
  (define (try guess)
    (let ((next (f guess)))
      (display next) ;ここで画面に印字。
      (newline)      ;そして改行。
      (if (close-enough? guess next)
	  next
	  (try next))))
  ;実際の処理を開始
  (try first-guess))

さて、次。


x├→ log(1000)/log(x)
をlambda式にするとこんな感じ。

(define (f x)
  (/ (log 1000) (log x)))

まずはこれで計算してみよう。


gosh> (fixed-point f 2)
9.965784284662087
3.004472209841214
6.279195757507157
3.759850702401539
5.215843784925895
4.182207192401397
4.8277650983445906
4.387593384662677
4.671250085763899
4.481403616895052
4.6053657460929
4.5230849678718865
4.577114682047341
4.541382480151454
4.564903245230833
4.549372679303342
4.559606491913287
4.552853875788271
4.557305529748263
4.554369064436181
4.556305311532999
4.555028263573554
4.555870396702851
4.555315001192079
4.5556812635433275
4.555439715736846
4.555599009998291
4.555493957531389
4.555563237292884
4.555517548417651
4.555547679306398
4.555527808516254
4.555540912917957
4.555532270803653
4.555532270803653
gosh>
検算してみると・・・

gosh> (/ (log 1000) (log 4.555532270803653))
4.555537970114198
gosh>
あってるねぇ。


さて、最後。
平均緩和について説明をしないと問題が解けないわけだが、しかしどーやってまとめたらいいのやら。。
この概念が必要な理由として、「振動が発生して近似値へ収束しないケース」が参考書に挙げられていた。

例えば y = x/y という方程式の場合。

最初の予測値を y_1 とすると、次の予測値 y_2 は y_2 = x/y_1、
じゃ、その次の予測値は、と考えると、 y_3 = x/y_2 = x/(x/y_1) = y_1 となる。
つまり、交互に、y_1, x/y_1, y_1, x/y_1, ...となり、いつまでたっても収束しないわけだ。
じゃ、どうするか。そこで出てくるテクニックがこれ。


y = x/y
y + y = x/y + y (両辺に y を足す)
2y = y + x/y
y = (y + x/y)/2
これで再度計算すると振動しないらしい。理屈は・・・ちょっとわからんかった。
というわけで、定義した f をこういう風に修正してみよう。

(define (f x)
  (/ (+ x (/ (log 1000) (log x)))
     2))

こいつを読み込ませて、もう一度実行してみた結果がこれ。


gosh> (fixed-point f 2)
5.9828921423310435
4.922168721308343
4.628224318195455
4.568346513136242
4.5577305909237005
4.555909809045131
4.555599411610624
4.5555465521473675
4.555537551999825
4.555537551999825
gosh>
ステップ数が減ってる。。う〜む。。

SICP § 1.3.3 一般的方法としての手続き(区間二分法による方程式の零点の探索)

あんまり意味ないかもしれないけど、一応学習したのでもったいないから晒す。

区間二分法による方程式の零点の探索】

二点 a, b が f(a) < 0 < f(b) なら、f は a と b の間に少なくとも一つの零点を持つ。
零点をみつけるのに x を a と b の平均とし、f(x) を計算する。
f(x) > 0 なら f の零点は a と x の間にある。
f(x) < 0 なら f の零点は x と b の間にある。
これを続け、f の零点のあるべき次第に狭い区間を見つけていく。

(define (search f neg-point pos-point)
  ;中間点を取得
  (let ((midpoint (average neg-point pos-point)))
    ;「十分近い」かどうか判定
    (if (close-enough? neg-point pos-point)
	midpoint
	(let ((test-value (f midpoint)))
	  ;区間二分法の核となるロジック
	  (cond ((positive? test-value)
		 (search f neg-point midpoint))
		((negative? test-value)
		 (search f midpoint pos-point))
		(else midpoint))))))

他に必要な手続きも書いとく。

;十分近いかどうかの判定手続きを記述。
(define (close-enough? x y)
  (< (abs (- x y)) 0.001))

;平均値の算出手続き
(define (average . options)
  (let loop ((values options)
	     (accum 0)
	     (count 0))
    (if (null? values)
	(if (= count 0) #f (exact->inexact (/ accum count)))
	(loop (cdr values) (+ accum (car values)) (+ 1 count)))))

指定された関数fが、a, b の位置で未定義な可能性や、(f a)、(f b)が同一符号の場合等は、
この区間二分法そのものが使用できないので、それらのケースを回避するラッパー手続きを定義。

;search手続きを使う前提をチェックする手続き。
(define (half-interval-method f a b)
  (let ((a-value (f a))
	(b-value (f b)))
    (cond ((and (negative? a-value) (positive? b-value))
	   (search f a b))
	  ((and (negative? b-value) (positive? a-value))
	   (search f b a))
	  ;(f a)と(f b)の値が互いに同符号だった場合はエラー
	  (else
	   (error "Value are not of oppoiste sign." a b)))))

これらを全て読み込ませて、


f = sin (※2.0〜4.0の区間)
f = x^3 - 2x - 3 (※1.0〜2.0の区間
の二つの零点を探すとこんな感じ。

gosh> (half-interval-method sin 2.0 4.0)
3.14111328125
gosh> (half-interval-method (lambda (x) (- (* x x x) (* 2 x) 3)) 1.0 2.0)
1.89306640625
gosh>

SICP 問題 1.35(不動点を求める方法を使って黄金律を算出する)

【問題】

(1.2.2節の)黄金比Φが、変換 x ├→ 1 + 1/x の不動点であることを示し、
この事実を使い、fixed-point 手続きによりΦを計算せよ。

【解答】

まずこの問題を解くに当たって必要な情報を記載する。


黄金比
Φ = (1 + √5) / 2
Φ^2 = Φ + 1

「├→」
写像する」と読む。
この表記は、lambdaの数学者流の記法。
「y ├→ x/y」は、「(lambda (y) (/ x y))」のことで、つまり y での値が x/y の関数のことである。

不動点
x が f(x) = x を満たす時、x を関数 f の不動点(fixed point)という。
ある関数 f について、最初の予測値から始め、f を値があまり変わらなくなるまで繰り返し作用させることで、不動点を見つけることができる。

もひとつ、使うことが前提になっている fixed-point 手続きを記載しよう。

;精度?
(define tolerance 0.00001)

;不動点を求める手続き
(define (fixed-point f first-guess)
  ;十分近いかどうか判定する手続き
  (define (close-enough? v1 v2)
    (< (abs (- v1 v2) ) tolerance))
  ;不動点を求める核となるロジック。再帰します。
  (define (try guess)
    (let ((next (f guess)))
      (if (close-enough? guess next)
	  next
	  (try next))))
  ;実際の処理を開始
  (try first-guess))

以上が材料。じゃ、実際に問題を解いていこう〜♪

黄金比Φの関係式から、


Φ^2 = Φ + 1
→ Φ = 1 + 1/Φ (但し、Φ ≠ 0)
これをlambda式で記述するとこんな感じ。

(define fai (lambda (x) (+ 1 (/ 1 x))))

さて、実際にΦの不動点をfixed-pointを使って求めてみる。


gosh> (fixed-point fai 1.6)
1.6180327868852458
gosh>
いーんでない?