前の例は簡単すぎて面白くないので、その次の課題のやつをやってみた。拡散方程式。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() }
まっくろくろすけみたいになってたり、軸の方向が明後日向いてたりするのは仕様です。