Rでカオス2

前の例は簡単すぎて面白くないので、その次の課題のやつをやってみた。拡散方程式。Rのインデックスが1から始まるとかいらんところではまった。

diffusion <- function(time,sep){
   #拡散方程式がどうなっているかを調べるための関数
   h <- 0.1
   k <- 0.001
   #とりあえずこれで与えられているとする
   diff <- array(dim=c((2.0/sep)+1,time))
   #初期値とかを与えとく
   diff[1,] <- 0
   diff[2/sep+1,] <- 0
   for(i in 2:(2/sep)){  #(2/sep)ってやっとかないととびかたがでかくなっちゃう
       diff[i,1] <- cos(pi/2*((i-1)*sep-1))  #対象にしておくためにインデックスをあわせておく
   }
   for(i in 2:time){  #
       for(j in 2:(2/sep)){  #縦
           diff[j,i] <- k/(h^2)*(diff[j-1,i-1]+diff[j+1,i-1])+(1-2*k/(h^2))*diff[j,i-1]
       }
   }
   return(diff)
}

difplot <- function(time,sep){
    png("difplot.png")
    x <- 1:((2/sep)+1)
    y <- 1:time
    persp(x,y,diffusion(time,sep),theta=120,main="拡散方程式")
    dev.off()
}


まっくろくろすけみたいになってたり、軸の方向が明後日向いてたりするのは仕様です。