数量は原則として SI 単位であらわすが、「年」や「℃」などの慣用の単位をつかうこともある。 その場合は、秒やKなどのSI単位とのあいだの換算をおこなう。 単位は、「#」で始まるコメントとして書きこんだが、これは人が見るためのもので、演算では使っていない。 単位のうちで、べき乗をあらわす上つき添え字を、上つきでない通常の字で代用している。
例文のプログラムは、MS Windows上の R バージョン 4.0.2 で動作確認した。
(作業用ディレクトリを用意したうえで、 ブラウザで、つぎのリンクを、右クリックから「名まえをつけてリンク先を保存」あるいは「対象をファイルに保存」で保存してほしい。 そのさきの動作は、[R言語を使ってみる]のページを参照してほしい。)
まず(A)のプログラムにそって、どこで何をしているかを説明する。
# # 定数とモデルパラメータ # sigmab <- 5.67e-8 # ステファン・ボルツマン定数 [W m-2 K-4] tzero <- 273.15 # 水の凝固温度 [K] solconst <- 1.36e3 # 「太陽定数」 [W m-2] alphap <- 0.3 # 惑星アルベド [無次元量] sphcapac <- 4.2e3 # 比熱容量 [J kg-1 K-1] (水を想定) dens <- 1.0e3 # 密度 [kg m-3] (水を想定) depth <- 500 # 層の深さ [m] (海洋の風成循環のおよぶ層を想定) capac <- sphcapac * dens * depth # 面積あたり熱容量 [J m-2 K-1] yearsecond <- 60*60*24*365.25 # 1年の秒数 [秒] tinityear <- 0 # 初期時刻 [年] tinit <- tinityear * yearsecond # 初期時刻 [秒] deltatyear <- 1.0 # 時間刻み [年] deltat <- deltatyear * yearsecond # 時間刻み [秒] nstep <- 200 # ステップ数 [無次元整数]
まず、物理定数や、このモデルにとっては定数であるもの、 変えられるが1回の実行のあいだは定数とみなすモデルパラメータなどの値を、 Rの変数に代入する。 モデルの説明とつきあわせれば意味は推測できると思う。
数値の、10の指数を使った表現、たとえば「1.36×103」が、 「1.36e3」のように「e」を使って表現されている。 これは多くのプログラム言語に共通する表現で、 e は「指数」にあたる exponent からきているが、 「10の何乗」であって「eの何乗」ではないことに注意してほしい。
Rをふくめて多くのプログラム言語では物理量の単位をあつかうことができない。 単位については人が判断して、演算式は数値の関係で書く必要がある。 SI単位にそろえれば、たいていの場合、換算定数は不要になるが、 人が使いたい単位とSI単位がちがうばあいは、換算が必要になる。
ここではとくに、時間の単位として「年」を使いたいので、 「秒」との換算をふくめた。 同じ時間刻みを、年単位であらわした値を deltatyear、秒単位であらわした値を deltat としてある。
数値実験の設定として、時間刻みを 1年、ステップ数を 200 としてみた。 モデルの世界の 200年のシミュレーションをすることになる。 これは実行ごとに変更しうるものである。 層の深さを小さくしたときは、時間刻みも小さくしないといけないかもしれない。
数量の単位は、コメント(計算処理には関係なく、人が読むための部分。 Rでは各行の「#」からあとの部分)に、角かっこでくくって書いておいた。 ただし、10の指数を上つきにすべきところをそのまま書いている。
# 時間の配列を用意する tim <- tinit + seq(0, nstep) * deltat # 時間 [秒] timyear <- tinityear + seq(0, nstep) * deltatyear # 時間 [年] # 結果の配列を用意する (timと同じ長さの配列、内容は 0) temp <- tim * 0 # 温度 [K] fin <- tim * 0 # エネルギー収入 [W m-2] fout <- tim * 0 # エネルギー支出 [W m-2] # 初期値 tdegcinit <- 50 # kt = 1, 時刻 = 0 での 温度 [℃] temp[1] <- tdegcinit + tzero # kt = 1, 時刻 = 0 での 温度 [K]
時間発展型のシミュレーションのプログラムのつくりかたとしては、 各ステップごとに計算結果をファイルに書き出すやりかたと、 全ステップの計算結果をプログラム内の配列変数にとっておくやりかたがある。 シミュレーション対象を記述する変数の個数が多いときは、ファイルに書き出さないといけない。 ここでは変数の数が少ないので、配列変数にとっておくやりかたにした。
各ステップの時刻がはいる配列変数として、 ここでは、秒単位のもの(tim)と年単位のもの(timyear)の両方をつくって数値を入れている。 実用的には、これはどちらか1つにして、一方から他方の数値を計算して利用したほうがよいだろう。
各ステップの計算結果の変数値がはいる配列変数を、 寸法(要素の個数のこと)を tim と同じとして用意した。便宜上、値としてゼロを入れておく。
配列変数では、その要素を、「添え字」とよばれる整数(Rのばあいは角かっこ内に書かれる)で指定する。 プログラム言語によって、添え字が 0 から始まるもの、1 から始まるもの、始まりの数値を選べるものがある。 Rのばあいは、添え字は 1 から始まる。そこで、初期値を配列 temp の 添え字 1 の要素に入れている。
シミュレーションの初期値として、時刻 0 での温度の値を与える。 モデル内の温度 temp は、単位 K (ケルビン)の温度(いわゆる絶対温度)を使うことにした。 初期値を、℃の温度(セルシウス度) で数値を与えたいばあい、 その値を tdegcinit に与えておき、tzero (上の定数のところで値を与えてある)をたして 単位 K になおしてから、temp の最初の時刻のところに代入している。 もし初期値を 単位 K で与えたいのならば、tdegcinit をやめて、 直接 temp[1] への代入文にすればよい。
シミュレーションの本体というべき部分である。
# 時間ステップを deltat ずつ進めて計算
for(kt in 1:(nstep+1)){
fin[kt] <- solconst * (1 - alphap) / 4
fout[kt] <- sigmab * temp[kt]^4
if(kt <= nstep) {
temp[kt+1] <- temp[kt] + deltat * (fin[kt] - fout[kt]) / capac
}
}
時間を、0 から nstep × Δt まで進める。 kt という変数を 1 から 1 ずつ変えて、時刻 (kt-1)×Δt での数値を計算することにする。 まず「kt = 1」のばあい、つぎに「kt = 2」のばあい、と、 「kt = nstep+1」のばあいまで、順に計算していくことにする。 それを実現するのに「for」文をつかっている。 回数をかぞえながらくりかえす構文はプログラム言語ごとに表現がちがい、 ここで使った表現には R 特有のところがある。
fin[kt], fout[kt]に、時刻 (kt-1)×Δt でのエネルギー収入とエネルギー支出がはいる。 このとき temp[kt] にその時刻の温度がはいっていることを前提としている。
temp[kt+1]に、時刻 kt×Δt での温度がはいる。 (kt = nstep+1 になったら、そのつぎのステップの値のための配列は用意されていないので、 この計算は、kt <= nstep という条件がなりたつとき実行するようにした。
[注] [表計算ワークシートの例文]で、 ステップ番号を kt と書いたが、その kt は 0 から始まるものである。 このページで紹介する R の例文の kt は、表計算ワークシートの例文の kt+1 にあたる。
まず、実際の ex1.R とはちがうが、画面に作図するばあいをしめす。
# 作図 # 温度の時系列グラフ plot(timyear, temp, las=1, xlab="時刻 [年]", ylab="温度 [K]") # エネルギーの出入りの量の時系列グラフをPNG画像ファイルに書き出す plot(timyear, fout, las=1, xlab="時刻 [年]", ylab="エネルギーの出入り [W / m2]") lines(timyear, fin)
まず、横軸に時刻(単位は年)、縦軸に温度(単位はK)をとったグラフを 「plot」ひとつでかいている。 それから、別の図として、 横軸に時刻(単位は年)、縦軸にエネルギーの出入り(単位はW/m2)をとったグラフを、 まずエネルギー支出のほうを「plot」でかき、エネルギー収入を「lines」でかきたす、 という形でかいている。 (折れ線でなく記号列でかきたいときは lines のかわりに pointsをつかう。) 横軸、縦軸のラベルを xlab, ylab でかいている。 (「W/m2」の上つき添え字を上つきで表示する方法はあるのだが、 わたしが対応できておらず、ここでは、m とならべて 2 をかく形になってしまっている。)
実際の ex1.R では、グラフを画像ファイルに書き出す。 画像ファイルの形式は PNG にした。 作図の前に「png("ファイル名")」、あとに「def.off()」を入れる。
# 作図
# 温度の時系列グラフをPNG画像ファイルに書き出す
# # png(...) と dev.off() をはずせばグラフは画面に表示される
png("ex1_temp.png")
plot(timyear, temp, las=1, xlab="時刻 [年]", ylab="温度 [K]")
dev.off()
# エネルギーの出入りの量の時系列グラフをPNG画像ファイルに書き出す
png("ex1_fluxes.png")
plot(timyear, fout, las=1, xlab="時刻 [年]", ylab="エネルギーの出入り [W / m2]")
lines(timyear, fin)
dev.off()
結果の図はつぎのようになる。
計算結果を数値としてファイルに書きだしておき、別のプログラムで使いたいことがある。 ファイルでの数値の表現はいろいろあるが、 ここでは、テキストファイルに数字の表で書きだしておくことにする。
時間ステップ計算のあいだは数値を配列変数にとっておき、 最後にまとめて書き出すことにする。
Rの「データフレーム」という機能を使い、 配列変数の内容をデータフレームにまとめてから、 Rの write.table という関数を使って書き出すことにする。
# データフレームにまとめる df <- data.frame(timyear, temp, fin, fout) # データフレームをファイルに書き出す write.table(df, "ex1output.txt", quote=F, header=T)
data.frame をこのように使うと、配列変数名がデータフレーム df の列名として使われる。
write.table のオプションの quote=F は、書き出す各項目を引用符でくくらないことを指定している。 row.names=F は、行に名まえをつけないことを指定している (これがないと行番号の列ができる)。
結果のファイル ex1output.txt [1は数字の1] には、 1行めに列名、2行めからは変数値が、スペース(空白)くぎりで書かれている。
これは、Excel でも、スペースくぎりのテキストファイルであることを指定すれば、 表に読みこむことができる。 (Excel 2013のばあい、「テキストファイルウィザード」で、 「カンマやタブなどの区切り文字によってフィールドに区切られたデータ」を選択し、 「次へ」進んで、「区切り文字」の「スペース」にチェックを入れてから次に進む。)
(B)で書き出したファイルを読んで、(A)と同様に作図する。
実際の ex1readf.R は グラフをPNG画像ファイルに書き出すようにしてあるが、 ここではPNG画像関係の処理を省略してしめす。これを実行するとグラフを画面に表示する。
# データをファイルから読みこむ
df <- read.table("ex1output.txt", header=T)
# 作図
# 温度の時系列グラフ
plot(df$timyear, df$temp, las=1, xlab="時刻 [年]", ylab="温度 [K]")
# エネルギーの出入りの量の時系列グラフ
plot(df$timyear, df$fout, las=1, xlab="時刻 [年]", ylab="エネルギーの出入り [W /m2]")
lines(df$timyear, df$fin)
write.table で書き出されたデータフレームは、read.table で読みこむことができる。 作図をする手順は (A) とほとんど同じだが、 縦・横どちらも、数値を与える変数が、(A)では配列変数だったが、 (C)ではデータフレームの列なので、「df$」が前についている。
この例題で、数値を変えることを想定しているモデルパラメータは、 「層の深さ」depth の値である。 これを変えると、「時間刻み(年)」 deltatyear の値も変える必要があるかもしれない。 層の深さを小さくすると面積あたりの熱容量も小さくなり、 温度の変化が速くなるので、時間刻みを短くする必要があるだろう。 シミュレーション対象となる時間は、時間刻みとステップ数 nstep の積である。 「初期時刻(年)」tinityear は、この例題では 0 のまま変えないのがよいだろう。 (別の例で、たとえば西暦2020年からのシミュレーションをすることを考えて、変えられるようにしておいた。)
また、モデルパラメータの「太陽定数」solconst と「惑星アルベド」alphap も、 例題としては一定値を与えたが、変更しうる量である。
時間ステップごとの計算の内容は、エネルギー保存の法則にもとづいて、 単位面積あたりのエネルギーのたまりの量の変化をもとめているのだが、 便宜上、「エネルギーのたまりの量」の変数をもちださずに、 温度の変化をもとめる形に書いた。 このモデルでは、単位面積あたりのエネルギーのたまりの量は、 温度 temp[kt] に、単位面積あたりの熱容量 capac をかけたものであり、 capac は時間とともに変わらないと仮定したから、このような単純化ができた。 (エネルギーのたまりの量がほしければ、temp[kt] に capac をかければよい。) 熱容量が時間とともにかわりうるようなモデルでは、 このような単純化ができず、 温度とエネルギーのたまりの量のそれぞれの変数をもつ必要があるだろう。
この例文プログラムで、複数の計算結果を比較するためには、 まずは、png("...") の ... のところのファイル名を毎回書きかえて実行し、 画像をならべて見るのがよい。 複数の時系列をひとつのグラフにまとめて表示したければ、 (C)のプログラムを改造して、plot を呼び出してひとつの時系列をかいたあとに、 lines または points を呼び出してほかの時系列をかけばよい。 [指数関数型の成長]のページの Rプログラム例を参考にしてほしい。
まず、モデルパラメータを変更せずに、温度の初期値 tdegcinit に入れる数値を なんとおりか変えて、結果を比較してみよう。
つぎに、層の深さ depth の数値を、なんとおりか変えて、 初期値のほうは変えないで計算してみよう。 ただし、層の深さを小さくすると、時間刻み (「年」単位) の値 deltatyear も 小さくする必要があるかもしれない。