数値計算,洲之内著,石渡改訂,サイエンス社を読む (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

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

演習問題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つのパターンに分かれる。

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

そして、いずれの場合も
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)!!")))

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

有限桁で終わらない場合（つまり循環小数の場合）があるので、対策が必要。安直だが、カウンターを用意して途中で打ち切ることにした。上のコードでは小数点以下50桁で打ち切る。
小数部の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))

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