第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=10^mの形のものに限る)は?

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)の誤差が10^-8以下になる最小のnは?

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..の誤差が10^-3以下になる最小のn(ただしn=10^mの形のものに限る)は?

100
1000
10000
100000
1000000

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

次の数は10進数で何?

0 01111110 00000000000000000000000

(単精度表現)

丸め誤差

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

進捗状況の確認

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