数値計算アルゴリズムの例: 二次方程式

◎二次方程式

数値計算が必要になる簡単な例として、二次方程式をとりあげてみよう。 たとえば、

• 直接測れない高さや距離を求める場合や、家具が戸口を通れるかどうか 調べる場合に、直角三角形の二辺の長さからもうひとつの辺の長さを求め たいということは日常にあると思う。ピタゴラスの定理(三平方の定理)

    A2 + B2 = C2
  

  (A: 底辺の長さ; B: 高さ; C: 斜辺の長さ)を使えばよいが、 このうちひとつ、たとえばBが未知だとすれば、これは一種の二次方程式 になる。
• 点光源から光が出ているとする(音でも基本的には同じ)。途中での減衰 (吸収や反射)が無視できるとすれば、ある距離の場所に届く、時間あたり、 (光源と受ける点を結ぶ直線に垂直な面の)面積あたりのエネルギー量は、 光源からの距離の2乗に反比例、言いかえれば(-2)乗に比例する。
  強さ(時間あたりに出すエネルギー量)の違う2つの光源からのエネルギー量が 等しくなるところは、どういう位置にあるだろうか。まず、2光源を結ぶ線上 での位置を求めてみよう。光源の位置を x₁とx₂、強さを Q₁, Q₂ とし、 受ける位置を x とすれば、

    Q1 ・ (x - x1)-2 = Q2 ・ (x - x2)-2
  

  これは(-2)乗を含んでいるので2次方程式ではないが、近似・変形すれば2次方程式 に持ちこめるのではないだろうか(未確認)。

◎二次方程式をとく

二次方程式については、解の公式がある。代数方程式(「xの多項式 = 0」 の形)のうち4次までのものについて公式が知られている。 それより高い次数のものについては、直接解を求める式はなく、 数値計算の問題としては、有限桁数での近似解を求める方法がいろいろと くふうされている。

二次方程式

   a x2 + b x + c = 0

    - b + sqrt ( b2 - 4 a c )
        -
   ---------------------------
                2 a

a, b, cを読みこんで、2つの解を計算するプログラムを書いてみよう。 まずは、公式をすなおにプログラム化すればよい。 実数の平方根を計算する関数は、Pascal, Fortran, Cのいずれでも sqrt である。(一部のBASICではsqrというが、Pascalのsqrは2乗を計算するので 区別に注意すること。Fortranではbの2乗は b**2 と書く。) 例文は nizi2.p, nizi2.f, nizi2.c にある。 ここでは(わざと)単精度を使っているが、実用には倍精度のほうがよい。

なお、 nizi.p, nizi.f, nizi.c も同じ計算だが、 もしこれをコンパイラがすなおに機械語に翻訳したとすると、 平方根の計算を2回することになるので、 nizi2.{p,f,c}よりも実行効率が悪い。 しかし、多くのコンパイラには、これを1回ですませるように 「最適化」するオプションがある。それを使うとすれば、 nizi.{p,f,c}のほうが単純なだけよいプログラムということになる。

この例文は結果が実数になる場合だけを想定している。つまり 引き数(かっこの中)が正または0でないと、動作しない(エラーメッセージ を出して、止まるかまたは無意味な結果を出すはずである)。したがって、 解が複素数になる場合までカバーしようとすれば、b² - 4 a c の値 によって場合分け(if文による条件分岐)をしなければならない。

Fortranでは複素数型を使うこともできるが、ここではそうしても 特に便利にならないので、実数型(倍精度を含む)でやってみてほしい。

直角三角形の2辺からもう1辺を求める問題は、この公式を使っても よいが、もっと簡単にできるので考えてみてほしい。

◎二次方程式と浮動小数点表現のまるめ誤差

公式をそのままプログラム化したものは、仮に計算機の実数型の 有効数字が無限にあれば、正しい結果を出すはずだが、実際には有限 であるために、誤差が生じる。これをまるめ誤差という。

数値計算法の古典的教科書の著者でもあるフォーサイス (G. E. Forsythe)先生が、計算をするときはこの誤差を考慮して 計算手順(アルゴリズム)を考えなければならないことを、二次方程式を例に 述べている。詳しくは本を見てほしい。

• フォーサイス、「コンピュータで二次方程式を解く」。 アメリカ数理科学研究委員会編、本間龍雄監訳、1972年、「数理科学の世界」、 講談社ブルーバックスB-198、363 -- 388.

この解説では、10進法8桁の浮動小数点表現を持ち、 指数の範囲は絶対値30までであるような計算機を想定している。 まず、さきに述べた公式をすなおにプログラム化してみる。

1. a = 6.0, b = 5.0, c = -4.0
   これはなんの問題もなく、sqrt関数がうまくできていれば、 計算機の精度いっぱいの答えがえられる。
2. a = 1.0, b = -1.0e5, c = 1.0
   これは精度落ちを起こす例である。 これを精度よく計算する方法としては、 (1/x)に関する二次方程式に変形してからプログラム化する方法がある。
3. a = 6.0e30, b = 5.0e30, c = -4.0e30
   これは(1)と同じ答えになるはずだが、途中で 浮動小数点の有効範囲を越えてしまう。ただし、係数を aで割ってから計算を始めれば問題はない。
4. a = 1.0e-30, b = -1.0e30, c = 1.0e30
   解の一方は1.0e60に近い値になるので計算できないが、もうひとつの 1に近い解を求めるアルゴリズムは考えられる。ただしくふうが必要である。
5. a = 1.0, b = 4.0, c = 3.9999999
   これは(2)とは別の意味の精度落ちを起こす例である。 中間結果(sqrtの引き数)の桁数をふやす(倍精度を使う)以外に よい対策はない。

浮動小数点のまるめ誤差の問題はいろいろなところで顔を出すが、 特に、「よく似た数値の引き算」 (あるいは、絶対値が近くて符号が反対の数のたし算)をすると 精度が落ちる。「桁落ち」というのはふつう、このことである。 たとえば、

       1.0000001
    -) 1.0000000
    ------------
       0.0000001

現実には実数の式通りにプログラムする ことが多いが、倍精度を使い、さらに結果のわかった問題で テストしてみて、おかしいと思ったら桁落ちが おきそうなところをチェックすることを勧める。

◎改訂版の例文

必ずしも完全ではないが、多くの場合に実用的な形にした 2次方程式を解くプログラムを nizi3.p, nizi3.f, nizi3.c に用意した。

次の点を考慮したものである。

• 複素数解を「99.999 +/- 99.999 i」のような形で 表示するようにした。 (重解の場合も場合分けするべきかもしれないが、ここでは 他の実数解の場合と同じ扱いにしている。)
• 実数型の変数や定数は倍精度を使った。ただし、このため、 Pascalの場合は、DEC Pascal独自の表現(DOUBLE型、指数のd0)を使っている ので、Pascal言語の標準から離れてしまった。
• 実数解の場合、b +/- sqrt(...) の加減による 桁落ちを防ぐため、絶対値の大きいほうの解を解の公式から求め、 もうひとつの解は「解と係数の関係」を使って求める形にした。 (Fortranの場合、IF...ELSEの代わりにsign関数を使うと プログラムを少し短くできるが、 他の言語の例文と対応する形にしておいた。)

[目次へ]