平方根算出のための級数展開と近似値計算

2015/10/25

    

 「数学という学問I 〜 概念を探る」(志賀 浩二 著、筑摩書房 2011年)を読むと、級数展開に関する箇所で、ニュートンが sqrt(1+x^2) を展開させたという記述が出てくる。なお、ここで sqrt() は平方根の意味。

    sqrt(1+x^2) = 1 + 1/2x^2 - 1/8x^4 + 1/16x^6 - 5/128x^8 ……

 なぜ上のような展開式になるのかの解説はない。私は、ここでつまずいた。

 しかも、この後、2項定理や三角関数の級数展開へと話が急展開する。とても追いつかない。そこで、平方根のところで立ち止まることにした。

    

1. 平方根算出のための級数展開

 平方根を算出するための展開式を求める時の鍵は、とりあえずその式を次のように仮定してみることにある。(1+x)の平方根を求めるケースを記す。

    sqrt(1+x) = a + bx + cx^2 + dx^3 + ex^4 + ……

 このように仮定しても大丈夫なのかどうか、私は確信を持てないが、多くの解説でこの仮定が採用されている。

 この式で両辺を2乗すると下のようになる。

    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=1, 2ab=1 というのはともかく、それ以外の (2ac + b^2) とか (2ad + 2bc) がすべて0にならなければならない。

 「……」のところが無限に続くということで、かろうじて帳尻合わせできるような気がするが、なんだか心許ない。

 ともあれ、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

    

 なかなか面倒な計算だ。ただ、救いなのは、1枚ずつ皮をはぐようにして値を計算できることである。

 aが分かればbが分かり、bが決まればcが定まる。こういう計算はパソコンにやってもらうのがいい。ということで、プログラムを組んでみることにした。プログラムを組むプロセスで、その仕組み・構造を知ることにつながるかなという思いもある。

 なお、一般化された2項定理は、まだ十分飲み込めていないので用いない。


2. 平方根の級数展開のためのプログラム

 ここでは下の二つを行うプログラムを取り上げる。計算よりも、先ずは方程式の表示に焦点を当てる。用いるプログラミング言語はruby。

    

    

 もちろん、式を表示するだけではパソコンを使う意味がないので、表示のためのプログラムを少しだけ書き換えて、平方根の近似値計算も行う。

 難しいプログラムにはならないが、ややこしいので順をおってかいてみたい。

    

