GNU octave でグラフを書いてみよう

冬なので湿度をネタにグラフを書いてみる。
暖房をかけると湿度が下がるのはなぜなのか？
octave 4.4.1 を使います。

湿度とは？

空気中に存在できる最大の水蒸気量を 100% としたときの現在の水蒸気量の比。
空気中に存在することのできる最大の水蒸気量は温度によって決まり、温度が高いほど多くの水蒸気を含むことができる。(ある温度で $1m^3$の大気が含むことのできる最大の水蒸気の量を飽和水蒸気量 $[g/m^3]$ と表す)
※正しくは水蒸気圧と飽和水蒸気圧の比らしいですが、おおまかに水蒸気量と飽和水蒸気量の比で近似さるので気にしない。例えば大気が存在しないと、その分水蒸気を多く含むことができますが、大気圧を凡そ 1とした時の近似。

飽和水蒸気量のグラフを書いてみる

とりあえず線を引く

Wikipediaによると飽和水蒸気量と温度との関係式が以下のように書けるらしい。

飽和水蒸気量

a(T)=\frac{\displaystyle 217 \times 6.1078 \times 10 ^ \dfrac{7.5T}{T+237.3}}{T + 273.15}

$T$:温度 [℃]
$a(T)$: $T$[℃]の時の飽和水蒸気量 [g/m^3]

この飽和水蒸気量を GNU octave でプロットしてみると

飽和水蒸気量

T = linspace(-10,120);
a = ( 217*6.1078*10.^((7.5*T)./(T+237.3)) ) ./ (T + 273.15);
plot(T, a);

T = linspace(-10,120) ⇒ -10～120 を 100分割する (=-10～120℃)
a = ... ⇒ 飽和水蒸気量を求める
plot(T, a) ⇒ 横軸 T、縦軸 a のグラフを描く

こんな感じでグラフができます。

ちょいちょい計算式に不思議なピリオドが入っていますが、octave は基本的に行列演算のソフトなので、スカラーであることを明示するためのものです。(と思ってる)

線だけだと寂しいので、通過点に × を入れてみる

軸に単位がなかったり、するのでちょっと追加。

飽和水蒸気量のグラフ改

T = linspace(-10,120);
a  = 217*6.1078*10.^((7.5*T )./(T +237.3)) ./ (T +273.15);

Tx = -10:10:120;
ax = 217*6.1078*10.^((7.5*Tx)./(Tx+237.3)) ./ (Tx+273.15);
plot(T, a,'-', Tx, ax , 'xb');
grid on
xlabel("Celsius [^\\circC]");
ylabel("saturated [g/m^3]");
title("temperature vs saturated")

直線のみのグラフだったので Tx と、ax で 10℃ 毎に x をプロット。
xlabel, ylabel, title で軸 & グラフに見出しを追加。
Tx = -10:10:120 ⇒ -10 から 120まで、10 単位で計算する。
ax = ... ⇒ は変わらず。
plot(T, a,'-', Tx, ax , 'xb')
　↓
横軸に T、縦軸に a, 直線 `-' でプロット
横軸に Tx、縦軸に ax、`x' でプロットし色は青

飽和水蒸気量のグラフ改_2回プロット


T = linspace(-10,40);
a  = 217*6.1078*10.^((7.5*T )./(T +237.3)) ./ (T +273.15);
plot(T, a,'-r; line ;');
hold on

T = -10:10:40;
a  = 217*6.1078*10.^((7.5*T )./(T +237.3)) ./ (T +273.15);
plot(T, a,'^g; point ;');

grid on;
xlabel("T Celsius [^\\circC]", 'FontSize',14);
ylabel("a(T) saturated [g/m^3]", 'FontSize',14);
title("Celsius T vs saturated a(T)",
    'FontSize',16,
    'FontName', 'LimoneFitR',
    'FontWeight', 'bold',
    'color', 'b');

plot を 2回使ったこちらでも同じグラフが表示されます。
hold on と宣言することで以降の plot は消去せずにそのまま書きます。
普段の気温は -10～40℃まであれば十分かなということで範囲を変更。

plot コマンドのオプション

plot の第三引数 (フォーマット指定)

線の種類
‘-’ 直線 (default)	‘--’ 破線	‘:’ 点線	‘-.’ 一点鎖線

マークの種類
‘+’ プラス	‘o’ 円	‘*’ アスタリスク
‘.’ 点	‘x’ バツ	‘s’ 四角
‘d’ ダイヤ
‘^’ 上向き三角	‘v’ 下向き三角	‘>’ 右向き三角
‘<’ 左向き三角	‘p’ 五芒星	‘h’ 六芒星

色
‘k’ blacK	‘r’ Red	‘g’ Green	‘b’ Blue
‘m’ Magenta	‘c’ Cyan	‘w’ White

";key;" で、凡例に表示されます。

