やりました

1.3.3 汎用手法としての手続き

区間二分法によって方程式の根を求める

　区間二分法についてはグラフを適当に書いてみるとわかりやすい。ここのように、プログラミングを用いた方程式の解の導出の題材として多く使われるようだ。

search を直接使うと、ちょっと面倒なことがあります。うっかり f の値が符号の条件を満たさないような点を与えてしまうかもしれないということです。その場合、間違った答えになってしまいます。そうする代わりに、以下に示す手続きを経由して search を使うことにしましょう。これは、どちらの端点が負の関数値を持ち、どちらが正の関数値を持つかを調べ、それに応じて search 手続きを呼ぶというものです。

(Page 71~72).

　起こりうるミスを未然に防ぐための関数。こういう関数を表す名前があると良いと思った。

関数の不動点を求める

関数を受け取ってそれの不動点を求める関数 (この本における fixed-point という関数) は不動点コンビネータと呼ばれるものである。

すなわち高階関数g が不動点コンビネータであるとは、任意の関数 f に対し、p = g(f) とすると, f(p) = p が成立する参考:wiki

　方程式を解くのにも応用できる。例えば単純な三角関数 sin(x) の解を求めてみる。既に私達は ±nπ (n ⊂ N) が解であることを知っているので、結果がこれの何かになればよい。

　sin(x) が 0 となるところが解であるので sin(x) = 0. 次にこの式の両辺に x を加えて sin(x) + x = x のようにすることで不動点が求められるかどうかを調べることができる形になる。というのはこの式の左辺の sin(x) が 0 になるような x 、つまり sin(x) の解であるときに限り x = x という関数と見ることができ、それの不動点は明かに x となる。以下で fixed-point を使って実際に調べている (また毎回の関数適用時の値を出力している)。

(fixed-point (^[x] (+ (sin x) x)) 1.0)
#?-    1.8414709848078965
#?-    2.80506170934973
#?-    3.135276332899716
#?-    3.141592611590653
#?-    3.141592653589793
3.141592653589793


gosh> (fixed-point (^[x] (+ (sin x) x)) -1.0)
#?-    -1.8414709848078965
#?-    -2.80506170934973
#?-    -3.135276332899716
#?-    -3.141592611590653
#?-    -3.141592653589793
-3.141592653589793

1.0 が初期値である例の実際の計算は以下のようになってる

f(x) = sin(x) + x
のとき
x = 1.0
f(x) = 0.8414709848078965 + 1.0                 = 1.8414709848078965

x = 1.8414709848078965
f(x) = 0.9635907245418334 + 1.8414709848078965  = 2.80506170934973

x = 2.80506170934973
f(x) = 0.3302146235499861 + 2.80506170934973    = 3.135276332899716

x = 3.135276332899716
f(x) = 0.006316278690937182 + 3.135276332899716 = 3.141592611590653

x = 3.141592611590653
f(x) = 4.19991402881363e-8 + 3.141592611590653  = 3.141592653589793

このように不動点が求められ、同時に sin(x) の一つの解の近似が (誤差は tolerance の値に依存) 3.141592653589793 であることがわかった。

　しかし不動点を求めるにあたって、初期値 (fixed-point における first-guess) をどのように決めるかが問題になる。本の中ではこれに関して言及されていないため謎だった。今回の場合は sin(x) が周期関数であるため多くの入力に対して解をみつけることができるが、単純な方程式 (例えば x^2 + 5x + 6 のようなもの) は入力に対して出力が大きすぎる上に解が -2, -3 だけであるため収束するのが難しい。

　上記の対策として挙げられているのが平均緩和法であるが、結局のところ良い入力というのが何であるかの判断というのは、関数が複雑になるにつれて非常に難しくなってしまう。初期値を決める問題はアルゴリズムなどのパラメータの決定と同じような問題になるのだろうか。

練習問題 1.35

黄金比 φ(1.2.2 節) が x → 1 + 1/x という変形の不動点であることを示し、そのことを利用して φ を fixed-point 手続きによって求めよ。

(Page 74).

　まず φ = (1 + √5) / 2 が x → 1 + 1/x の不動点であることを示す。

f:id:lhcpr:20170705120257p:plain:w300

次に fixed-point へは単に 1 + 1/x を与えればよい。

gosh> (fixed-point (^[x] (+ 1 (/ 1 x))) 1.0)
#?-    2.0
#?-    1.5
#?-    1.6666666666666665
#?-    1.6
#?-    1.625
#?-    1.6153846153846154
#?-    1.619047619047619
#?-    1.6176470588235294
#?-    1.6181818181818182
#?-    1.6179775280898876
#?-    1.6180555555555556
#?-    1.6180257510729614
#?-    1.6180371352785146
#?-    1.6180327868852458
1.6180327868852458

練習問題 1.36

　既に表示はしているので log の計算のみ。

gosh> (fixed-point (^[x] (/ (log 1000) (log x))) 1.1)
#?-    72.47657378429035
#?-    1.6127318474109593
#?-    14.45350138636525
#?-    2.5862669415385087
#?-    7.269672273367045
#?-    3.4822383620848467
#?-    5.536500810236703
#?-    4.036406406288111
#?-    4.95053682041456
#?-    4.318707390180805
#?-    4.721778787145103
#?-    4.450341068884912
#?-    4.626821434106115
#?-    4.509360945293209
#?-    4.586349500915509
#?-    4.535372639594589
#?-    4.568901484845316
#?-    4.546751100777536
#?-    4.561341971741742
#?-    4.551712230641226
#?-    4.558059671677587
#?-    4.55387226495538
#?-    4.556633177654167
#?-    4.554812144696459
#?-    4.556012967736543
#?-    4.555220997683307
#?-    4.555743265552239
#?-    4.555398830243649
#?-    4.555625974816275
#?-    4.555476175432173
#?-    4.555574964557791
#?-    4.555509814636753
#?-    4.555552779647764
#?-    4.555524444961165
#?-    4.555543131130589
#?-    4.555530807938518
#?-    4.555538934848503
4.555538934848503

36ステップだった。

　次に平均緩和法を利用した式の変形と変形後の結果は以下のようになる。

f:id:lhcpr:20170705123143p:plain:w300

gosh> (fixed-point (^[x] (/ (+ (/ (log 1000) (log x)) x) 2)) 1.1)
#?-    36.78828689214517
#?-    19.352175531882512
#?-    10.84183367957568
#?-    6.870048352141772
#?-    5.227224961967156
#?-    4.701960195159289
#?-    4.582196773201124
#?-    4.560134229703681
#?-    4.5563204194309606
#?-    4.555669361784037
#?-    4.555558462975639
#?-    4.55553957996306
#?-    4.555536364911781
4.555536364911781