グロブ

1.1.7 Newton 法による平方根

問題　1.7

さて、Newton 法 も何となくわかった気になったが、本書に書かれているロジック とは何か違う。

前回は なんや、微分して、x軸との交点を求めて・・・ みたいなイメージだったのだが、本書では、予測値で、平方根を出したい値を割って、さらにその平均をとることを繰り返している。

√2 を計算する

うーん 違和感を感じる。何故平均でいいんだ？

最近、「数学ガール」 を読了して、数式をみて思考停止しなくてすむようになったのと、ちょっと数式を展開するのがおもしろく感じたこともあるし、今年は、中学程度も怪しい数学レベルをせめて高校程度までもっていこうという野望を抱いているので、ここはちょっとこだわっておこう。そのうち息子に算数を教えることができるようにしとくためにも・・・

まず、上記で行っている作業は、平方根を求める対象を b として、予測値を a とし、素直に以下の式とする。

で、一方、Newton 近似の 一般式は、

  でした。

本書で言う、予測値の部分が、 箇所なので、やはり平方根を求める対象を b とすると、

となるので、

2a で通分して、

分子を計算して、

分母のa を消すと・・・

お、一致しましたね。

これでようやく、本書で実装されたところの意味を理解することができました。

一人読書会 「計算機プログラムの構造と解釈」 (2)

1.1.7 Newton 法による平方根

Newton法 により平方根を求める手続きを実装する。

問題　1.7

まず、残念ながら、Newton 法 を知らないので、今年の目標は数学の勉強をすることなので、読み流さずにWikipediaで調べることとする。

重ねて残念ながら、Wikipedia の説明では、知識不足で、わかったような、わからないような。

偶然、所有している本「微分・積分のしくみ」に、ニュートン近似として、丁寧な説明があったので、まずそちらを参考にして、Newton法による平方根の近似とはどんなものかの理解を試みる。

たとえば、√5 の値を調べるとする。

2次関数のグラフを書く

まず、上の図の2次関数のグラフは、y = x^2 – 5 のグラフ。

√5  を 満たすx、　x = √5  は、 x^2 = 5 であり、 x^2 –5 = 0 となる。

Python で、

exp  = lambda x : x**2 - 5

と定義してみる。

適当な開始位置を決める

求めたいのは、上記の2次関数のグラフで、y=0 となるときの　x の値。グラフを見ると、0 ～ 5 の間にあるので、とりあえず開始位置を x = 10 としてみる。すると、x = 10 の場合、y = 95 に値が定まる。(上図①)

接線を求め、x軸との交点を求める

y = x^2 –5 のグラフ上の、上記で決めた (10, 95) での接線が、x1 で、y = 0 となる接線の方程式を求める。

y = f(x) = x^2 – 5  を 微分すると  f'(x) = 2x で、f’(10) = 20 が 接線の傾きになる。

とすると、

接線の方程式 : y – f(a) = f’(a)(x – a)

y – 95 = 20x – 200

y = 20x – 105

( y = 20x + b 、 95 = 200 + b 、　b = –105 )

この接線が、y = 0 となるのは、x = 5.25 となる (上図②)

これで、みてわかるように、最初の値、x = 10 より、y = x^2 – 5 の x軸との交点により近づいた。

ここで取得した値、5.25 を初期値として、上記の処理を繰り返すことで、√5 の値に近似していく。

というのが ニュートン近似だと。

上図の ① を ② を 、として、一般化すると、 となる。

計算結果を出力してみる。

ようするに、n をどんどん増やしていくと、徐々に近づいていっているのがわかる。とりあえず、富士山麓に鸚鵡が鳴く条件は満たした。

グラフのコードなど

Python で グラフを書くための準備については、以下を参照

Python で グラフをかけるようにしとくと、こういった勉強の時に役立つなぁ。と。簡単にかけるし。Python イイ!

http://typea.info/tips/wiki.cgi?page=Python+NumPy

http://typea.info/tips/wiki.cgi?page=Python+matplotlib

上記のコード

# -*- encoding:utf-8 -*-
from pylab import *
from numpy import *
from numpy.ma.core import arange

# 二次関数 f(x) y = x^2 - 5
exp  = lambda x : x**2 - 5
# 二次関数の微分 f'(x) y = x^2 - 5 
expd = lambda x : 2 * x
        
def f1(x_range, intarval):
    '''
        関数expに対して、xの値 x_range[0] ～ x_range[1] 間隔 intarval に対応する y を取得
    '''
    x = arange(x_range[0], x_range[1], intarval)
    y = [exp(tx) for tx in x]
    return (x, y)

