第6章の練習，投票

台形公式

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

練習

以下の積分を台形公式で近似して円周率を求めよ。n を引数として円周率を返す関数を定義せよ。

trapezoid.rb の f を定義しなおし、

 def pi(n)
   4.0 * trapezoid(....)
 end

進捗状況の確認

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

正しく動いたら、n の値を変えて、真の値に近づく様子を観察しよう。さらに時間があれば次のシンプソン法にチャレンジ。

練習

以下の積分をシンプソン公式で近似して円周率を求めよ。n を引数として円周率を返す関数を定義せよ。

ヒント

def f(x)
  # に対応する式を書く
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) + ..... # 穴埋め
  end
  deltax * sum / 3.0
end

以下を確かめる

f(0.0) -> 1.0
f(0.5) ->  0.8
f(1.0) -> 0.5
4.0*simpson(0.0,1.0,100) -> 3.14159265….

練習

円周率を求める式として以下の式を用いた時の，台形公式，シンプソン公式での収束の様子を観察しなさい．

以下を確かめる

f(0.0) -> 1.0
f(0.5) -> 0.866025...
f(1.0) -> 0.0
4.0 * simpson(0.0,1.0,10)
4.0 * simpson(0.0,1.0,100)
4.0 * simpson(0.0,1.0,1000)
4.0 * simpson(0.0,1.0,10000)
4.0 * simpson(0.0,1.0,100000)

円周率の近似計算

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

def mcplot(n)
  a = make2d(500,500)
  for i in 1..n
    x= rand() # random number in [0,1)
    y= rand()
    if x*x + y*y < 1.0
      a[y*500][x*500] = 1.0
    else
      a[y*500][x*500] = 0.5
    end
  end
  a
end

練習

4*montecarlo(n) を色々な引数で実行し、円周率が近似できることを確かめよ。
show(mcplot(n)) を実行して、点がばらまかれる様子を観察せよ。
- 配列の引数に浮動小数点数を与えると切り捨てられる。

次の数は10進数で何?

0 01111111 10000000000000000000000

(単精度表現)

ヒント
IEEE 754の実数表現では単精度表現32ビットで

s e₁e₂...e₈ f₁f₂...f₂₃

と表現される数は，

(-1)^s × (1.f₁f₂...f₂₃)₍₂₎ × 2^{(e₁e₂...e₈)₍₂₎-127}

である(例外として，指数部と仮数部がすべて0の実数は0になる)．Ruby言語では2進数の整数は先頭に「0b」をつけて表現できるので，

(-1)**(s) * (0b1f₁f₂...f₂₃ * 2**(-23)) * 2**(0be₁e₂...e₈ - 127)

を計算すれば良い．今の例では，

(-1)**(0) * (0b110000000000000000000000 * 2**(-23)) * 2**(0b01111111 - 127)

となる．

次の数は10進数で何?

0 10000000 10000000000000000000000

(単精度表現)

次の数は10進数で何?

0 01111110 00000000000000000000000

(単精度表現)

次の数は10進数で何?

0 00000000 00000000000000000000000

(単精度表現)

ヒント
例外にあたる．

1/3 は 2 進数では？

0.01
0.011111111111..
0.010101010101..
0.010010010010..
0.011011011011..

ヒント
手計算で求めても良いが，Rubyを使って求めることもできる

1.0/3.0 => 0.333333333333333

だが，小数点以下8桁を求めるために，2**8=256倍してみる．

(1.0/3.0)*256 => 85.3333333333333

整数を2進数に直すには，「.to_s(2)」をつける．

85.to_s(2) => "XXXXXXX"

したがって，

(1.0/3.0) = 0.0XXXXXXX....

となる．

丸め誤差

0.1 * 3 == 0.3

0.1 * 3 - 0.3

2 ** -54

桁落ち誤差

include(Math)
def abs (x)
  if x >=0 
    x
  else
    -x
  end
end
def f(x) 
  log(x)
end
def cancellations(x,h_digits) 
  errors=Array.new(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)

情報落ち誤差

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