[教]練習問題（第六章）

教科書「第六章・数値計算」の練習問題の解答例です。


数学が終わっている私は、泣きそうになりながら解いてます。


練習6.1 （台形公式の正確さ）
略




練習6.2 （Simpson公式による積分）
a) def  simpson(xs,xe,n)
       deltax=(xe-xs)/(2.0*n)
       sum=f(xs)+4*f(xs+deltax)+f(xe)
       for  i  in 1..(n-1)
         sum=sum+2*(xs+2*i*deltax)+4*f(xs+(2*i+1)*deltax)
       end
       deltax*sum/3.0
     end

b) 略

練習6.3 （Monte Carlo法の正確さ）

略

練習6.4 （乱数列となる条件）　※ここの解答は今までで一番適当です

a) 地域によって、夏や冬などある程度規則があるので、言えない、と思う＞＜

b) 言える。言えなきゃ宝くじにならないぜ！

c) 各数字が4枚までと決まっているので、乱数とは言えない。

練習6.5 （浮動小数点数の表現）

a) ※(2)は、二進数であることを示す

　符号 s = 1

　指数 d'1..d'8 = 11111110

　仮数 d1..d23 = 11000000000000000000000

また、単精度のバイアス値 b = 127 より

　(-1)^s × 1.d1..d23(2) × 2^(d'1..d'8(2) - b)

 = (-1)^1 × 1.11000000000000000000000(2) × 2^(11111110(2) - 127)

 = (-1) × { 2^0 + 2^(-1) + 2^(-2) } × 2^(254 - 127)

 = -2.977471 × 10^38

b) (a) s = 0 , d'1..d'8 = 01111111, d1..d23 = 00000000000000000000000

     (b) 3145728 = 2^20 × 3 = (-1)^0 × 1.5 × 2^21 = (-1)^0 × 1.5 × 2^(148 - 127)　より

           s = 0 , d'1..d'8 = 10010100 , d1..d23 = 10000000000000000000000

     (c) 5/65536 = 5 / 2^16 = (-1)^0 × 1.25 × 2^(-14) = (-1)^0 × 1.25 × 2^(113 - 127)　より

           s = 0 , d'1..d'8 = 01110001, d1..d23 = 01000000000000000000000

練習6.6 （不動小数点の大小）

えっ

符号のとこ見るといっぱつで反証になっちゃうんだけど

x < 0 < y のとき

sx = 1 かつ sy = 0 であるから

(xの浮動小数点数) > (yの不動小数点数）

・・・あれー？

一応符号を無視してやってみます。
でも、私数学苦手ですから。もっとすっきり証明もとむ。

十進数で表されたx,y（0<x<y）はそれぞれ
　x = p × 2^u
　y = q × 2^v
とできる。ただし、
　1 ≦ p,q < 2　かつ　-127 ≦ u,v ≦ 128　ー※
である。
I) u=vと仮定すると（つまり指数部が一致すると仮定）
　x<yより、
　　p<q ⇔ p-1<q-1 ⇔ dx1..dx23(2)<dy1..dy23(2)
　であるから、
　（xの浮動小数点数）<（yの浮動小数点数）
II) u>vと仮定すると
　　y = q × 2^v < 2^(v+1) ≦ 2^u < q × 2^v = x　（∵※）
　よって矛盾。
III) u<vと仮定すると
　　dx'1..dx'8(2) < dy'1..dy'8(2)
　であるから、

　（xの浮動小数点数）<（yの浮動小数点数）

以上より、証明できた。




練習6.7 （丸め誤差）
0.1*5-0.5 = 0.0
となる。これは、　　※(36)は36桁分の省略があることを示す
　0.1(10) = 1.100110011..(36)..0011010(2) × 2^(-4)(10)
　　5(10) = 1.010000000..(36)..0000000(2) × 2^2(10)
　0.5(10) = 1.000000000..(36)..0000000(2) × 2^(-1)(10)
であることを考えると、
　0.1(10) × 5(10) = 10.00000000000..(36)..0000010(2) × 2^(-2)(10)
　　　　　　　　　　　　　　　　　　　　　　　　　桁落ち
　　　　　　　　　 = 1.000000000..(36)..0000000(2) × 2^(-1)(10)
となり、0.5(10)と一致することから成り立つ。




練習6.8 （倍精度表現の精度）
　　仮数部のビット数：52
10進数有効桁数　(約)：log10{2^(52+1)} ≒ 16
　　指数部のビット数：11
　　　　　最小の指数：2^(-1022) ≒ 6.6 × 10^309
　　　　　最大の指数：2^(1023) ≒ 8.9 × 10^307
　　　　　　　最大値：2 × 8.9 × 10^307 ≒ 1.8 × 10^308