if __name__ == '__main__':
    grid(True)
    
    # 二次関数のグラフを書く
    (x, y) = f1((-11.0, 11.0),0.1)
    plot(x, y)

    # 適当な値を設定
    tmp_x2 = 10.0
    
    repeat = 5
    for i in range(repeat):
        print '%d : %s' % (i+1, tmp_x2)
        tmp_x1 = tmp_x2
        
        x1 = (tmp_x2, tmp_x1) 
        y1 = (0, exp(tmp_x1))
        plot(x1, y1, 'k--')
    
        tmp_x2 = tmp_x1 - (exp(tmp_x1) / expd(tmp_x1))
        x1 = (tmp_x2, tmp_x1) 
        y1 = (0, exp(tmp_x1))
        plot(x1, y1)
      
    show()

ようやく、問題に取りかかる準備ができた。

このペースだと読了するのに1年ではきかんな～。

一人読書会 「計算機プログラムの構造と解釈」 (2)

1.1.7 Newton 法による平方根

Newton法 により平方根を求める手続きを実装する。

問題　1.6

以下、特殊形式の if を 通常の手続きである、new-if を定義し、平方根を Newton法 で求める処理を書き直した。

(define (square x) (* x x))
(define (sqrt-itr guess x)
  (new-if (good-enugh? guess x)
      guess
      (sqrt-itr (improve guess x)
                x)))
(define (improve guess x)
  (avarage guess(/ x guess)))
(define (avarage x y)
  (/ (+ x y) 2))
(define (good-enugh? guess x)
  (< (abs (- (square guess) x)) 0.001))
(define (sqrt2 x)
  (sqrt-itr 1.0 x))

(define (new-if predicate then-clause else-clause)
  (cond (predicate then-clause)
        (else else-clause)))

実際に動かしてみると、スタックオーバーフローだろう、Out Of Memory が発生した。

特殊形式 if を使った場合には、想定通りの動きをしていた。ぱっと見、等価に思えるが。

うーーん 。。。 なぜだろう 。。。

・・・

特殊形式・・・ define は 特殊形式で、作用は行わない。。。 if も特殊形式・・・

通常評価は、部分式を評価し、作用(演算子を適用)させる・・・

おお！そうか、特殊形式ではないから、被演算子がすべて評価されてしまうからか！、sqrt-itr の else 節 で、再帰しているのが、終了条件に合致しても、else節 がなんと、評価されてしまうのだ！

if は特殊形式で、条件に合致した then または、 else のどちらかしか評価しないのだ。

なるほど。

一人読書会 「計算機プログラムの構造と解釈」 (1)

1.1.3 組合せの評価

組合せは、以下のように評価される.

1. 組合せの部分式を評価
2. 最左の演算子を被演算子に作用させる

これを 作用的順序の評価(applicative-order evaluation) というが、もう一つ、組合せを完全に展開してから演算を行う評価方法を、正規順序の評価(normal-order evaluation) という.

問題　1.5

>(define (p) (p))
>(define (test x y)
  (if (= x 0)
      0
      y))
>(test 0 (p))

if がショートサーキットであるとするならば、解釈系が、作用的順序であれば、(test 0 (p)) で、最初に x = 0 が評価されて、結果が 0 となり、test の第2引数 (p) は評価されない。

正規順序の評価であれば、まず組合せが展開される段階で、test の第2引数 (p) が無限ループ(?) に陥ってしまう。

・・・ のではないか。

ちなみに、DrScheme で実行すると、結果が返ってこなかった。

「計算機プログラムの構造と解釈」 を購入。しっかりと読み込みたいので、1年くらいかけて一人読書会としゃれ込もうと思う。

・・・電車で読むだけではもったいない。

1.手続きによる抽象の構築

まず、全編を通して、プログラムの説明には、C や Java や Python ではなく、Lisp の方言である、Scheme (スキーム) を使用する。

そもそも、Scheme の設計者である、ジェラルド・ジェイ・サスマン が、本書の著者。

amazon のレビューをみると、訳に対する苦情が多いようだが、自分も最初のページで違和感を感じた。デバッグのことを「虫とり」とか。

ただ、訳者が、1931年生まれの先生ということを知ると、訳の雰囲気の違和感も哲学書のような雰囲気を醸し出してきた。自分はこういう訳も嫌いではない。

で、本書では、 cheme がわからないと話にならないので、最初の章は、文法および考え方の説明に費やされている。

こちらに簡単にまとめた。さらにまとめると、

確認用の処理系には、DrScheme を使用することとする。

式

手続きの作用を表現。括弧で囲んで、組合せと呼ぶ。

> (+ 1 2)
3

この例だと、加算(+) が、手続きを表す演算子(オペレーター)、1、2 が作用される被演算子(オペランド)

演算子を左端に書く記法を、前置記法という。利点として、任意個数の引数をとることができる、組合せを入れ子にできる。

特殊形式(special forms)

引数に作用しない場合、組合せとは呼ばず、特殊形式という。たとえば、define

> (define x 2)

これで、変数 x に 2 をセットする。

合成手続き

define をつかって、演算に名前をつけることができる。自乗を行う演算に名前をつけてみる。

> (define (square x) (* x x))
> (square 2)
4

条件式と述語

cond ・・・ 場合分け

if ・・・ 場合が2つの場合

