数値解析法及び演習 第八回

誤差と数値積分
  1. 数の表現
  2. 台形公式
  3. シンプソン公式
  4. オイラー法
  5. 今日の宿題

[1] 数の表現

[1.1] 整数の内部表現

練習問題1

計算機はミクロにはon,offの2状態しか扱うことができない. したがって計算機が計算している「数」は実際にはすべて2進数である. そこで計算機の内部で数がどう表されているか考えてみよう. ここではまず整数を考える.

簡単のために3ビットの計算機を考えよう. つまりこの計算機で取り扱えるのは3個の0, 1の組である. 3ビット計算機は次の表に表されるように整数を取り扱う.

整数の内部表現（3ビット計算機）
符合のための1bit 数値のための2bits 表される整数
0 11 3
0 10 2
0 01 1
0 00 0
1 11 -1（1に対する2の補数）
1 10 -2（2に対する2の補数）
1 01 -3（3に対する2の補数）
1 00 -4=2の（3-1）乗

この表から推測されるように, nビット計算機で取り扱える数は2の(n-1)乗-1からマイナス2の(n-1)乗までとなる. 次のプログラムを実行してみよ. 

      implicit none
      integer n,m,k
      n=2147483647
      m=n+1
      k=m+n+1
      write(6,*) n,m,k
      end

このプログラムはどうなるであろうか？ 

      implicit none
      integer n,m,k
      n=2147483648
      m=n+1
      k=m+n+1
      write(6,*) n,m,k
      end

[1.2] 実数の内部表現

現在流通している多くのコンピューターは32ビット計算機である．これは，内部で扱う数値を32個の0と1で表現しているということである．

単精度実数の内部表現（32ビット計算機）
仮数部の符号1bit 指数部の8bits 仮数部（23bits）
0 00000000 00000000000000000000000
倍精度実数の内部表現（32ビット計算機）
仮数部の符号1bit 指数部の11bits 仮数部（52bits）
0 00000000000 00000000000000000000000..

1bitはプラスかマイナスの符号につかう．8bitsを指数部に使う（倍精度の場合11bits）．8個の0と1で-126から＋127までの数を表す（倍精度の場合-1022から1023まで）. ここからフォートランで用いることができる数値は2の127乗=10の38乗まで（倍精度の場合10の308乗まで）ということになる．広く採用されているIEEE規格では基数を2ととることに注意．

残りの23個（倍精度の場合52個）を使って仮数部(0.234E5 だったら 0.234の部分)を表す．小数第1桁目は1ときまっているので表現はしないため実質2進数で24桁（倍精度の場合53桁）ということになる. 24桁なので仮数部で表せるもっとも小さな数は2の-24乗=0.6*10^(-7)（倍精度の場合0.1*10^(-15)）となる．このことから数値的な誤差が発生する． 次のプログラムを実行させてみよう． 
      implicit none
      real*8 a,b,c,d,e
      a=1.D-15
      b=1.D-16
      c=1.D0
      d=c+a
      e=c+b
      write(6,"(3E26.18)") a,b,c,d,e
      end
結果はどうなったであろうか? これが誤差を生み出す原因である．ここで見ら れる誤差は二つあって である．

[2] 台形公式

関数f(x)のある区間[x0, xn]の積分を考える．この区間をn個に分割すると積分は


\begin{displaymath} I=\int_{x_0}^{x_n}f(x)dx=\sum_{i=0}^{n-1}\int_{x_i}^{x_{i+1}}f(x)dx \end{displaymath} (3)

となる．ここでf(x)よりも下の面積を台形で近似すると


\begin{displaymath} \int_{x_i}^{x_{i+1}}f(x)dx={h\over 2}(f(x_i)+f(x_{i+1})) \end{displaymath} (4)

となる．ここから, 求めたかった積分は次のように書ける.


$\displaystyle I$ $\textstyle =$ $\displaystyle {h\over 2}(f(x_0)+2f(x_1)+2f(x_2)+...+2f(x_{n-1})+f(x_n))$  
  $\textstyle =$ $\displaystyle \sum_{i=0}^{n-1}{h\over 2}(f(x_i)+f(x_{i+1}))$ (5)

