11/22 第7回 数値解析

前回までの補足

• 講義のスライド(PDF形式)は，cfiveの「情報科学(水4)」のコースの「教材」のところからダウンロード可能にしてある．講義受講者以外への再配布は許可されないので注意すること．

スライドの補足

• 数値解析には汎用のプログラミング言語を使わずに，数式処理系を使う方法もある．教育用計算機システムでは，Mathematica(HWB:Mathematica入門)が使える．
  • Mathematica で微分方程式を解く

10/25の課題について

• 締切は，今日(11/22)の21時. 11/22 14:00 時点で提出者は11名．コメントを付けておいたので今後の参考にするように．
• コメント欄に書ききれなかった質問に対して回答する．

  Q.

  配列の初期化をするには，

  station=Array.new(174)
  station[0]=Subway.new("銀座線","渋谷",0.0,0,"")
  station[1]=Subway.new("銀座線","表参道",1.3,2,"")
  ..
  

  のようにする必要があるのですか?

  A.

  スライドの例ではそのようになっていましたが，

   
  station=[Subway.new("銀座線","渋谷",0.0,0,""),
           Subway.new("銀座線","表参道",1.3,2,""),
           ...
          ]
  

  のような書き方ができます．配列の大きさは station.size で得られます．

  ---

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

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"}