読者です 読者をやめる 読者になる 読者になる

数値計算,洲之内著,石渡改訂,サイエンス社 を読む (1)

このエントリーは、完全に自分用の勉強ページです。

このエントリーでは、
数値計算,洲之内著,石渡改訂,サイエンス社 (2005年4月10日新訂第4刷)
を読み、例題や演習問題を解く。

解くとは言っても、本には解答が付いている。今回、私は実際にパソコンを使って解答を確認していくことにした。原則としてschemeを使い、グラフ描画にはgnuplotを使う。

例題1.1

 \displaystyle
f(x) = \frac{x}{\sqrt{x^2 + 1} + x} \quad \quad (1.2)

を計算するとき、分母を有理化して,

 \displaystyle
f(x)= x \sqrt{x^2 + 1} - x^2 \quad (1.3)

とした方が除算がないので良さそうに見える. しかし,コンピュータで計算するアルゴリズムとしては(1.2)と(1.3)のどちらがよいか?

解答:xが大きいときに(1.2)の方が良い。(1.3)では桁落ちが起きる。

まず、以下の極限値が得られる。
 \displaystyle
\lim_{x\to \infty} f(x)
= \lim_{x\to \infty}
   \frac{1}{\sqrt{1+ \frac{1}{x^2}} +1}
=\frac{1}{2}


ここで実際に、xとして大きな値を選び、コンピュータで計算してみる。

;;式(1.2)
(define (f1 x)
  (/ x (+ (sqrt (+ (* x x) 1)) x)))

;;式(1.3)
(define (f2 x)
  (- (* x (sqrt (+ (* x x) 1))) (* x x)))


;;xy-datは、2次元グラフ用のデータ x, f(x) を
;;標準出力するための関数。
;;interval : x軸の区間
;;step-size : 刻み 
(define (xy-dat f interval step-size)
  (define (xy-dat-p f interval step-size)
    (if (> (car interval) (cdr interval))
      (newline)
      (begin
        (display (car interval))
        (display ",  ")
        (display (f (car interval)))
        (newline)
        (xy-dat-p f (cons (+ (car interval) step-size) (cdr interval)) step-size))))
  (begin
    (newline)
    (display f)
    (newline)
    (xy-dat-p f interval step-size)))

横軸の区間を[10000,100000000]として、グラフをプロットする。
Schemeで次のように実行し、

(xy-dat f1 (cons 10000 100000000) 10000)
(xy-dat f2 (cons 10000 100000000) 10000)

出力されたデータを、それぞれgnuplotで描画した。


式(1.2)
f:id:nibosiiwasi:20170311230104p:plain


式(1.3)
f:id:nibosiiwasi:20170311230109p:plain

考察:式(1.2)はf(x)の値が0.5に収束している。
式(1.3)では、誤差が発生している。
式(1.3)では、引き算で桁落ちが起き、x>=67110000でf(x)=0となった。

疑問点:式(1.3)のグラフで途中、値が振動している部分がある。なぜ?
(やばい、しょっぱなから分からない現象が起こってる)

なぜ振動するのか?
f(x)を別の式で表現してみる。
\displaystyle
f(x) = x\left( \sqrt{x^2 + 1} - x \right)
これを関数f2_ としてschemeでコードを書き、
gnuplotでグラフを描画した。

(define (f2_ x)
  (* x (- (sqrt (+ (* x x) 1)) x)))

(xy-dat f2_ (cons 10000 100000000) 2000)

f:id:nibosiiwasi:20170319230359p:plain

私は自分の目を疑った。このグラフは3つのグラフの中で1番ひどい。
不連続点を含む振動はもちろん真値ではあり得ない。コンピュータ上の計算の誤差によるものだ。しかし、なぜこんな振動が起こる??

演習問題1.3


\hat{x}, \hat{y} をそれぞれ真の値x,y の近似値とし,

誤差および相対誤差をe(\hat{x}), e_R(\hat{x})と表すとき,
次のことを示せ.