練習6.9 （Pivotingの完成）
load("./abs.rb")
def  maxrow(a,k)
  s=k

  t=abs(a[k][k])

  if  k!=a.length()-1

    for  y  in k..(a.length()-2)

      if  t < abs(a[y+1][k])

        s=y+1

        t=abs(a[y+1][k])

      end

    end

  end

  s

end

gjp(a)は教科書P137にあるので省略しますよ。

練習6.10 （部分の通過距離）
load("./abs.rb")
include(Math)
def  penetration(theta,d,x,y)
  q = abs(x*cos(theta) + y*sin(theta) - d)
  if  q<0.5
    sqrt(1-4*q**2)
  else
    0.0
  end
end


練習6.11 （係数行列の作成）

def  coefficients(l,n)

  m = make2d(l*n,l**2)

  for  k  in 0..(n-1)

    for  d  in  0..(l-1)

      for  y  in  0..(l-1)

        for  x  in  0..(l-1)

          m[l*k+d][l*y+x] = penetration(3.1415926535897932384*k/n,d,x,y)

        end

      end

    end

  end

  m

end

練習6.12 （係数行列と定数ベクトルの合成）

def  make_equations(l,n,m,s)
  a=make2d(l*n,l**2+1)
  for j in 0..(l*n-1)

    for i in 0..(l**2-1)

      a[j][i]=m[j][i]

    end
    a[j][l**2]=s[j/l][j%l]
  end

  a
end

練習6.13 （画像の復元）

def  equations_to_image(l,solved)
  a=make2d(l,l)

  for j in 0..(l-1)
    for i in 0..(l-1)
      a[j][i]=solved[l*j+i][l**2]
    end
  end
  a
end

練習6.14 （解析的な積分が難しい関数の数値積分）

まずはg(x)を定義。ちなみにsinやlogが出てくるので、数学関数を含むことを宣言しましょう。

Include(Math)

def g(x)

  sin(x)*1.0/log(x)

end

a) def trapezoid_sinlog(xs,xe,n)

      deltax = (xe-xs)*1.0/n

      sum = (g(xs)+g(xe))/2.0

      for i in 1..(n-1)

        sum = sum + g(xs+i*deltax)

      end

      deltax * sum

    end

b) def simpson_sinlog(xs,xe,n)

      deltax = (xe-xs)/(2.0*n)

      sum = g(xs)+4*g(xs+deltax)+g(xe)

      for i in 0..(n-1)

        sum = sum + 2*g(xs+2*i*deltax)+4*g(xs+(2*i+1)*deltax)

      end

      deltax*sum/3.0

    end

練習6.15 （Monte Carlo法による３次元の積分）

include(Math)　　　←乱数rand()は数学関数に含まれる
def  montecarlo3d(n)
  m=0
    for  i  in 1..n
      x=rand()
      y=rand()
      z=rand()
      if  x*x+y*y+z*z<1.0
        m=m+1
      end
    end
  m*8.0/n
end

誤差は自分で比べてみましょう。


練習6.16 （Monte Carlo法による積分）
include(Math) 

def f(xs,x)
  (x+xs)*1.0/((x+xs+1)*(x+xs+2))
end

def montecarlo_f(xs,xe,n)
  m=0
    for i in 1..n
      x=rand()*(xe-xs)
      y=rand()*(3-2*sqrt(2))
      if y<f(xs,x)
        m=m+1
      end
    end
  m*(xe-xs)*(3-2*sqrt(2))/n
end

[ス]練習問題（第六章）

授業スライド「06数値計算」の練習問題の解答例です。 

えーーーーーーー
皆さん今日は。S23-7 情報科学シケタイです。
突然ですが私数学がとっても苦手でして、この章に突入してからあわあわな感じです。
それを考慮の上ご覧ください＞＜

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

trapezoid.rb の f を定義しなおし、
def pi(n)
4.0 * trapezoid(....)
end

def  f(x)
  1.0/(1+x**2)
end


def  pi(n)
  4.0*trapezoid(0,1,n)
end

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

def  f(x)
  4.0/(1+x**2)
end


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

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

略

[管理情報] 教科書第四章の解答について

教科書第四章の練習問題の解答で重大なミスが発見されました。

練習4.6、練習4.7で、「繰り返しを用いて」と指定されているにも関わらず、

再帰を用いた解答を掲載していました。

今現在訂正を急いでおります。


12/04 0:12　練習4.6を訂正いたしました。
12/04 0:47　練習4.7を訂正いたしました。

[基]ターミナルへの命令

Rubyを起動する前、ターミナルへの命令をいくつかあげます。

pwd　　　　今プログラムが働いている階層（フォルダ）を確認する

cd ruby　　今いる階層内のフォルダ"ruby"に移動する

cd ~　　　　最初にプログラムが働いていた階層に戻る

is　　　　　今いる階層内のファイル名を表示する

cat bmi.rb　今いる階層内のファイル"bmi.rb"の内容を表示する

覚えましょう。