論理演算子

and、or、not

用例は、こちら を参照。

問題

問題1.1

入力して確認。

問題1.2

これでよいと思うけど・・・

> (/ (+ 5 4 (- 2 (- 3 (+ 6 (/ 4 5))))) 
     (* 3 (- 6 2) (- 2 7)))
-37/150

問題1.3

いままで経験した言語との考え方の違いを思い知る。

今の知識、思考方法では、こんな感じにしかかけない・・・

変数や配列を使わずに、2番目に大きい値を取得するって困難。

> (define (smaller x y) (if (< x y) x y))
> (define (last x y z) (smaller x (smaller y z)))
> (define (square x) (* x x))
> (define (q1_3 x y z) (+ 
                        (square (if (> x (last x y z)) x (last x y z)))
                        (square (if (> y (last x y z)) y (last x y z)))))
> (q1_3 4 3 2)
25

問題1.4

引数 b が 負の値の場合、符号を反転し絶対値を取得し、a と加算する

おーなんと、被演算子に演算子を置くと、その演算子が作用するのか！

> (define (a-plus-abs-b a b)
    ((if (> b 0) + -) a b))
> (a-plus-abs-b 2 -3)
5

おー結構しんどいなー。今日はここまで。

ゆっくりでもいいから何とか続けていきたい。。。

先週、毎年恒例の麻雀旅行に有志と行ってきた。4回目くらい。

いままでは温泉にでもいって本当に麻雀やっておしまいだったが、前回金沢で兼六園いったらやはり観光すると非常におもしろかったので、今回は結構観光にウェイトをおいた日程とした。

せっかくなので、簡単に旅行記などを。

といっても、そもそもカメラ自体持って行ってないし、携帯の充電器ももってきてないので、携帯のカメラもあまり使えない。もう少し考えればよかったよ。

大神神社

まず、大神神社。

個人的には神社巡りが好きで、家族旅行ではここのところ、伊勢神宮、諏訪大社、出雲大社と神社づいているのだが、次は大神神社行きたいねぇと嫁といっていたのだが、抜け駆け。

熱田神宮にはよく行くのだが、それに比べると、プリミティブという感じの参道。

参道を登り切ると社が。その後ろにご神体である三輪山がそびえている。というか、自分のイメージだと鳥居と三輪山しかないかと早合点していたのだが、ちゃんと社があるのね。

お供えに卵がいっぱいあるので、なんでかなぁとおもったら、そうか卵は蛇の好物だからね。納得

纒向遺跡

大神神社のすぐ近くに纒向遺跡ってのがある。ヤマト政権発祥の地、邪馬台国の有力候補地、前方後円墳発祥の地。

遺跡といっても、町中に分散しているというか、町自体が遺跡というか、そこら中に前方後円墳が・・・

奈良すごい。時間の都合上、箸墓遺跡を遠目に見るくらいしかできなかった。是非もう一度きたいところ。

非常にロマンがかき立てられる地だ。大神神社のすぐ隣に、埋蔵文化財センターという発掘されたものを展示しているところがあり、そこでにわか考古学談義。みんな中途半端な知識であーでもないこーでもないと話し合う。こんなすごい遺跡なのに客が我々の他には怪しげなおじさん一人。

後ほど京都へいくのだが、京都にくらべて素朴で奈良はいいな。学生時代に大阪にすんでたので、もっときとけばよかった。4年間で数回しかきてなかったが残念。

石舞台

まぁ今回は自分の好みで行き先を決めた感があり申し訳なかったが、石舞台。

息子の名前に一文字いただいている手塚治虫大先生の代表作の一つである火の鳥のヤマト編に登場する。みてみたいと思い続けておおかた30年。

やっとみることができました。

なんと、この巨石の下が、石室になっていて入れるのですよ。しらんかった。

それにしてもこれほどの名所の割には人が少ない。修学旅行でないとこんな感じか。そこがいい！

これまた手塚治虫大先生の代表作の一つ「三ツ目がとおる」に出てきた酒船石がすぐ近所にあるのだが、時間の都合上よることができなかった。残念。

高松塚古墳

石舞台からちょっと行くと、飛鳥美人壁画の高松塚古墳。現状は、掘り起こした後に再現した塚となっている。

ここでも、園内に、みかんが無造作におかれており200円と書かれている。すてきだ。

現物は、修復中なので、現地には壁画を再現した高松塚古墳壁画館 があるが、これがかなり圧巻。ここでも古代人に思いを馳せて、各自適当なことを言いまくる。

すぐそばには、文武天皇陵。森の中になにかクリーチャーがいても不思議ではなさそう。

ここらで奈良タイムアップ。宿が京都のため京都へ向かわねば。

奈良主体で宿をとってもよかったなぁというのは、みんなの感想。

毎回、おっさんの中の唯一の若者(といっても三十路近いが)が宿の手配やら、車の運転やらしてくれるので、この場を借りて感謝。

京都旅行記は、次回気が向いたら書こう。