第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)の誤差が10-8以下になる最小のn(ただしn=10mの形のものに限る)は?
  1. 100
  2. 1000
  3. 10000
  4. 100000
  5. 1000000

練習

Simpson公式に従って関数fを積分するプログラムの穴埋めをして完成させなさい.
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

進捗状況の確認

  1. 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)の誤差が10-8以下になる最小のnは?
  1. 10
  2. 21
  3. 42
  4. 100
  5. 1000
  6. 10000
  7. 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..の誤差が10-3以下になる最小のn(ただしn=10mの形のものに限る)は?
  1. 100
  2. 1000
  3. 10000
  4. 100000
  5. 1000000
乱数を使っているので他の人と違う答えが出ることはある.

次の数は10進数で何?

0 01111110 00000000000000000000000
(単精度表現)
  1. 0.0
  2. 0.1
  3. 0.2
  4. 0.3
  5. 0.5
  6. 1.0
  7. 2.0
  8. 3.0
  9. 5.0

丸め誤差

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

情報落ち誤差

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

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)

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するかしないかで、結果が変わるか観察せよ。