第6章の練習，投票

台形公式

include Math
def f(x)
  x/((x+1.0)*(x+2.0))
end

def trapezoid(xs,xe,n)
  deltax = (xe-xs)*1.0/n
  sum = (f(xs)+f(xe))/2.0
  for i in 1..(n-1)
    sum = sum + f(xs+i*deltax)
  end
  deltax * sum
end

def g(x)
  2*log(2+x)-log(1+x)
end

trapezoid(0,1,n)の誤差は?

trapezoid(0,1,n) とg(1)-g(0)の(絶対)誤差がいくつになるか，n=10^mを小さい方から試して，はじめて誤差が10^-8以下になったのは?

100
1000
10000
100000
1000000

練習

Simpson公式に従って関数fを積分するプログラムの穴埋めをして完成させなさい．

load("./trapezoid.rb")
def f(x)
  x/((x+1.0)*(x+2.0))
end

def simpson(xs,xe,n)
  deltax = (xe-xs)*0.5/n
  sum = f(xs)+f(xe)+4*f(xs+deltax)
  for i in 1..(n-1)
    sum = sum + ... * f(xs+2*i*deltax) + 4 * f( ... ) # 穴埋め
  end
  deltax * sum / 3.0
end

進捗状況の確認

simpson(xs,xe,n)ができた時点で投票してください．

simpson(0,1,100)がg(1)-g(0)=0.1177803565... と近い値になっていることを確認してください．

simpson(0,1,n)の誤差は?

simpson(0,1,n) とg(1)-g(0)の(絶対)誤差がいくつになるか，nを小さい方から試して，誤差がはじめて10^-8以下になったのは?

10
21
42
100
1000
10000
100000

Monte Carlo法による積分

def montecarlo(n)
  m=0
  for i in 1..n
    x=rand() # random number in [0,1)
    y=rand()
    if x*x + y*y < 1.0
      m = m + 1
    end
  end
  m*1.0/n
end

montecarlo(n)の誤差は?

montecarlo(n)とπ/4 = 0.785398163397448…の(絶対)誤差がいくつになるか， n=10^mを小さい方から試して，誤差がはじめて10^-3以下になったのは?

100
1000
10000
100000
1000000

乱数を使っているので他の人と違う答えが出ることはある．

次の数は10進数で何?

0 01111110 00000000000000000000000

(単精度表現)

ヒント

["00111111000000000000000000000000".to_i(2)].pack("L").unpack("e")

丸め誤差

0.1 * 3 == 0.3

0.1 * 3 - 0.3

2 ** -54

桁落ち誤差

def abs(x)
  if x < 0
    -x
  else
    x
  end
end
def f(x)
    log(x)
end
def cancellations(x,h_digits)
    errors=make1d(h_digits+1)
    for i in 0..h_digits
        h = 0.1**i
        df=(f(x+h)-f(x))/h
        errors[i] = abs(df-1.0/x)
    end
    errors
end

cancellations(2.0,15)

式変形による桁落ちの回避

include Math
def f(x)
 1-sqrt(1-x)
end
def f1(x)
 x/(1+sqrt(1-x))
end

irb(main):055:0> require 'bigdecimal'
irb(main):056:0> BigDecimal::new("1.0",100)-BigDecimal::new("0.99999999",100).sqrt(100)
=> #
irb(main):024:0> f(1e-8)
=> 5.00000008063495e-09
irb(main):025:0> f1(1e-8)
=> 5.0000000125e-09

式変形による桁落ちの回避(12/13追加分)

12/13 追加分のプログラムは avoid_error.rbからダウンロード可能にしておきます．
式変形による桁落ちの回避の別の例を紹介する．絶対値の小さいxに対して， \[g(x) = 1 - \frac{1}{\sqrt{1+x}} = \frac{\sqrt{1+x}-1}{\sqrt{1+x}} \] を計算する．x=0.0000001に対して正しい値は，およそ　

=> 4.9999996250000312e-08

となる．これを

include Math
def g1(x)
  1-1/sqrt(1+x)
