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

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

[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）乗

補数，とは，元の数と補数を足した時に桁が上がる数のうち最小のものをいう．試しに，上の表をみて1(001)と1に対する2の補数（111)を足してみよ．111+001なので，桁が一つあがって1000となる．このような数を（2進数の場合）1に対する2の補数という．

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

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

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

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

[1.2] 実数の内部表現

単精度実数は，実数を32個の0と1で表現する．

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

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

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

[2] 関数副プログラム

フォートランにはデフォルトでいくつかの関数が用意されているが，ユーザー が独自に関数を定義して用いたい場合は「関数副プログラム」というものをつかって関数 を定義することが出来る．使い方は


t function f(a1,a2,....aN)
......
f=.....
return
end
function文 f=.... 関数副プログラムの定義をおこなう
t 型宣言：double precision, integerなど
a1,a2,...aN 仮引数：型宣言をすること
プログラム中で用いる時は, 

b=f(c1,c2,....cn)
のように引数を代入して行う．

たとえば, 1からnまでの総和を求める関数副プログラムは 

     
      integer function souwa(n)
      implicit none
      integer i,n

      souwa = 0  !  ここで(n)はつかない
      do i = 1, n
         souwa = souwa + i   ! ここでも(n)はつかない
      enddo
      
      end function

      program sum
      implicit none
      integer souwa,n  ! souwaを宣言する必要あり．(n)はつかない
      read(*,*) n
      write(*,*) souwa(n)  ! ここではnを代入するので(n)がつく
      end program sum

注意点

[3] 台形公式

関数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，．．．と変化させ，誤差がどう変化するか観察せよ．

[4] シンプソン公式

シンプソン公式は，台形公式の精度をあげたものである．台形公式の精度を計算すると，幅Δ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分割したときの誤差がほぼ同程度となり，シンプソン則の効率のよさが分かる．

[5] オイラー法

一階の常微分方程式


\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)

と比較せよ．積分区間(2.D0)を分割する数nを20, 200, 2000と変化させ，（時間刻みdtはそれぞれ1.D-1, 1.D-2, 1.D-3となる）誤差の変化を観察せよ．

[5] 今日の宿題（1点）

練習問題4のプログラムを結果（３とおり）と共に送ること. 〆切は12/12.


日程表へ戻る <<