スケールファクタのシミュレーション　バージョン1

バージョン1は、単一のパラメータの組み合わせについて、シミュレーション結果をグラフ化する

数値積分のためのパッケージ導入

install.packages(“cubature”, dep=T) #パッケージインストール library(cubature) #パッケージの呼び出し

密度パラメータの定義

R<-5*10^(-5)　　　　　　　　　　　　　 #輻射の密度パラメータ M<-0.30 #全物質（バリオン+ダークマター）の密度パラメータ L<-0.70 #宇宙定数の密度パラメータ O<-R+M+L #曲率のパラメータK=1-O=1-(R+M+L)

数値積分の実行

a<-seq(0,3,0.01) 　　　　 #スケールファクタを0.01刻みで0から3まで定義 B<-seq(0,3,0.01) 　　　　 #数値積分の積分値を格納する箱 HT<-seq(0,3,0.01) 　　 #数値積分結果をNumericとして格納する箱 dHt<-seq(0,3,0.01) #現在値との差分を格納する箱 f<-function(x) x/(R+Mx+Lx^4+(1-O)x2)(1/2)　#被積分関数の定義 for(i in 1:301){B[i]=adaptIntegrate(f,0,a[i])[1]}　#301個のaに対して数値積分しBに格納 HT<-as.numeric(B) #Numericに変換して格納 HT0<-adaptIntegrate(f,0,1)[1] #現在のHt_0を計算 HT0<-as.numeric(HT0) #数値化 for(i in 1:301){dHt[i]<-HT[i]-HT0} #現在値との差分を格納

プロット H(t-t0)を横軸、aとしたグラフを作成

plot(0,0,type=“n”,xlim=c(min(dHt),max(dHt)), ylim=c(0,3), xlab=“H(t-t0)”, ylab=“a”)　#枠作成 lines(dHt,a,lty=1,col=1) #プロット

スケールファクタのシミュレーション バージョン2　

バージョン2は、密度パラメータの組み合わせに対して、シミュレーション結果をグラフ化することが可能

数値積分のためのパッケージ導入

install.packages(“cubature”, dep=T) #パッケージインストール library(cubature) #パッケージの呼び出し

密度パラメータの定義

BR<-c(510^(-5),510(-5),510^(-5),510(-5))　#輻射の密度パラメータベクトル BM<-c(1,0.30,0.30,0.3) #全物質（バリオン+ダークマター）の密度パラメータベクトル BL<-c(0,0.70,1.3,1.7) #宇宙定数の密度パラメータベクトル BO<-BR+BM+BL #曲率のパラメータK=1-O=1-(R+M+L)ベクトル

数値積分の実行

スケールファクタの範囲、計算結果格納ボックスの定義

a<-seq(0,3,0.01) 　　　　#スケールファクタを0.01刻みで0から3まで定義 B<-seq(0,3,0.01) 　　　　#数値積分の積分値を格納する箱 HT<-seq(0,3,0.01) 　　 #数値積分結果をNumericとして格納する箱 dHt<-seq(0,3,0.01) #現在値との差分を格納する箱

密度パラメータのパターンセット毎の計算結果の格納ボックスの定義

DM<-matrix(1:301,nrow=4,ncol=301)　　　#4パターンのパラメータセットで積分値を格納するための箱

フリードマン方程式の積分形の被積分関数の定義

f<-function(x) x/(R+Mx+Lx^4+(1-O)*x2)(1/2)　#被積分関数の定義

繰り返しの数値積分の実行

for(j in 1:4){
R<-BR[j] #輻射密度パラメータベクトルのj成分をRにセット M<-BM[j] #全物質密度パラメータベクトルのj成分をMにセット L<-BL[j] #宇宙定数密度パラメータベクトルのj成分をLにセット O<-BO[j] #曲率密度パラメータベクトルのj成分をOにセット for(i in 1:301){ B[i]=adaptIntegrate(f,0,a[i])[1] 　 #301個のaに対して数値積分しBに格納 HT<-as.numeric(B) #Numericに変換して格納 HT0<-adaptIntegrate(f,0,1)[1] #現在のH*t_0を計算 HT0<-as.numeric(HT0) #数値化 dHt[i]<-HT[i]-HT0 DM[j,]<-dHt #行列Mのj行目に値をセット } #現在値との差分を格納 }

プロット H(t-t0)を横軸、aとしたグラフを作成

枠作成

plot(0,0,type=“n”,xlim=c(min(DM),max(DM)), ylim=c(0,3), xlab=“H(t-t0)”, ylab=“a”)　#枠作成は一回のみ

グラフの色を定義

cols<-c(1,2,4,6)

2. グラフ作成

for(k in 1:4){
lines(DM[k, ],a,lty=1,col=cols[k]) #プロット }