end
def g2(x)
  (sqrt(1+x)-1)/sqrt(1+x)
end

のどちらで計算しても，

irb(main):009:0* g1(0.0000001)
=> 4.9999996254435075e-08
irb(main):010:0> g2(0.0000001)
=> 4.9999996307948274e-08

と，かなり誤差を含んでいることがわかる． g1では「1」と「1/sqrt(1+x)」というほぼ等しい値の減算が生じているし，g2でも「sqrt(1+x)」と「1」の減算が生じている．
g(x) の分子と分母に\(\sqrt{1+x}+1\)をかけて， \[ \begin{eqnarray} g(x) & = & \frac{(\sqrt{1+x}-1) * (\sqrt{1+x}+1)}{\sqrt{1+x} * (\sqrt{1+x}-1)} \\ & = & \frac{x}{(1+x)+\sqrt{1+x}} \end{eqnarray} \] と変形した関数g3(x)

include Math
def g3(x)
  x/(1+x+sqrt(1+x))
end

を使うと

irb(main):036:0> g3(0.0000001)
=> 4.99999962500003e-08

となり，誤差がほとんどないことがわかる．

表現範囲を超えるエラーの回避(12/13追加分)

複素数の a + b iの絶対値は教科書通りに実行すると，\( \sqrt{a^2 + b^2}\) となる．

include(Math)
def abs_complex(a, b)
  return sqrt(a**2 + b**2)
end

これを使って，\(3 * 10^{200} + 4 * 10^{200} i \) の絶対値を計算すると，計算結果は \(5 * 10^{200}\)で倍精度表現で表現できる範囲(\(1.7 * 10^{308}\)よりも小さい)なのに結果は以下のようにInfinity (正の無限大)を表す値になってしまう．

irb(main):007:0> abs_complex(3e200,4e200)
=> Infinity

a**2を計算した時点で\(9 * 10^{400}\)となるが，この時点で倍精度表現で表現できる範囲を超えてしまっているから，それ以降の演算は正しくおこなえない．
これを避けるには以下のように書く．

def abs_complex_r(a, b)
  s = max(abs(a), abs(b))
  if b == 0
    s
  else
    x = a / s
    y = b / s
    s * sqrt(x**2 + y**2)
  end
end

同じ計算をしてみると，

irb(main):019:0> abs_complex_r(3e200, 4e200)
=> 4.9999999999999995e+200

と正しく計算される．

情報落ち誤差

# simulate 32 bit floating point representaiton 
def coerce32(f)
  [f].pack("f").unpack("f")[0] 
end

def sum(digits) 
  n=10**digits 
  sum = 0.0 
  d = coerce32(1.0 / n) 
  for i in 1..n
    sum = coerce32(sum + d) 
  end
  sum 
end

sum(6)
sum(7)
sum(8)

誤差の回避(12/13追加分)

同じ位の大きさのからなる長いリストの総和を計算する時は，誤差を回避するために，同じ位の長さのリスト2つに分割してそれぞれの総和(同じ位の長さのリストの総和は同程度の桁になることが期待される)を求めて足し合わせるというテクニックが実際のプログラムでも良く用いられている．これを見てみる．以下の関数make_listは 1.00, 1.01, 1.02, ..., 1.09 を繰り返した長さnのリストを作成する．

load("./make1d.rb")

def make_list(n)
  a = make1d(n)
  for i in 0..n-1
    a[i] = 1.0 + 0.01 * (i % 10)
  end
  a
end

これを普通に足し合わせてみる．

def sum_array(a)
  r = 0
  for i in 0..a.length()-1
    r = r + a[i]
  end
  r
end

これに n = 10000000 (1e7) で計算してみると，

irb(main):047:0> sum_array(make_list(10000000))
=> 10449999.99945153

となり，有効数字が11桁程度まで落ちてしまっていることがわかる．以下のように，長さが2以上の時には半分に分けて計算する関数を書いてみる．