練習問題2


\begin{displaymath} I=\int_{0}^{1}e^xdx=e-1=1.718281828... \end{displaymath} (6)

を台形公式により求めよ．区間の分割の数nを10，100，1000，．．．と変化させ，誤差がどう変化するか観察せよ．またなぜそのnで誤差が変化するのか考察せよ．

[3] シンプソン公式

シンプソン公式は，台形公式の精度をあげたものである．台形公式の精度を計算すると，幅Δxの2乗に比例することを示すことができる．
台形公式ではf(x)上の二つの点を直線でむすんだ．よりよく数値積分を行うためには次の精度，つまりf(x)を二次曲線で近似すればよい．台形公式では2 点を用いることで直線を決定した．なので今度は3点を用いれば二次曲線を決定することができる．シンプソン公式は, 誤差がΔxの4乗に比例する.
並んだ二つの区間x1-x2,x2-x3を考える．それぞれ関数fの値はf1,f2,f3とする．この区間において関数fが二次式


\begin{displaymath} f(x)=ax^2+bx+c \end{displaymath} (8)

であったとすると，この区間の積分は次になることを示すことができる．


\begin{displaymath} \int_{x_1}^{x_3}f(x)dx={h\over 3}(f_1+4f_2+f_3) \end{displaymath} (9)

この等式を示すにはまず, 式（8）の積分を行なうと


$\displaystyle \int_{x_1}^{x_3}f(x)dx$ $\textstyle =$ $\displaystyle {a\over 3}(x_3^3-x_1^3)+{b\over 2}(x_3^2-x_1^2)+c(x_3-x_ 1)$  
  $\textstyle =$ $\displaystyle (x_3-x_1)\left[{a\over 3}(x_1^2+x_1x_3+x_3^2)+{b\over 2}(x_1+x_3)+c\right]$ (10)

一方, 式（9）の右辺を計算すると


$\displaystyle {h\over 3}(f_1+4f_2+f_3)$ $\textstyle =$ $\displaystyle a(x_1^2+4x_2^2+x_3^2)+b(x_1+4x_2+x_3)+c(1+4+1)$  
  $\textstyle =$ $\displaystyle 2h[{a\over 3}(x_1^2+x_1x_3+x_3^2)+{b\over 2}(x_1+x_3)+c]$ (11)

x3-x1=2hであるのでこれら二式は等しいことがわかる.

練習問題3

シンプソン公式を用い，先ほど台形公式を用いて計算した積分を行い，台形公式の場合と結果を比較せよ．積分区間を10, 100, 1000等分して誤差を比較すること．台形則の誤差は刻みの数の二乗に反比例し，シンプソン則は4乗に反比例する．台形則で1000分割したときの誤差とシンプソン則で10分割したときの誤差がほぼ同程度となり，シンプソン則の効率のよさが分かる．

[4] オイラー法

一階の常微分方程式


\begin{displaymath} {dy\over dt}=f(t,y) \end{displaymath} (1)

を初期条件 t=t0 のとき y=y0 のもとで解くことを考える．tをtn=t0+n△tと△tづつ増加させるとき，n回目の値 yn=y(tn)を順に求めていく．
最も単純な方法は，


\begin{displaymath} y_{n+1}=y_n+\Delta t y_n'=y_n+\Delta t f(t_n,y_n) \end{displaymath} (2)

としてどんどん足し算していく方法である．これをオイラー法と呼ぶ．

練習問題4

次の常微分方程式
\begin{displaymath} {dy\over dt}=-y \end{displaymath} (3)

を初期条件t=0のときy=1のもとでt=2まで解き，解析的に得られる解


\begin{displaymath} y=\exp(-t) \end{displaymath} (4)

と比較せよ．時間刻みを変化させて(1.D-1, 1.D-2, 1.D-3)誤差の変化を観察せよ．

[5] 今日の宿題

練習問題4のプログラムを結果と共に送ること. 〆切は11/25.


日程表へ戻る <<