{\rm (i)}\quad
e(\hat{x}\pm \hat{y}) = e(\hat{x}) \pm e(\hat{y})


{\rm (ii)}\:\:\:
e(\hat{x} \hat{y})
 \fallingdotseq
\hat{y} e(\hat{x}) + \hat{x} e(\hat{y})

\quad\quad
e_R(\hat{x} \hat{y})
 \fallingdotseq
e_R(\hat{x}) + e_R(\hat{y})

\displaystyle
{\rm (iii)}\:\:
e\left(\frac{\hat{x}}{\hat{y}}\right),
e_R\left(\frac{\hat{x}}{\hat{y}}\right)
 についてはどうなるか.


 (注){\rm (ii)}の近似計算( \fallingdotseq ) には、「e(\hat{x}), e(\hat{y}) , e_R(\hat{x}), e_R(\hat{y})
の絶対値は十分小さく、それらの2次以上の項は無視してよいとする。」
という仮定が必要。問題文が舌足らずである。


解答:
まず、用語の定義を確認する。P4より
誤差:=近似値-真の値
相対誤差:=(近似値-真の値)/真の値


\begin{eqnarray}
{\rm (i)}\:
e(\hat{x}\pm \hat{y})
&=&
(\hat{x}\pm \hat{y})
-
(x\pm y)\\
&=&
(\hat{x}-x) \pm (\hat{y} -y) \\
&=&
e(\hat{x}) \pm e(\hat{y})
\end{eqnarray}


\begin{eqnarray}
{\rm (ii)}\:
e(\hat{x} \hat{y})
&=&
\hat{x}\hat{y} -xy \\
&=&
\hat{y}(\hat{x}-x) + x\hat{y} -xy \\
&=&
\hat{y}e(\hat{x}) + x e(\hat{y}) \\
&=&
\hat{y}e(\hat{x}) + \hat{x} e(\hat{y}) -(\hat{x}-x) e(\hat{y}) \\
&=&
\hat{y}e(\hat{x}) + \hat{x} e(\hat{y}) -e(\hat{x}) e(\hat{y}) \\
&\fallingdotseq &
\hat{y}e(\hat{x}) + \hat{x} e(\hat{y})
\end{eqnarray}


 (注)誤差(の絶対値)が大きいとき、この\fallingdotseq は成り立たない。

\displaystyle
\begin{eqnarray}
e_R(\hat{x} \hat{y})
&=&\frac{\hat{x}\hat{y} -xy}{xy} \\
&=&\frac{\hat{y}(\hat{x}-x)+x\hat{y} -xy}{xy} \\
&=&
\frac{\hat{x}-x}{x} \frac{\hat{y}}{y} + \frac{\hat{y}-y}{y} \\
&=&
e_R(\hat{x})\frac{\hat{y}}{y}+ e_R({\hat{y}}) \\
&=&
e_R(\hat{x})+ e_R({\hat{y}}) + e_R(\hat{x}) \frac{\hat{y}-y}{y} \\
&=&
e_R(\hat{x}) + e_R(\hat{y}) + e_R(\hat{x})e_R(\hat{y}) \\
&\fallingdotseq&
e_R(\hat{x}) + e_R(\hat{y})
\end{eqnarray}


 (注)相対誤差(の絶対値)が大きいとき、この\fallingdotseq は成り立たない。


(iii) 考え中。「どうなるか?」というザックリした問いに対する適切な答えは何だろうか?
誤差の定義と、(ii)まで解いた感じから、微分とほとんど同じことが言えるはずだと予感している。

演習問題1.4

x,y の四捨五入した近似値を3.32, 5.39 とするとき,
(i) x+y (ii) x+0.1y (iii) x+0.01y
の値を, 小数点以下2桁まで求め, その誤差限界を求めよ.
なお, このx,y の値は

\sqrt{11}, \: \sqrt{19}
の近似値である. 今, 有効桁数を倍の長さにとると、

x\fallingdotseq 3.31662, \quad
y\fallingdotseq 5.38516
となる. このことを用いて, 上に求めた近似値と比較してみると, どのようなことがわかるか.

