11/22 第7回 数値解析


前回までの補足


スライドの補足


10/25の課題について


スライドに使われたプログラム

Jacobi法の計算例

ma=[[5.0,0,1.0],[3.0,4.0,0.0],[1.0,1.0,3.0]]
b=[3,2,1]
mb=[[0,0,0],[0,0,0],[0,0,0]]
mc=[[0,0,0],[0,0,0],[0,0,0]]
md=[[0,0,0],[0,0,0],[0,0,0]]
mdi=[[0,0,0],[0,0,0],[0,0,0]]
c=[0,0,0]
for i in 0..2
  for j in 0..2
    if (i==j)
      md[i][i]=ma[i][i]
      mdi[i][i]=1.0/ma[i][i]
    else
      md[i][j]=mdi[i][j]=0
    end
    mc[i][j]=ma[i][j]-md[i][j]
  end
end
for i in 0..2
  for j in 0..2
    for k in 0..2
      mb[i][j]+=(-mdi[i][k]*mc[k][j])
    end
    c[i]+=(mdi[i][j]*b[j])
  end
end
iter=15
xs=[0.5763,0.0678,0.1186]
x=[0,0,0]
iter.times{ 
  nx=[0,0,0]
  for i in 0..2
    for j in 0..2
      nx[i]+=(mb[i][j]*x[j])
    end
    nx[i]+=c[i]
  end
  x=nx
  err=Math.sqrt((x[0]-xs[0])**2+(x[1]-xs[1])**2+(x[2]-xs[2])**2)
  print "k= ",k,"\t","x=[ ",x[0],",",x[1],",",x[2],"] err= ", err, "\n"
}

円周率の近似計算

n=1000000
m=0
n.times{
  x=rand
  y=rand
  if (x*x+y*y)<1.0 
    m+=1
  end
}
puts 4*Float(m)/n

Random Walk

m=30000
n=100
hist=Array.new(n*2,0)
m.times{
  x=0
  n.times{ |i|
    p=rand
    if (p>0.5)
      x+=1
    else
      x-=1
    end
  }
  hist[n+x]+=1
}
hist.size.times{ |i| print i-n,"\t",hist[i],"\n"}