「数学という学問I 〜 概念を探る」(志賀 浩二 著、筑摩書房 2011年)を読むと、級数展開に関する箇所で、ニュートンが sqrt(1+x^2) を展開させたという記述が出てくる。なお、ここで sqrt() は平方根の意味。
sqrt(1+x^2) = 1 + 1/2x^2 - 1/8x^4 + 1/16x^6 - 5/128x^8 ……
なぜ上のような展開式になるのかの解説はない。私は、ここでつまずいた。
しかも、この後、２項定理や三角関数の級数展開へと話が急展開する。とても追いつかない。そこで、平方根のところで立ち止まることにした。
平方根を算出するための展開式を求める時の鍵は、とりあえずその式を次のように仮定してみることにある。(1+x)の平方根を求めるケースを記す。
sqrt(1+x) = a + bx + cx^2 + dx^3 + ex^4 + ……
このように仮定しても大丈夫なのかどうか、私は確信を持てないが、多くの解説でこの仮定が採用されている。
この式で両辺を２乗すると下のようになる。
1 + x = (a + bx + cx^2 + dx^3 + ex^4 + ……)^2
ここで右辺を展開すると次のような式になる。
a^2 + 2abx + (2ac + b^2)x^2 + (2ad + 2bc)x^3 + ……
これが 1+x と一致するように a, b, c, d, …… の係数を定めれば、平方根を求める級数展開の式が得られると考える。
といっても、a^2=１, 2ab=１ というのはともかく、それ以外の (2ac + b^2) とか (2ad + 2bc) がすべて０にならなければならない。
「……」のところが無限に続くということで、かろうじて帳尻合わせできるような気がするが、なんだか心許ない。
ともあれ、a, b, c, d, e くらいまでを具体的にみてみよう。
a^2 = 1 → a = 1
2ab = 1 → 2b = 1 → b = 1/2
2ac + b^2 = 0 → 2c + 1/4 = 0 → c = -1/8
2ad + 2bc = 0 → 2d + -1/8 = 0 → d = 1/16
2ae + 2bd + c^2 = 0 → 2e + 1/16 + 1/64 = 0 → e = -5/128
なかなか面倒な計算だ。ただ、救いなのは、１枚ずつ皮をはぐようにして値を計算できることである。
aが分かればbが分かり、bが決まればcが定まる。こういう計算はパソコンにやってもらうのがいい。ということで、プログラムを組んでみることにした。プログラムを組むプロセスで、その仕組み・構造を知ることにつながるかなという思いもある。
なお、一般化された２項定理は、まだ十分飲み込めていないので用いない。
ここでは下の二つを行うプログラムを取り上げる。計算よりも、先ずは方程式の表示に焦点を当てる。用いるプログラミング言語はruby。
もちろん、式を表示するだけではパソコンを使う意味がないので、表示のためのプログラムを少しだけ書き換えて、平方根の近似値計算も行う。
難しいプログラムにはならないが、ややこしいので順をおってかいてみたい。
最も簡単な (a + bx)^2 を考えてみよう。これは、a^2 + 2ab + b^2x^2 である。
すぐに矛盾が露呈してしまうが、この式が 1+x と等価だとすると、a^2=1, 2ab=1 であると見ることができる。すると、a=1, b=1/2 が導き出される。
もちろん、これだと最後の項の b^2x^2 が消えない。b^2=1/4 であり、０ではない。1/4x^2 という、つじつまの合わない項が残ってしまう。
しかし、この矛盾には目をつぶることにする。無限の級数であれば、この矛盾は解決できるはずだと思い込むわけだ。
ところで、a, b の出現の仕方はどうだろうか。aa, ab, ba, bb の４通りあるが、abとbaは同じ組合せとみなして 2ab と表現する。
この出現の仕方は、a, b の２つの球が入った袋から１個の球を取り出して戻す作業を２回やる場合の、すべてのパターンを記述するのと同じことである。なので、パターンは 2*2=4通りある。
球が a, b, c の３個の場合は、3*3=9通りになる。
３個の球の各々に０〜２の番号を割り当てた場合、その組合せを出力するプログラムは下のようになる。
for i in 0..2
for j in 0..2
p [i,j]
end
end
上が今回のプログラムの最も基本となる部分である。
ただし、これだと [0, 1] と [1, 0] が別々のものとして扱われてしまう。そこで、ソート(整列)の処理を加えて、どちらも [0, 1] にした上で、その数をカウントする。
カウントする場合、rubyのHashを利用する。Hashは、perlでは連想配列、pythonでは辞書といわれているようだが、通常の配列が要素に整数しか使えないのに対し、Hashは文字列等を要素にすることができる。
たとえば、リンゴの値段を記録したい時に、hs[“リンゴ”] = 200 のように使える。
そして、rubyではHashの要素に文字列でなく配列を指定することもできる。
hs[[0,1]] = 2 とすれば、[0,1] という組合せが２回あったことを記録できる。
３個の球を２回取り出すケースについて、その組合せのパターンと出現回数を記録するプログラムを掲げてみよう。
hs = Hash.new
n = 3
for i in 0...n # i: 0, 1, 2
for j in 0...n # j: 0, 1, 2
comb = [i, j].sort
hs[comb] = 0 unless hs[comb] # 未定義なら０に初期化
hs[comb] += 1 # １を加えてカウント
end
end
hs.each do |comb, count| # Hashに記録した内容を表示
printf("%s = %d\n", comb.inspect, count)
end
上のプログラムを実行すると、すべて数字で出力されるので分かりにくいが、[0, 0] が a^2、[0, 1] がab、[1, 2] ならbcを意味する。
また、[0, 0] とか [1, 2] の各々について、第１要素と第２要素を加算すると、それがxの何乗かの指数を示すことになる。
[0, 1] であれば abx、[1, 2] なら bcx^3 である。
以上がプログラムの骨子。あとは、丁寧に処理を繰り返す仕組みを書けばいい。
プログラムそのものは、少し長くなるのでここには掲げないが、sqrt01.rb としてアップロードしておくので興味のある人はダウンロードしてほしい。
計数を a〜e の５個にした時のプログラムの出力結果を下に掲げる。一つの項を１行に書き出す形になっている。
(a + bx + ... + ex^4)^2 =
a^2 +
2abx +
(2ac + b^2)x^2 +
(2ad + 2bc)x^3 +
(2ae + 2bd + c^2)x^4 +
(2be + 2cd)x^5 +
(2ce + d^2)x^6 +
2dex^7 +
e^2x^8
xが何乗かを示す指数に着目してみよう。
a〜eの五つの係数に割り当てられる番号０〜４を二つ組み合わせるため、そのパターンは [0, 0] 〜 [4, 4] の25通りある。そして、組の二つの数を足すと指数になる。結局、指数は０〜８のどれかになる。
指数０〜８の各々が何パターンから構成されるかをみると、最も少ないのが [0, 0] と [4, 4] で、それぞれ一つしかない。
最も多いのは中央に位置する４で、[0, 4], [4, 0], [1, 3], [3, 1], [2, 2] の５パターンだ。
このパターン数を指数０〜８の順で並べると、1, 2, 3, 4, 5, 4, 3, 2, 1 となる。階段を上って、頂点までいったら下りてくるかたちである。頂点のところのパターン数は、係数の個数になる。
頂点まで上るプロセスで、a〜eの係数が１回は出現するため、係数の値を求める時は上り階段の方にだけ注目すればよい。
下り階段の部分は、無限級数でなく有限になった途端に「帳尻の合わない部分」、つまり誤差になってしまう。で、係数の個数が多いほど、この誤差は小さくなる。
無限級数は、ひたすら階段を上り続ける形になるので、実質的に誤差が０になる。
なお、sqrt01.rb で係数の個数を27以上に設定して実行すると、アルファベットの個数を超えるためa〜zでない文字が割り当てられる。そのため意味不明な式になってしまう。
係数をa〜zでなく、a0, a1, a2 などにすれば、そうした制限はなくなる。ただ、数時がやたらと出てくる式になってしまう。それでもよければ sqrt01b.rb の方を実行してみてほしい。
今度は、sqrt(1+x) を求める級数展開の式を表示するプログラムを考える。無限個を扱うのは無理なので、係数の個数は、あらかじめ変数nに代入しておくことにする。
ここでは sqrt(1+x) を取り上げるので、２乗した時の展開式の x^0 と x^1 にかぶさる数がそれぞれ１ということになる。そして、x^2, x^3 , x^4 …… にかぶさる数がすべて０になるものと仮定する。
x^0にかぶさる数a^2は1であり、x^1にかぶさる数2abも1だ。それ以外は０になるよう c, d, e, …… を定めればいい。
そして、都合のいいことに、これら係数の値は、一つずつ算出できる。一度に二つの係数を算出する必要はない。なので、連立方程式を解くアルゴリズムを盛り込まなくても済む。
ということで、前述の２乗の展開式を表示する sqrt01.rb を素材にして、簡単な計算処理を組み込めばプログラムができる。
プログラムを書く時のポイントを、係数が a, b, c の三つのケースでみてみよう。
sqrt(1+x) = a + bx + cx^2
1+x = (a + bx + cx^2)^2
ここで、左辺の x^0 と x^1 にかぶさる数(それぞれ１)を配列 original_vals にセットする。具体的には下のとおり。
original_vals = [Rational(1), Rational(1)]
また、右辺を展開すると x^0 にかぶさる数が a^2 であり、a=1 なので、この値もセットしてしまう。係数 a, b, c の値を配列 coeff_vals に記録するものとすると、下のように書ける。
coeff_vals = [Rational(1)]
プログラムを走らせる過程で b, c の値が判明したら、その都度、この coeff_vals に値を記録していけばいい。
Rational というのは ruby に予め組み込まれているクラスで、有理数を簡単に扱える。内部で分母と分子が保持されており、表示の時に分数形式にできる。もちろん、計算処理の対象として用いることもできる。
Rational(1) は、一種の省略記法で、ほんとは Rational(1, 1) と書く。１番目が分子、２番目が分母を示す。２番目の分母を省略すると１になる。
このクラスには、分数形式の文字列(string)に変換するための to_s、浮動小数点数(float)に変換するための to_f が用意されている。たとえば次のように書く。
rr = Rational(1, 2)
p rr.to_s # -> "1/2"
p rr.to_f # -> 0.5
係数 b, c の算出に話を進めよう。
まだ値が判明していない係数bが出てくるのは、2abx だ。x^1 にかぶさる数は１なので(これは original_vals に既に記録済み)、2ab = 1 となるが、aの値は１であるから(coeff_vals に記録済み) 2b = 1 という式になる。これを解くと b = 1/2 だ。
一つの係数を算出する時の式を一般化して書くと、下のように書ける。未知の係数をxと表記する。
αx + β = γ
この式で、α, β, γ が分かれば、xが機械的に算出される。x = (γ - β) / α である。
係数bを求める場合は、2b + 0 = 1 となるので、b = (1 - 0) / 2 = 1/2 という結果を得る。
係数cについては、(2ac + b^2)x^2 という項だ出てくるので、2c + (1/2)^2 = 0 となり、c = -1/8 を導き出せる。
x^2 にかぶさる数が original_vals には記録されていないが、記録されていないものは、今回、０とみなす。
こうした処理を具体化したのが sqrt02.rb というプログラムである。長くなるのでここには掲げないが、必要ならダウンロードしてほしい。
参考まで、係数が10個の場合の出力結果を下に掲げてみよう。プログラムは、一つの項を１行にして出力するが、ここでは２行の形にして掲げる。
sqrt(1+x) =
1 + 1/2x + -1/8x^2 + 1/16x^3 + -5/128x^4 + 7/256x^5 +
-21/1024x^6 + 33/2048x^7 + -429/32768x^8 + 715/65536x^9
蛇足ながら、sqrt01.rb で係数の個数を10にして実行すると、x^18 まで表示される。
それに対し、sqrt02.rb の方では x^9 で終了している。これは、「上り階段」だけに注目するからである。x^10 以降は、帳尻の合わない誤差として切り捨てる。
誤差がどのような様相を呈するかも興味深いが、とりあえず忘れることにする。
級数展開の式が得られれば、実際に計算するのは容易だ。プログラムを書く場合も、sqrt02.rb を僅かに書き換えてやれば目的のプログラムになる。
ただ、xにどんな値を代入するかが問題になる。展開式を見て想像できるとおり、xに１以上の値を代入しても、あまりうまくいきそうにない。
たとえば、２の平方根を求める場合、sqrt(1+x) の式なのだから x=1 とすればいいようなものだが、１は何乗しても１なので、級数展開の式がなかなか終息しない。
３の平方根を算出するために x=2 とすれば、結果としてとんでもない数になってしまう。終息どころか、大波の嵐だ。できれば、xには -1 < x < 1 の範囲の値を入れたい。
この問題を解決する最も簡単な方法は、逆数を用いるというものだろう。
sqrt(2) を算出するのに sqrt(1/2) を計算する。そして、得られた結果を逆数にすればいい。1+x=1/2 になるようにするには、x=-1/2 とする。その値を式に代入して計算し、得られた結果を逆数にする。
この方法で計算するのが sqrt03.rb というプログラムである。
プログラムの第３行目に square = 5 というのがある。これは sqrt(5) を算出するという指定になる。これを２にすれば sqrt(2) が算出される。
その次の行に n = 50 とあるが、これは係数の個数を指定するものである。この数を増やすほど計算の精度が上がる。ただ、実際のパソコンでの処理で、どのくらいまで有効かは確認していない。
プログラムの最後では、rubyに予め用意されている組み込み関数 sqrt() を使った計算結果も出力するようにした。この本当の値とどれくらい違うかを見ることができる。
今回の話題については、この辺でおわりにしたい。平方根等の近似値計算ではニュートン法という有名な方法があるが、別の機会に触れられればと思う。
平方根を求めるための級数展開ができれば、三角関数の級数展開も行えるようになるみたいだが、まだ大雑把にしか咀嚼できていない。次回、書きたいと思うのだが、どうなることやら。
〜 「平方根算出のための級数展開と近似値計算」おわり 〜
Copyright (C) T. Yoshiizumi, 2015 All rights reserved.