def sum_rec(a, low, high)
  if high - low == 1
    a[low]
  else
    mid = (high + low) / 2
    sum_rec(a, low, mid) + sum_rec(a, mid, high)
  end
end

def sum_array_r(a)
  if a.length() == 0
    return 0
  end
  return sum_rec(a, 0, a.length())
end

今度は，

irb(main):049:0> sum_array_r(make_list(10000000))
=> 10450000.0

となり，誤差が微小(ある可能性もあるが，表示できる範囲では区別がつかない)だとわかる．「同じような数の和」という制約を外して適用可能なKahan summation algorithmというアルゴリズムもある．

誤差一般に関する補足(12/13追加)

「誤差」という言葉から「測定誤差」を思い浮かべる人も多いと思う．その類推で，「桁落ち」を「同じ計算をしても毎回答えが違う」と思ってしまった人もいるかもしれないが，そのようなことはない．計算機は与えられた入力に対して，最善の努力で計算結果を表現可能な範囲で計算するだけである．

しかし，「桁落ち」，「情報落ち」等の名前をつけるのは，計算中でポイントとなる点を探して，それを避ける方法がある場合に全体の計算結果を改善できるケースがあるからである．

たとえば，同じような数 x, yがあった時に x-y を計算して結果を得る必要がある場合，得られた結果の有効桁数が落ちてしまうのは当然であり，「桁落ち」が生じてしまうは回避できない．しかし，全体の計算の一箇所で「桁落ち」が生じる場所が分かるならば，それを回避できるケースがあるのは例で見たとおりである．

「情報落ち」について

Gauss-Jordan法

def gj(a) # Gauss-Jordan method without pivoting
  row = a.length()
  col = a[0].length()
  for k in 0..(col-2)
    akk = a[k][k]
    for i in 0..(col-1)       # normalize row k
      a[k][i]=a[k][i]*1.0/akk
    end
    for i in 0..(row-1)      # eliminate column k
      if i != k              # of all rows but k
    aik = a[i][k]
    for j in k..(col-1)
          a[i][j] = a[i][j] - aik * a[k][j]
    end
      end
    end
  end
  a
end

a=[[1,1,-1,2],[3,5,-7,0],[2,-3,1,5]]
gj(a)

Gauss-Jordan法(Pivoting)

ファイル abs.rb を新しく作る．

def abs (x)
  if x >=0 
    x
  else
    -x
  end
end

#a[i][k]の絶対値が最大となるi>=kを探す
def maxrow(a,k)
  max_i=k
  for i in (k+1)..(a.length()-1)
    if abs(a[max_i][k]) < abs(a[i][k])
      # 穴埋め
    end
  end
  max_i
end

def swap(a,i,j)
  tmp=a[i]
  a[i]=a[j]
  a[j]=tmp
end
def gjp(a) # Gauss-Jordan method WITH pivoting
  row = a.length()
  col = a[0].length()
  for k in 0..(col-2)
    max=maxrow(a,k) # find absolute-maximal coeff. 
    swap(a,k,max)   # swap rows
    akk = a[k][k]
    for i in 0..(col-1)       # normalize row k
      a[k][i]=a[k][i]*1.0/akk
    end
    for i in 0..(row-1)          # eliminate column k
      if i != k          # of all rows but k
    aik = a[i][k]
    for j in k..(col-1)
      a[i][j] = a[i][j] - aik * a[k][j]
    end
      end
    end
  end
  a
end

b=[[1,-50,-3,-90],[-85,2,-25,-6],[79,5,30,-1]]
gj(b)

一度解を求めると b が書き変わってしまうので，もう一度解く前にbを設定し直す必要がある．

b=[[1,-50,-3,-90],[-85,2,-25,-6],[79,5,30,-1]]

maxrow(b,0)
=> 1
maxrow(b,1)
=> 2

gjp(b)

練習

以下の連立1次方程式をGauss-Jordan法で解いてみよ。Pivotingするかしないかで、結果が変わるか観察せよ。

進捗状況の確認

maxrowができた時点で投票してください
値がおかしい
プログラムが動かない