コマンドウィンドウでは日本語を入力してはいけません。入力すると

って感じになるのでやっちゃだめです。使うときはエディタで使いましょう。でも UTF で保存しちゃダメですよ。

温度と湿度から、水蒸気量を求める

計算式の確認

現在の温度 $T_{now}$ と、湿度 $H_{now}$ から現在の水蒸気量 $\rho_{now}$ は
飽和水蒸気量 $a(T_{now})$ x 湿度 $ H_{now}÷100$ で近似されるので

現在の水蒸気量

>> T = 20.0;
>> H = 45.0;
>> a = 217*6.1078*10.^((7.5*T )./(T +237.3)) ./ (T +273.15);
>> a * H / 100
ans =  7.7883

良く分からないけど、たぶんそんなもんだろう。

水蒸気量を求める式を使って湿度を計算しグラフにする

温度vs相対湿度_(水蒸気量一定)

1;
T_now = 20.0; # 現在の温度 (℃)
H_now = 40.0; # 現在の湿度 (%)

T_min = -40 # 計算する 最低温度
T_max = 120 # 計算する 最高温度
T_step = 5  # プロット間隔

T_disp_min = -10 # 表示する最低温度
T_disp_max = 40  # 表示する最高温度

a = 217*6.1078*10.^((7.5*T_now)./(T_now+237.3)) ./ (T_now+273.15); # 飽和水蒸気量を求める
R = a * H_now / 100 # 水蒸気量を求める

# グラフに線を引く
T = linspace(T_min,T_max);
H = 100 * R ./ (217*6.1078*10.^((7.5*T )./(T +237.3)) ./ (T +273.15));
plot(T, H,'-b; point ;');
hold on

# 5℃毎にグラフにマークをつける
T = T_min:T_step:T_max;
H = 100 * R ./ (217*6.1078*10.^((7.5*T )./(T +237.3)) ./ (T +273.15));
plot(T, H,'xb');

plot(T_now, H_now,'^r; given value ;'); #入力値を赤の三角でプロット

grid on;
xlabel("Celsius [^\\circC]", 'FontSize',14);
ylabel("R.H. [%]", 'FontSize',14);
title(["Celsius T vs Relative Humidity",
    sprintf("\n(%d^oC %d%% : %.2f g/m^3)", T_now, H_now, R)],
    'FontSize',16,
    'FontName', 'LimoneFitR',
    'FontWeight', 'bold',
    'color', 'b');

ylim([0 100]); # 湿度の表示範囲を 0～100 にする
xlim([T_disp_min T_disp_max]); # 温度の表示範囲を指定

このグラフから 20℃ 湿度40% の部屋が 25℃まで室温が上がると湿度が 30% に、15℃まで下がると 50% を超えてしまうことが分かります。

おまけ

水蒸気量一定としたときの温度-湿度グラフを描く

温度vs相対湿度グラフ

1;
T_min = -20; # 計算する 最低温度
T_max = 120; # 計算する 最高温度

T_disp_min = -10; # 表示する最低温度
T_disp_max =  40; # 表示する最高温度

T_and_H = [-30,-20,-10, 0,10,20,30,40,50,60,70,80,90,100;# 温度
            50, 50, 50,50,50,50,50,50,50,50,50,50,50,50] # 湿度

for i = T_and_H
  T_now = i(1)
  H_now = i(2)

  a = 217*6.1078*10.^((7.5*T_now)./(T_now+237.3)) ./ (T_now+273.15); # 飽和水蒸気量を求める
  R = a * H_now / 100 # 水蒸気量を求める

  # グラフに線を引く
  T = linspace(T_min,T_max,1000);
  H = 100 * R ./ (217*6.1078*10.^((7.5*T )./(T +237.3)) ./ (T +273.15));
  plot(T,H,'-');
  hold on
endfor

grid on;
xlabel("Celsius [^\\circC]", 'FontSize',14);
ylabel("R.H. [%]", 'FontSize',14);
title("Celsius T vs Relative Humidity",
    'FontSize',16,
    'FontName', 'LimoneFitR',
    'FontWeight', 'bold',
    'color', 'b');

ylim([0 100]); # 湿度の表示範囲を 0～100 にする
xlim([T_disp_min T_disp_max]); # 温度の表示範囲を指定

実行させるとこんな感じ。温度―湿度等水蒸気量曲線

間違えてオートスケールを押してしまうと

湿度 3万%とか意味わかりません。

T_and_H の行列から一発で水蒸気量 R が出るような計算を思いつけば速度出るんだろうけどねぇ。そして、その行列 R を使って H を求められたら...そこまで考える頭が(^^; のんびりやればできるのかなぁ。

paraarrayfun というのを使うと、octave のプロセスを別に立ち上げて並列に処理できるらしい。並列化 | Octaveやってみる！

GNU Octave (前編), (後編) 鷲沢嘉一