解答:

x,yを四捨五入した近似値を
\hat{x} = 3.32, \: \hat{y} = 5.39 \:
と書く。

このとき、

\begin{eqnarray}
 -0.005 &\le x-\hat{x} \le 0.005 , \\
 -0.005 &\le y-\hat{y} \le 0.005
\end{eqnarray}
である。

(i)

\begin{eqnarray}
 -0.01 &\le& (x+y)-(\hat{x}+\hat{y}) \le 0.01 \\
 (\hat{x}+\hat{y})-0.01 &\le& x+y \le (\hat{x}+\hat{y}) + 0.01 \\
 8.71 - 0.01 &\le& x+y \le 8.71 + 0.01
\end{eqnarray}
したがって、
x+y の小数点以下2桁までの近似値は8.71
誤差限界は0.01である。

(ii)

\begin{eqnarray}
 -0.005 &\le& x-\hat{x} \le 0.005 , \\
 -0.0005 &\le& 0.1(y-\hat{y}) \le 0.0005
\end{eqnarray}
であるから

\begin{eqnarray}
 (\hat{x}+0.1 \hat{y}) -0.0055 &\le& x+0.1y \le (\hat{x}+0.1 \hat{y}) + 0.0055 , \\
 3.859 -0.0055 &\le& x+0.1y \le 3.859 + 0.0055
\end{eqnarray}
したがって
x+0.1y の小数点以下2桁までの近似値は3.86
誤差限界は0.0055+0.001=0.0065である。

(iii)

\begin{eqnarray}
 -0.005 &\le& x-\hat{x} \le 0.005 , \\
 -0.00005 &\le& 0.01(y-\hat{y}) \le 0.00005
\end{eqnarray}
であるから

 3.3739-0.00505 \le x+0.01y \le 3.3739 + 0.00505
したがって
x+0.01y の小数点以下2桁までの近似値は3.37
誤差限界は0.00505+0.0039=0.00895である。

(後半は省略)

演習問題1.6

10進数(decimal number)で表された数を2進数(binary number) で表すアルゴリズムを考えよ. たとえば, 11.3125 は2進数で表すとどうなるだろうか.

解答: 11.3125の例を使って考える。
まず、整数部11と小数部0.3125に分けて考える。最後に合わせればよい。

10進数の非負の整数Dは次のように書くことができる。

D = f_n 2^n + f_{n-1} 2^{n-1} + f_{n-2} 2^{n-2} + \cdots + f_1 2 + f_0,\:

ただし \forall i\in\{0,1,\dots,n\},\; f_i \in \{0,1\}.


このとき、10進数の非負の整数Dを表す2進数は、
f_n f_{n-1}f_{n-2}\dots f_1 f_0 である。


どうやって、これらのf_i を計算すればよいか?
例に戻って考えてみる。
以下のように11を2で割り、商をさらに2で割るという計算を繰り返す。

\begin{eqnarray}
11÷2&=& \fbox{5}余り1 \\
\fbox{5}÷2&=& \fbox{2}余り1 \\
\fbox{2}÷2&=& \fbox{1}余り0 \\
\fbox{1}÷2&=& 0余り1
\end{eqnarray}
最後の行で、商は0、余りは1になる。


このとき、余りは上から順に、f_0, f_1, f_2, f_3
となる。

\begin{eqnarray}
11÷2&=& 5余り \fbox{1} \leftarrow f_0 \\
5÷2&=& 2余り \fbox{1} \leftarrow f_1 \\
2÷2&=& 1余り \fbox{0} \leftarrow f_2 \\
1÷2&=& 0余り \fbox{1} \leftarrow f_3
\end{eqnarray}

確認すると、たしかに
11= 1\times 2^3 + 0\times 2^2 + 1\times 2^1 +1 となっている。

アルゴリズムとしては(上の方法を11に限定せず)、一般の非負の10進数の整数に対して行えばよい。


 f_n 2^n + f_{n-1} 2^{n-1} + f_{n-2} 2^{n-2} + \cdots + f_1 2 + f_0
