数値計算と誤差

計算機が扱う数は桁数が有限であるために、数値計算を行う上で必ず誤差が生じます。今回は、数値計算の誤差に関する簡単な演習を行います。


2次方程式の解

まず、2次方程式の解の公式を数値計算することにより生じる数値誤差とその対策法をかんがえてみましょう。2次方程式には解の公式があるのでプログラムを作らずとも解を求めることはできるのですが、数値的に扱いたい場合も多々あります。2次方程式

ax2 + bx + c = 0 (a ≠ 0)

の解は

x = {-b ± (b2 - 4ac)1/2} / (2a)

でしたので、ますはこの公式をそのまま使ってプログラムを書いてみましょう。

     1    program quad
     2      implicit none
     3    
     4      real a,b,c,x1,x2
     5    
     6      write(*,'(a,$)') 'input coefficients a, b, c: '
     7      read(*,*) a,b,c
     8    
     9      x1 = ( -b + sqrt( b*b - 4.0*a*c ) ) / ( 2.0*a )
    10      x2 = ( -b - sqrt( b*b - 4.0*a*c ) ) / ( 2.0*a )
    11    
    12      write(*,'(2e15.6)') x1,x2
    13    
    14    end program quad

○課題○ 入力する係数の値をいろいろ変えて、プログラムを実行してみましょう。


このプログラムでは虚数解が扱えないので、すぐに計算が破たんします。虚数解を扱うプログラムを作る前に、このプログラムで数値誤差に関する実験をおこないましょう。次のような方程式を考えます。

(x − 1.23456 ∗ 103) ∗ (x − 1.23456 ∗ 10-3) = 0

すなわち

(x2 − 1234.56123456 ∗ x + 1.5241383936) = 0

です。

とりあえず、 a = 1.0、b = − 1234.56123456、c = 1.5241383936 として実行してみましょう。結果はどうなりましたか？小さいほうの解が正しく求まっていませんね。理由を考えてみましょう。

まず大きいほうの解は、

x1 = {-b + (b2 - 4ac)1/2} / (2a)

で計算されます。また小さいほうの解は、

x2 = {-b − (b2 - 4ac)1/2} / (2a)

で求められます。ここで右辺の分子に着目しましょう。今回の問題では b の値は負ですから、x1 の右辺の分子は 「(正の数 + 正の数)」になっています。また x2 の方は、「(正の数 − 正の数)」になっています。ところで今回の問題の方程式の係数では、 − b と (− b2 − 4ac)1/2 がほとんど同じ値になっています (b2 に比べて 4ac は無視できるくらい小さくなっています)。

このように2つの数の値がほとんど同じ場合、それらの間で引き算すると、有効桁数が著しく失われてしまいます。例えば 9.87654321 から 9.87654312 を引くと、結果は 0.00000009 になります。はじめ9桁の有効桁を持っていた数が、互いに引き算をした結果、有効桁がわずか1桁になってしまいました。これでは計算の精度が保てません。こういった理由で小さいほうの解 (x2) が正確に計算できなかったわけです。これを桁落ちと言います。

では桁落ちを避けるためにはどうしたらよいでしょうか。引き算を含む方の解(今回の問題では x2。x1とx2のどちらになるかは当然 b の符号によります)を別の方法で求めればよいわけです。ここでは解と係数の関係を使って x1 から

x2 = c / a / x1

として求めることにします。 b の符号が正の場合は先に x2 を求め、x1 を上式と同様に x2 から求めます。以下にプログラム例を示します。 
     1    program quad1
     2      implicit none
     3    
     4      real a,b,c,x1,x2
     5    
     6      write(*,'(a,$)') 'input coefficients a, b, c: '
     7      read(*,*) a,b,c
     8    
     9      if ( b .ge. 0.0 ) then
    10        x2 = ( -b - sqrt( b*b - 4.0*a*c ) ) / ( 2.0*a )
    11        x1 = c / a / x2
    12      else if ( b .lt. 0.0 ) then
    13        x1 = ( -b + sqrt( b*b - 4.0*a*c ) ) / ( 2.0*a )
    14        x2 = c / a / x1
    15      end if
    16    
    17      write(*,'(2e15.6)') x1,x2
    18    
    19    end program quad1

プログラムを入力して実行してみてください。今度はx1、x2 両方の解が正確に求まりました。

○課題○　各自、実験の方程式のような桁落ちの可能性のある例を考えて、修正前と修正後のプログラムで実行し、結果を比較しなさい。

○課題○　判別式 D = b2 - 4ac の値によって処理の制御をすることにより、虚数解も扱えるプログラムを作成しなさい。

虚数解も扱える2次方程式の解のプログラム例をこちらに置いておきます。


多項式の和

次に以下のような多項式

(1) f (x) = a0 + a1 x + a2 x2 + … + an xn

