数学が終わっている私は、泣きそうになりながら解いてます。
練習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
練習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 （係数行列の作成）
練習6.12 （係数行列と定数ベクトルの合成）
練習6.13 （画像の復元）
for i in 0..(l-1)
a[j][i]=solved[l*j+i][l**2]
end
end
a
end
a) def trapezoid_sinlog(xs,xe,n)
練習6.15 （Monte Carlo法による３次元の積分）
練習6.16 （Monte Carlo法による積分）
include(Math)
練習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)
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[j][l**2]=s[j/l][j%l]
end
a
end
end
練習6.13 （画像の復元）
def equations_to_image(l,solved)
a=make2d(l,l)
for j in 0..(l-1)a=make2d(l,l)
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
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
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
誤差は自分で比べてみましょう。 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
(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