を2で割る。


このこき、商は
 f_n 2^{n-1} + f_{n-1} 2^{n-2} + f_{n-2} 2^{n-3} + \cdots + f_1,\quad 余りは f_0.

更にその商を2で割ると、商は

 f_n 2^{n-2} + f_{n-1} 2^{n-3} + f_{n-2} 2^{n-4} + \cdots + f_2,\quad 余りは f_1.

...
これを繰り返すと最終的に、

 商が0, \quad 余りがf_n
となる。
したがって、上の余りを並べていけば、求めたかった2進数表示

f_n f_{n-1} \dots f_1 f_0
が得られる。

これをSchemeで書いてみる。

;; di : 10進数の非負の整数とする
;; di->binary : 10進数の非負の整数を2進数に変換する
(define (di->binary di)
  (define (di->binary-p di tmp)
    (if (= di 0)
      (string->number tmp)
      (di->binary-p (quotient di 2)
                 (string-append (number->string (remainder di 2))
                                tmp))))
  
  (if (and (>= di 0) (integer? di))
    (if (= di 0)
      0 ;;string->numberで "" -> #fとなる。仕方なく0は例外として扱う
      (di->binary-p di ""))
    (error "di should be a non-negative integer!!")))

動作確認

(map di->binary '(0 1 2 3 4 5 6 7 8))

;;実行結果
(0 1 10 11 100 101 110 111 1000)


10進数の非負の小数部d (0<=d<1) は次のように書くことができる。

d = f_{-1} \frac{1}{2} + f_{-2} \frac{1}{2^2} + \cdots + f_{-n}\frac{1}{2^n},\:

ただし \forall i\in\{-1,-2,\dots,-n\},\; f_i \in \{0,1\}.


このとき、10進数の非負の小数部dを表す2進数は、
f_{-1} f_{-2} \dots f_{-n} である。

\displaystyle
f_{-1} \frac{1}{2} + f_{-2} \frac{1}{2^3} + \cdots + f_{-n}\frac{1}{2^n}
に2を掛けると
\displaystyle
f_{-1} + f_{-2} \frac{1}{2} + \cdots + f_{-n}\frac{1}{2^{n-1}} となる。
\displaystyle
この整数部は f_{-1}, 小数部はf_{-2} \frac{1}{2} + \cdots + f_{-n}\frac{1}{2^{n-1}}.

更にその小数部に2を掛けると、
\displaystyle
f_{-2}+ f_{-3} \frac{1}{2} + \cdots + f_{-n}\frac{1}{2^{n-2}} となる。

この整数部は f_{-2}, 小数部はf_{-3} \frac{1}{2} + \cdots + f_{-n}\frac{1}{2^{n-2}}.

...
これを繰り返すと、
2つのパターンに分かれる。

  1. 有限回で小数部が0になり、終了する。
  2. 循環小数となり同じ計算を繰り返す。

そして、いずれの場合も
10進数の非負の小数部dを表す2進数は、

f_{-1} f_{-2} \dots f_{-n} であり、
上の整数部を並べることで得られる。


これをSchemeで書いてみる。

;;d_: 10進数の小数。0<= d_ <1 とする。
;;d_->string-binary : 0以上1未満の、10進数の小数を2進数に変換する。出力形式は文字列
(define (d_->string-binary d_)
  (define (d_->string-binary-p d_ tmp count)
    (if (or (= d_ 0) (= count 50)) ;; 小数点以下50桁まで表示
      tmp
      ;;      (string->number tmp) ;;<- これはダメ。理由については(注)参照
      (d_->string-binary-p (- (* d_ 2)
                              (floor (* d_ 2)))
                           (string-append tmp
                                          (number->string (inexact->exact (floor (* d_ 2)))))
                           (+ count 1))))
  
  (cond ((= d_ 0) "")
        ((and (< 0 d_) (< d_ 1)) (d_->string-binary-p d_ "." 0))
        (error "d_ should be a number in [0,1)!!")))

(注)整数部の変換とほぼ同じだが、いくらかややこしい点を含んでいる。

  1. 有限桁で終わらない場合(つまり循環小数の場合)があるので、対策が必要。安直だが、カウンターを用意して途中で打ち切ることにした。上のコードでは小数点以下50桁で打ち切る。
  2. 小数部の2進数は文字列で出力することにした。

Schemeの非正確数は倍精度であり、本来10進数を扱うためのものである。上の小数部の変換コードで倍精度の非正確数を使うと16桁しか表示できない。仕方なく、文字列を使った。


もちろん0と1を扱うために正確数や文字列を使うのは、どう考えてもメモリーを浪費しているはずで、もっと上手いプログラムがあると思う。

ともかく、これで整数部,小数部に関する10進数→2進数の変換コードが書けた。
これらを合わせてコードを完成させる。

(define (decimal->string-binary d)
  ;;内部定義

  ;; di : 10進数の非負の整数とする
  ;; di->binary : 10進数の非負の整数を2進数に変換する
  (define (di->binary di)
    (define (di->binary-p di tmp)
      (if (= di 0)
        (string->number tmp)
        (di->binary-p (quotient di 2)
                      (string-append (number->string (remainder di 2))
                                     tmp))))
    
    (if (and (>= di 0) (integer? di))
      (if (= di 0)
        0 ;;string->numberで "" -> #fとなる。仕方なく0は例外として扱う
        (di->binary-p di ""))
      (error "di should be a non-negative integer!!")))
  
  ;;d_: 10進数の小数。0<= d_ <1 とする。
  ;;d_->string-binary : 0以上1未満の、10進数の小数を2進数に変換する。出力形式は文字列
  (define (d_->string-binary d_)
    (define (d_->string-binary-p d_ tmp count)
      (if (or (= d_ 0) (= count 50)) ;; 小数点以下50桁まで表示
        tmp
        (d_->string-binary-p (- (* d_ 2)
                                (floor (* d_ 2)))
                             (string-append tmp
                                            (number->string (inexact->exact (floor (* d_ 2)))))
                             (+ count 1))))
                
    (cond ((= d_ 0) "")
          ((and (< 0 d_) (< d_ 1)) (d_->string-binary-p d_ "." 0))
          (error "d_ should be a number in [0,1)!!")))
  
  ;;ここから本体
  (let ((d+->string-binary (lambda (x) (string-append (number->string (di->binary (inexact->exact (floor x))))
                                                      (d_->string-binary (- x (floor x)))))))
    (cond ((>= d 0) (d+->string-binary d))
          ((< d 0) (string-append "-"
                                  (d+->string-binary (* -1 d))))
          (else (error "decimal->string-binary: wrong type argument!!")))))

設問を解いてみる。11.3125(10進数)を2進数に変換する。

(decimal->string-binary 11.3125)

;;実行結果
"1011.0101"

うまくいった。

演習問題2.1

次の連立1次方程式をガウスの消去法で解け.


\left\{
\begin{eqnarray}
x_1 + 2 x_2 - 12 x_3 + 8 x_4 &=& 27\\
5 x_1 + 4 x_2 + 7 x_3 - 2 x_4 &=& 4 \\
 -3 x_1 + 7 x_2 + 9 x_3 + 5 x_4 &=& 11 \\
6 x_1 - 12 x_2 - 9 x_3 + 3 x_4 &=& 49
\end{eqnarray}
\right.

解答:
Schemeで正確数を使って計算した結果x1

((811/267) (-22703/10947) (3773/3649) (55444/10947))

Schemeで非正確数を使って計算した結果x2

((3.0374531835206025) (-2.073901525532115) (1.0339819128528416) (5.064766602722216))

最後に誤差を計算してみる。
x2-x1

((3.1086244689504383e-15) (-5.773159728050814e-15) (5.10702591327572e-15) (8.881784197001252e-15))

この例では、誤差はほとんどない。