(1) 2乗してできる展開式

 最も簡単な (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 であり、0ではない。1/4x^2 という、つじつまの合わない項が残ってしまう。

 しかし、この矛盾には目をつぶることにする。無限の級数であれば、この矛盾は解決できるはずだと思い込むわけだ。

    

 ところで、a, b の出現の仕方はどうだろうか。aa, ab, ba, bb の4通りあるが、abとbaは同じ組合せとみなして 2ab と表現する。

 この出現の仕方は、a, b の2つの球が入った袋から1個の球を取り出して戻す作業を2回やる場合の、すべてのパターンを記述するのと同じことである。なので、パターンは 2*2=4通りある。

 球が a, b, c の3個の場合は、3*3=9通りになる。

 3個の球の各々に0〜2の番号を割り当てた場合、その組合せを出力するプログラムは下のようになる。

    

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] という組合せが2回あったことを記録できる。

 3個の球を2回取り出すケースについて、その組合せのパターンと出現回数を記録するプログラムを掲げてみよう。

    

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]  # 未定義なら0に初期化
    hs[comb] += 1  # 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] の各々について、第1要素と第2要素を加算すると、それがxの何乗かの指数を示すことになる。

 [0, 1] であれば abx、[1, 2] なら bcx^3 である。

 以上がプログラムの骨子。あとは、丁寧に処理を繰り返す仕組みを書けばいい。

    

 プログラムそのものは、少し長くなるのでここには掲げないが、sqrt01.rb としてアップロードしておくので興味のある人はダウンロードしてほしい。

 計数を a〜e の5個にした時のプログラムの出力結果を下に掲げる。一つの項を1行に書き出す形になっている。

    

    (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〜4を二つ組み合わせるため、そのパターンは [0, 0] 〜 [4, 4] の25通りある。そして、組の二つの数を足すと指数になる。結局、指数は0〜8のどれかになる。

 指数0〜8の各々が何パターンから構成されるかをみると、最も少ないのが [0, 0] と [4, 4] で、それぞれ一つしかない。

 最も多いのは中央に位置する4で、[0, 4], [4, 0], [1, 3], [3, 1], [2, 2] の5パターンだ。

 このパターン数を指数0〜8の順で並べると、1, 2, 3, 4, 5, 4, 3, 2, 1 となる。階段を上って、頂点までいったら下りてくるかたちである。頂点のところのパターン数は、係数の個数になる。

 頂点まで上るプロセスで、a〜eの係数が1回は出現するため、係数の値を求める時は上り階段の方にだけ注目すればよい。

 下り階段の部分は、無限級数でなく有限になった途端に「帳尻の合わない部分」、つまり誤差になってしまう。で、係数の個数が多いほど、この誤差は小さくなる。

 無限級数は、ひたすら階段を上り続ける形になるので、実質的に誤差が0になる。

    

 なお、sqrt01.rb で係数の個数を27以上に設定して実行すると、アルファベットの個数を超えるためa〜zでない文字が割り当てられる。そのため意味不明な式になってしまう。

 係数をa〜zでなく、a0, a1, a2 などにすれば、そうした制限はなくなる。ただ、数時がやたらと出てくる式になってしまう。それでもよければ sqrt01b.rb の方を実行してみてほしい。


(2) 平方根を求める級数展開の式の表示

 今度は、sqrt(1+x) を求める級数展開の式を表示するプログラムを考える。無限個を扱うのは無理なので、係数の個数は、あらかじめ変数nに代入しておくことにする。

 ここでは sqrt(1+x) を取り上げるので、2乗した時の展開式の x^0 と x^1 にかぶさる数がそれぞれ1ということになる。そして、x^2, x^3 , x^4 …… にかぶさる数がすべて0になるものと仮定する。

 x^0にかぶさる数a^2は1であり、x^1にかぶさる数2abも1だ。それ以外は0になるよう c, d, e, …… を定めればいい。

 そして、都合のいいことに、これら係数の値は、一つずつ算出できる。一度に二つの係数を算出する必要はない。なので、連立方程式を解くアルゴリズムを盛り込まなくても済む。

 ということで、前述の2乗の展開式を表示する sqrt01.rb を素材にして、簡単な計算処理を組み込めばプログラムができる。

    

 プログラムを書く時のポイントを、係数が a, b, c の三つのケースでみてみよう。

sqrt(1+x) = a + bx + cx^2
1+x = (a + bx + cx^2)^2

 ここで、左辺の x^0 と x^1 にかぶさる数(それぞれ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) と書く。1番目が分子、2番目が分母を示す。2番目の分母を省略すると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 にかぶさる数は1なので(これは original_vals に既に記録済み)、2ab = 1 となるが、aの値は1であるから(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 には記録されていないが、記録されていないものは、今回、0とみなす。

    

 こうした処理を具体化したのが sqrt02.rb というプログラムである。長くなるのでここには掲げないが、必要ならダウンロードしてほしい。

 参考まで、係数が10個の場合の出力結果を下に掲げてみよう。プログラムは、一つの項を1行にして出力するが、ここでは2行の形にして掲げる。

    

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 以降は、帳尻の合わない誤差として切り捨てる。

 誤差がどのような様相を呈するかも興味深いが、とりあえず忘れることにする。


(3) 平方根を求める級数展開の式による近似値計算

 級数展開の式が得られれば、実際に計算するのは容易だ。プログラムを書く場合も、sqrt02.rb を僅かに書き換えてやれば目的のプログラムになる。

 ただ、xにどんな値を代入するかが問題になる。展開式を見て想像できるとおり、xに1以上の値を代入しても、あまりうまくいきそうにない。

 たとえば、2の平方根を求める場合、sqrt(1+x) の式なのだから x=1 とすればいいようなものだが、1は何乗しても1なので、級数展開の式がなかなか終息しない。

 3の平方根を算出するために x=2 とすれば、結果としてとんでもない数になってしまう。終息どころか、大波の嵐だ。できれば、xには -1 < x < 1 の範囲の値を入れたい。

    

 この問題を解決する最も簡単な方法は、逆数を用いるというものだろう。

 sqrt(2) を算出するのに sqrt(1/2) を計算する。そして、得られた結果を逆数にすればいい。1+x=1/2 になるようにするには、x=-1/2 とする。その値を式に代入して計算し、得られた結果を逆数にする。

    

 この方法で計算するのが sqrt03.rb というプログラムである。

 プログラムの第3行目に square = 5 というのがある。これは sqrt(5) を算出するという指定になる。これを2にすれば sqrt(2) が算出される。

 その次の行に n = 50 とあるが、これは係数の個数を指定するものである。この数を増やすほど計算の精度が上がる。ただ、実際のパソコンでの処理で、どのくらいまで有効かは確認していない。

 プログラムの最後では、rubyに予め用意されている組み込み関数 sqrt() を使った計算結果も出力するようにした。この本当の値とどれくらい違うかを見ることができる。

    

 今回の話題については、この辺でおわりにしたい。平方根等の近似値計算ではニュートン法という有名な方法があるが、別の機会に触れられればと思う。

 平方根を求めるための級数展開ができれば、三角関数の級数展開も行えるようになるみたいだが、まだ大雑把にしか咀嚼できていない。次回、書きたいと思うのだが、どうなることやら。

    

〜 「平方根算出のための級数展開と近似値計算」おわり 〜

Copyright (C) T. Yoshiizumi, 2015 All rights reserved.


「数学散歩の小道」のページへ

トップページへ