の値を求めるプログラムを作成しながら、数値誤差とサブルーチンの作り方を勉強します。

ところで多項式のプログラムを作る前に、以下の簡単なプログラムを実行してみてください。

     1    program jyoho
     2      implicit none
     3    
     4      real a1,a2,b1,b2
     5      integer i
     6    
     7      a1 = 1.0e7
     8      a2 = 1.0e7
     9      b1 = 1.0
    10      b2 = 0.1
    11    
    12      do i = 1, 10000
    13        a1 = a1 + b1
    14        a2 = a2 + b2
    15      end do
    16    
    17      write(*,'(2f9.0)') a1,a2
    18    
    19    end program jyoho

さて結果はどうなりましたか？予想される結果は

a1 = 10010000.

a2 = 10001000.

ですが、実際は

a1 = 10010000.

a2 = 10000000.

となってしまいました。なぜ 0.1 は加算されなかったのでしょうか。これを正確に理解するにはコンピュータ内部での実数の表現方法を知る必要がありますが、ここでは「実数は有限の有効桁しか持っていない」ということを理解しておけば十分です。単精度の実数の有効桁数は、10進法で約8桁です。そうすると、10000000 に 0.1 を加えた場合、10000000.1の下線の部分は9桁目になるため、無視されてしまいます。このように絶対値の非常に大きな値に、絶対値の非常に小さな値を加えたり引いたりした場合、小さいほうの数が計算結果に全く反映されないことが生じます。これを情報落ちとか積み残しと言います。

○課題○　上記のプログラムで、a1、a2、b1、b2 を倍精度で宣言した場合、結果がどう変わるか確かめましょう。単精度は「real」で宣言しましたが、倍精度は「double precision」で宣言します。


多項式の和を計算する場合、式(1)をそのまま使ったのでは、xの値や係数によっては、この情報落ちが生じる可能性があります。そこで式(1)を以下のように変形します。

(2) f (x) = (… ((an x + an-1) x + … + a1) x + a0

このようにすると、情報落ちをある程度避けることができます。次のプログラムは、

(3) f (x) = 0.10204 + 1.203 x + 0.2335 x2 + 1.786 x3 + 2.334 x4 + 3.4478 x5

の x = π のときの値を求めるものです。

     1    program poly
     2     implicit none
     3    
     4      double precision a(0:5),x,fx
     5    
     6      data a / 0.10204d0, 1.203d0, 0.2335d0, 1.786d0, 2.334d0, 3.4478d0 /
     7    
     8      x = acos( -1.0d0 )
     9      fx = (((( a(5)*x + a(4) )*x + a(3) )*x + a(2) )*x + a(1) )*x + a(0)
    10    
    11      write(*,'(2e15.6)') x,fx
    12    
    13    end program poly

ここではπの値は acos(-1.0d0)で与えています。

○課題○　係数 a や x、fx を単精度で与えて式(3)のように計算すると、例えば x の値が非常に小さい、あるいは大きいとき、情報落ちすることを確かめなさい。


このプログラムは式(3)の多項式にしか使えません。汎用性のあるサブルーチン(関数)をつくってみましょう。Fortran 90プログラミングのウェブページの8章「副プログラム」も参照してください。ここではサブルーチンの名前は poly とでもしましょう。

fortranのメインプログラムで、以下のようにサブルーチンを呼びます。

call poly ( a, n, x, fx)

サブルーチン側では以下のようにします。

subroutine poly ( a, n, x, fx)

  implicit none
　double precision a(0:n), x, fx
　
　　多項式の和の計算

end subroutine

肝心の多項式の和を計算する部分はどうすればよいでしょうか。式(2)の特徴をつかむ必要があります。一番内側の括弧の中から順に計算を進めてみましょう。

　　　第1段　　an x + an-1

　　　第2段　　(an x + an-1) x + an-2

　　　第3段　　((an x + an-1) x + an-2) x + an-3

この式からわかるように、下線部が前段の計算結果になっています。すなわち前段の計算結果に x をかけて、次の係数を加えるという構造になっています。これをプログラム風に書いてみると

　　　第0段　　fx = a(n)

　　　第1段　　fx = fx ∗ x + a(n-1)

　　　第2段　　fx = fx ∗ x + a(n-2)

　　　第3段　　fx = fx ∗ x + a(n-3)

　　　……

　　　第n段　　fx = fx ∗ x + a(0)

といった具合になります。第1段以降は同様な手続きの繰り返しですから、doループの中に入れることができます。すなわち、

     fx = a(n)
     do k = n-1, 0, -1
        fx = fx * x + a(k)
     end do

のようにできます。

○課題○　上記の解説を参考に、多項式の和を求めるサブルーチンを完成させなさい。

上の課題のプログラム例はこちらにありますが、まずは各自でつくってみましょう。