拡散項のみ

• 拡散のシミュレーション
  • このときはランダムウォークから拡散方程式まで
  • このときは解析的に解いたあとでグラフを書いたもの
• 今回は拡散方程式を解析的にではなく数値的に計算していく
• いろいろな項を考えているがまずは拡散項のみ
  • 独立変数：
  • 従属変数：
  • のつもり
• 計算にあたり
  • の値を array U[i,j,t]に入れてある
• 1. 計算について
  • 区間：の内点について計算
  • つまり の範囲で計算
  • のときのarray U の値は境界条件を表すようにする
• 2. 境界条件について
  • Dirichlet、Neumann 両方入れてみた
  • Dirichlet条件
    • 0に固定することに
    • については何も計算しないので
    • となっている
  • Neumann条件にはなっていないかも
    • U[1,,t]=U[2,,t]　としたので
    • にはなっているが
    • 今回の計算の仕方ではfluxが(x=1)のところで0になっている訳ではない？
    • 追記
      • 保存されるようにしたのがこちら
• 3. 初期条件については
  • 今回はとりあえずパルス状に発生させておいた

N<-30
M<-40
T<-50

d<-1    #キョリの次元
dt<-1   #時間の次元
D<-0.12  #拡散係数

U<-array(0,c(N,M,T))

U[2:10,2:10,1]<-1    #初期条件

#拡散
for(t in 1:(T-1)){

for(i in 2:(N-1)){
for(j in 2:(M-1)){
du<-D*(sum(U[(i-1):(i+1),(j-1):(j+1),t])-9*U[i,j,t])/(d^2)*dt        #uの変化量duを計算
U[i,j,t+1]<-U[i,j,t]+du
}
}
U[,1,t+1]<-U[,2,t+1]  #Neumann条件
U[1,,t+1]<-U[2,,t+1]

}
min(U)

for (t in 1:T){
persp(U[,,t],col=3,theta=20,phi=30,zlim=c(0,1))
#image(U[,,t],col=topo.colors(10))
}