gnuplot によるグラフ作成 - 弘前大学 Home Sweet Home
葛西 真寿 弘前大学大学院理工学研究科
はじめに¶
この Notebook では，gnuplot_kernel (A Jupyter/IPython kernel for Gnuplot) を使って，関数およびデータのプロットプログラムである gnuplot の使用例を紹介する。
セクションの構成は，テキスト「wxMaxima による数式処理とグラフ作成」の「第2章 wxMaxima によるグラフ作成」や，「Sympy による数式処理とグラフ作成」の 「2 Sympy や Matplotlib によるグラフ作成」に準じているので，Jupyter Notebook の活用例としてだけでなく，wxMaxima や SymPy，Matplotlib などのグラフ作成方法を比較するのに役立つかもしれない。
演算子・定数・関数¶
gnuplot では以下のような（FORTRAN と同じような... と言っても，もはや通じない）演算記号を使うことができる。
足し算 +，引き算 -，かけ算 *，割り算 /，べき乗 **
print コマンドを使えば，以下のように簡易電卓として使うこともできる。
print 12**2
数の表現方法も（FORTRAN と同様であり... と言っても通じない），整数と実数がある。
print 1/3
print 1./3
非常に大きい数や，とても小さい数を表すには...
（例: $M = 5.972 \times 10^{24}, q = 1.0.6 \times 10^{-10}$）
M = 5.972e24
q = 1.6e-19
gnuplot には定数として円周率 pi がはじめから組み込まれている。
print pi
また、組み込み関数としてよく使われるものをあげると，
sqrt(x)= $\sqrt{x}$ |
abs(x)= ｜$x$｜ |
|
exp(x)= $e^x$ |
log(x)= $\log_e x$ |
log10(x) = $\log_{10} x$ |
sin(x) |
cos(x) |
tan(x) |
sinh(x) |
cosh(x) |
tanh(x) |
asin(x)=$\arcsin x$ |
acos(x)=$\arccos x$ |
atan(x)=$\arctan x$ |
gnuplot には，このほか，ユーザが自分で定数や関数を定義できる機能もある。定数については，上で定義した M や q の例がすでにあるので...
print M
関数を定義する例をあげる。組み込み関数 sin(x) に入れる x の値は radian 単位である。
print sin(30)
$\displaystyle \sin \frac{\pi}{6} = ...$
print sin(pi/6)
以下のような関数 degsin(x) を定義してやると...
degsin(x) = sin(x*pi/180)
x = 30
print degsin(x)
ユーザが定義した定数や関数は，gnuplot を終了すると忘れ去られてしまう。次回もまた使いたいときは，
save "defs.gpt"
などとして現在の定義をファイルに保存し，次に gnuplot を使うときに，
load "degs.gpt"
として読み込んでやればよい。
なお，このページのように Jupyter notebook を使えばユーザが定義した定数や関数，そして描画したグラフもすべて notebook に保存されるので，とても便利。
練習問題: 重力加速度の大きさ¶
gnuplot を電卓として用い，重力加速度の大きさ $g$ を求める。 地球の質量を$M$、半径を$R$、万有引力定数を$G$とすると、
$$\displaystyle \ g = \frac{GM}{R^2} $$である。
$M = 5.972 \times 10^{24} \mbox{kg}$, $R = 6.378 \times 10^6 \mbox{m}$, $G = 6.672 \times 10^{-11} \mbox{N} \mbox{ m}^2/\mbox{kg}^2$ のとき、$g$はいくらか？
gnuplot によるグラフ作成¶
まず，ターミナルを開いて，シェルのプロンプトで gnuplot と入力し，Enter キーを押すと以下のようになる。
$ gnuplot
G N U P L O T
Version 5.4 patchlevel 0 last modified 2020-07-13
Copyright (C) 1986-1993, 1998, 2004, 2007-2020
Thomas Williams, Colin Kelley and many others
gnuplot home: http://www.gnuplot.info
faq, bugs, etc: type "help FAQ"
immediate help: type "help" (plot window: hit 'h')
Terminal type is now 'qt'
gnuplot>
以下では，Jupyter kernel for gnuplot を使っているので gnuplot のプロンプト gnuplot> は表示されない。
2次元グラフ plot の例¶
2次元グラフ $y = f(x)$ を $a \le x \le b$ の範囲で描く表記例は，plot [a:b] f(x) です。
例として，$y = \sin x$ のグラフを $-2\pi \le x \le 2 \pi$ の範囲で描きます。基本的な定数の一つである円周率 $\pi$ は gnuplot では pi と書きます。
reset
# 念の為`set` コマンドで定義できるグラフに関する
# 全てのオプションをデフォルトの値に戻します。
plot [-2*pi:2*pi] sin(x)
いくつかオプションを指定して描く例です。
以下では，座表軸を表示し，線の色（lc）を設定し，凡例（title）を表示させ，x軸とy軸にラベルを表示させて $y = \cos x$ を描いています。
set zeroaxis
set xlabel 'x'
set ylabel 'y'
plot [-2*pi:2*pi] cos(x) title '余弦関数' lc 2
$x$ の範囲と $y$ の表示範囲を両方指定してプロットする例。
plot [-pi:pi] [-2:2] tan(x)
3次元グラフ splot の例¶
3次元グラフ $z = g(x, y)$ を描く一般的な表記例は splot [a:b] [c:d] g(x, y) です。
例として，$z = x \sin y$ のグラフを $-3 \le x 3, -4 \le y \le 4$ の範囲で描きます。
# 上記で設定したオプションをデフォルトに戻す
reset
splot [-3:3] [-4:4] x*sin(y) notitle
デフォルトでは，splot は透明な曲面を描きます。曲面の重なった部分を非表示にしたい場合は，set hidden3d です。（解除するには，unset hidden3d です。）
曲面を滑らかにするには，isosamples の数値を増やします。
set hidden3d
show isosamples
set isosamples 50, 50
replot
数値データのグラフ作成例¶
配列の定義と利用¶
gnuplot では，数値データの配列をつくってグラフを描くことができます。
以下のように array で1次元配列を定義し，値を入れます。配列の番号は 1 からはじまります。
array X[6] = [0, 1, 2, 3, 4, 5]
array Y[6] = [0, 1**2, 2**2, 3**2, 4**2, 5**2]
do for ループ¶
以下のように，do for を使ったループで配列に値を入れることもできます。
array X[6]
array Y[6]
do for [i=1:6]{
X[i]=i-1
Y[i]=(i-1)**2
}
do for [i=1:6]{print X[i], Y[i]}
配列 X[i] を $x$，配列 Y[i] を $y$ の値として $(x, y)$ の点をプロットする例です。
reset
plot Y using (X[$1]):2
いくつかオプションを指定してプロットする例。ここでは，w lp (with linespoints の省略形) で各点を線でつないで描きます。
plot Y using (X[$1]):2 w lp pt 7
丸点の大きさが，svg にしたせいか大きめです。pointsize (ps) を 0.6 にしてみます。
set pointsize 0.6
replot
複数のグラフを重ねて表示¶
gnuplot で複数のグラフを重ねて表示する例を示します。
$y = x^2 - 1$ と $y = 4 x - 5$ の2つのグラフを重ねて描きます。$x$ の範囲は $-5 \le x \le 5$ です。
reset
plot [-5:5] x**2 - 1, 4*x - 5
いくつかオプションを設定して描く例です。
以下の例では，set zeroaxis で x 軸と y 軸を表示させ，
set grid; set xtics 1; set ytics 1 で x, y それぞれ 1 ごとに格子線を表示させ，
横軸 (x 軸)，縦軸 (y 軸）の表示範囲をそれぞれ
$-1 \le x \le 4$，
$-2 \le y \le 6$ に制限して，$x^2 - 1$ グラフの線を lw (linewidth) 2 で少し太めにして表示します。
$x=2, y=3$ あたりで両線が接している様子が見えてきます。
set zeroaxis
set grid
set xtics 1
set ytics 1
set key left top
plot [-1:4] [-2:6] x**2 - 1 title "x^2 - 1" lw 2, 4*x - 5 title "4 x - 5"
数値データを読み込む¶
gnuplot あらかじめ作成された数値データファイルを読み込んでグラフを描くこともできます。
以下のような内容のファイル test.dat がカレント・ディレクトリにあるとします。
# x y
0 0
1 1
2 4
3 9
4 16
5 25
このデータをプロットするには，以下のようにファイル名を指定して plot します。
pt (pointtype) 7 で点種を，lc (linecolor) 7 で色を設定しています。
reset
set pointsize 0.6
plot 'test.dat' pt 7 lc 7
数値データと理論曲線を重ねて表示¶
前節の数値データファイル test.dat と理論曲線 $y = x^2$ の2つのグラフを重ねて表示してみます。
また，グラフの線と重ならないように，set key left top で凡例の表示位置を左上にします。
set key left top
plot 'test.dat' pt 7 lc 7, x**2 title "x^2"
媒介変数表示の2次元グラフ¶
半径1の円の方程式は $x^2 + y^2 = 1$ です。
この円を描くには，$y = \pm\sqrt{1-x^2}$ として $y = f(x)$ の形にするよりも，以下のような媒介変数表示にしたほうが簡単でしょう。
$$ x = \cos t, \quad y = \sin t, \quad(0 \le x \le 2\pi)$$このような媒介変数表示の2次元グラフを gnuplot で描くには以下のようにします。
最初に set parametric とします。
reset
set parametric
plot [0:2*pi] cos(t), sin(t) title "半径 1 の円"
これは媒介変数表示の円のグラフですが，横につぶれて楕円のように見えます。せっかくですから，グラフの縦横の比を$1:1$にして円らしく見えるようにします。
set xrange [-1.2:1.2]
set yrange [-1.2:1.2]
set zeroaxis
set key samplen 1
set size square
replot
非媒介変数モードに戻るには，unset parametric と入力します。
なお，set parametric としたときの媒介変数は，デフォルトでは t であるが，どうしてもそれ以外の文字，例えば phi を使いたのであれば以下のようにする。
set xrange [-1.2:1.2]
set yrange [-1.2:1.2]
set parametric
set size square
plot [phi=0:2*pi] cos(phi), sin(phi) notitle
（練習） 楕円のグラフ（楕円中心を原点として）
同様にして，楕円のグラフを描くこともできます。原点を中心とし，長半径 $a$，短半径 $b$ の楕円の式は
$$ \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1 $$です。媒介変数表示では，
$$ x = a \cos t, \quad y = b \sin t \quad (0 \le t \le 2\pi)$$と書けます。$a, b$ に適当な値を入れて楕円のグラフを描きなさい。
海王星と冥王星の軌道¶
太陽からの万有引力を受けて運動する惑星（惑星だけでなく，準惑星や小天体も含みます）は，太陽を焦点とした楕円軌道を描きます。焦点を原点とし，長半径 $a$，離心率 $e$ の楕円の方程式は，極座標 $r, \phi$ を使って以下のように表すことができます。
$$ r = \frac{a(1-e^2)}{1 + e\cos \phi}$$さて，かつては第9惑星，現在では準惑星の一つである冥王星も楕円軌道を描きます。冥王星の軌道長半径 $a_P = 39.767 \,\mbox{au}$，離心率 $e_P = 0.254$ を使って冥王星の軌道を描きます。
まず，極座標表示の楕円の式を関数として定義します。
r(a, e, phi) = a*(1-e**2)/(1 + e*cos(phi))
極座標表示の2次元グラフ¶
極座標表示で描く場合は，set polar とします。（解除は unset polar です。）
reset
set zeroaxis
set xrange [-50:50]
set yrange [-50:50]
set polar
set size square
aP = 39.767
eP = 0.254
plot [phi=0:2*pi] r(aP, eP, phi) title '冥王星'
では，この冥王星の軌道と，その内側をまわる海王星の軌道を重ねて描いてみましょう。
海王星の軌道長半径は $a_N=30.1104\,\mbox{au}$，離心率は $e_N=0.0090$ と小さいので簡単のために $e_N = 0$ として扱います。実際の2つの天体の軌道は同じ平面上にありませんが，太陽からの距離のみをみるために，ここでは同一平面上に描きます。
以下の例では，凡例の位置を set key outside で図の外に設定しています。
set key outside
set grid
aN = 30.1104
eN = 0
replot r(aN, eN, phi) title '海王星'
月別気温のグラフと関数フィット¶
以下のような内容のファイル aomori.dat がカレント・ディレクトリにあるとします。
# 青森の月別気温 「理科年表 2019」(丸善)より。
# 月別 日最高 日最低
# 平年 気温の 気温の
# 月 気温 月 別 平 年 値
1 -1.2 1.6 -3.9
2 -0.7 2.3 -3.7
3 2.4 6.3 1.3
4 8.3 13.5 3.7
5 13.3 18.4 8.9
6 17.2 21.7 13.5
7 21.1 25.4 18.0
8 23.3 27.7 19.8
9 19.3 24.0 15.1
10 13.1 18.0 8.6
11 6.8 10.9 3.0
12 1.5 4.6 -1.4
1列目の「月」を x，2列目の「月別平年気温」を y としてプロットする例。
reset
set pointsize 0.6
set grid
plot [1:12] "aomori.dat" using 1:2 pt 7 title '平年気温'
グラフをみると，月別平年気温は12ヶ月周期の正弦関数または余弦関数のように見えます。では，以下のように関数フィットをしてみましょう。
以下のような関数でフィットしてみます。
$$f(x) = a - b \cos\left(\frac{2\pi (x-c)}{12}\right)$$a = 1
b = 1
c = 1
f(x) = a - b*cos( 2*pi*(x-c)/12 )
fit f(x) "aomori.dat" using 1:2 via a, b, c
print a, b, c
最小二乗法によって求めた $a, b, c$ の値をもった $f(x)$ のグラフを重ねて描きます。
set key Right at 5, 24
replot f(x)
日最高気温，日最低気温の月別平年値も合わせてプロットしてみてください。
stats コマンドによる統計情報¶
ちなみに，関数フィットした際の定数 $a$ は月別平年気温の年平均値に相当します。平均値を求めるだけなら，最小二乗法などやらなくても以下のようにして求めることができます。
* COLUMN:
Mean: 10.3667stats "aomori.dat" using 2
gnuplot を使ってもう少し...¶
gnuplot を使って，さらにいくつか問題に取り組んでみます。
波動¶
x 方向に進む波¶
$\sin(x - v t)$ は$x$の正の方向に速度$v$で進む正弦型の波動をあらわす。 波が進む様子を描いてみよう。
簡単のために，$v = 1$ とする。時刻 $t=0$ では，この波は $\sin x$。少し時間がたった後，たとえば $t = 0.5$ では，$\sin (x - 0.5)$，さらに，$\sin (x - 1.0)$，，$\sin (x - 1.5)$ など...
reset
set key outside
set key samplen 1
set samples 500
set zeroaxis
set xrange [-10:10]
set yrange [-1.1:1.1]
plot sin(x) lw 2 lt 7 title "sin(x) "
set key outside で凡例をグラフ外に表示させると，グラフの横サイズが狭くなります。Jupyter Notebook 内でグラフサイズを変更する例です。
%gnuplot inline svg size 640, 384 font "NotoSansMonoCJKjp-Regular.otf,14"
replot
# gnuplot 内では以下のように。
# gnuplot> set term qt size 640, 384
plot sin(x) title "sin(x) " dt 3, \
sin(x-0.5) lw 2 lc 7
plot sin(x) title "sin(x) " dt 3, \
sin(x-0.5) dt 3, \
sin(x-1.0) lw 2 lc 7
plot sin(x) title "sin(x) " dt 3, \
sin(x-0.5) dt 3, \
sin(x-1.0) dt 3, \
sin(x-1.5) lw 2 lc 7
while ループと GIF アニメ¶
Jupyter Notebook 環境では，簡単にアニメーションを作ることはできないが，以下のようにすれば見れそう。
set xrange [-10:10]
set yrange [-1.1:1.1]
i = 0
imax = 20
while(i < imax){
plot sin(x - 0.628*i) notitle
pause 0.5
i = i + 1
}参考：以下のようにすれば，GIF アニメを作れる。
set terminal gif animate delay 10 optimize size 512,384
set output "hadou.gif"
set samples 512
set zeroaxis
set xlabel ' >>> +x '
set xrange [-10:10]
set yrange [-1.1:1.1]
i = 0
imax = 20
while(i < imax){
plot sin(x-0.314*i) notitle
i = i + 1
}
unset output
うなり¶
振動数が少しだけ異なる2つの波を重ね合わせると「うなり」が起こることを示せ。
%gnuplot inline svg size 512, 384 font "NotoSansMonoCJKjp-Regular.otf,14"
reset
set key below
set samples 500
set zeroaxis
set yrange [-2.1:2.1]
plot [-100:100] sin(x + 0.05*x) + sin(x - 0.05*x)
$\sin(a+b) + \sin(a-b) = 2 \sin a \cos b$ ですから，うなりとなる波形を描いてみます。
replot 2*abs(cos(0.05*x)) lw 2 lc 7 notitle
つまり，$\omega_1 = a+b, \omega_2 = a-b$ の波を重ねると， 振動数の差の半分 $\displaystyle b = \frac{\omega_1 - \omega_2}{2}$ の コサイン波 $2 \cos b t$ が生じるように思えるが...
正確には，その絶対値 $2 | \cos b t|$ なので，振動数は実質その2倍，つまり うなりの振動数は $2 b = \omega_1 - \omega_2$ となり，
「振動数のわずかに異なる2つの波を重ねると，それらの振動数の差の振動数をもつうなりが発生する」
ということになる，という理解でどうでしょう。
定常波¶
波長・振幅が等しい正弦波が直線上を逆向きに同じ速さで進んでいると，重ね合わせによってどのような波が生じるか。
$$ \sin (x- v t) + \sin(x + v t)$$reset
set zeroaxis
set samples 500
set yrange [-2.1:2.1]
plot sin(x - 0) + sin(x + 0) notitle lw 2 lt 7
plot sin(x - 0) + sin(x + 0) notitle dt 3, \
sin(x - 0.5) + sin(x + 0.5) notitle lw 2 lt 7
plot sin(x - 0) + sin(x + 0) notitle dt 3, \
sin(x - 0.5) + sin(x + 0.5) notitle dt 3, \
sin(x - 1.0) + sin(x + 1.0) notitle lw 2 lt 7
plot sin(x - 0) + sin(x + 0) notitle dt 3, \
sin(x - 0.5) + sin(x + 0.5) notitle dt 3, \
sin(x - 1.0) + sin(x + 1.0) notitle dt 3, \
sin(x - 1.5) + sin(x + 1.5) notitle lw 2 lt 7
以下のような GIF アニメを作ってみてください。
一様重力場中の質点の運動¶
水平な地上で，水平となす角 $\theta$ で斜め上方に初速度の大きさ $v_0$ でボールを投げ出す。 $v_0$ を一定としたとき，水平到達距離を最大にする $\theta$ は？ またそのときの水平到達距離は？
無次元化¶
$t = 0$ で $x(0) = 0, y(0) = 0$ とすると，
$$ x = v_0 t \cos \theta , \quad y = v_0 t \sin\theta - \frac{1}{2} g t^2$$となるが，$g$ や $v_0$ に実際の次元のついた量を入れるのではなく，まずはこの式を，この系に特徴的な量で無次元化することを考える。
今，$\displaystyle \theta=\frac{\pi}{2}$ として垂直に投げ上げると，鉛直方向の $y$ は
$$ y = v_0 t - \frac{1}{2} g t^2 = \frac{v_0^2}{2 g} - \frac{1}{2} g \left( t - \frac{v_0}{g}\right)^2 $$となることから，特徴的な量として，$y$ が最大値をとる時刻 $\displaystyle t_m \equiv \frac{v_0}{g}$ と，その時の $y$ の最大値 $\displaystyle y_m \equiv \frac{v_0^2}{2g}$ を使って，以下のような無次元量をつくる。
$$ \bar{x} \equiv \frac{x}{y_m}, \quad \bar{y} \equiv \frac{y}{y_m}, \quad \bar{t} \equiv \frac{t}{t_m}$$これらの無次元量を使って軌道をあらわすと，
$$ \bar{x} = 2 \bar{t} \cos\theta, \quad \bar{y} = 2 \bar{t} \sin\theta - \bar{t}^2 $$以後は簡単のために $\bar{\ }$（バー）を省略する。
$\theta$ の値を変えて，$t$ を媒介変数としたグラフをいくつか描いてみる。
reset
x(theta,t) = 2*t*cos(theta)
y(theta,t) = 2*t*sin(theta) - t**2
set parametric
set zeroaxis
plot [0:2] x(pi/4, t), y(pi/4, t) title "θ=π/4"
媒介変数 $t$ の範囲を適当に [0:2] とした図が上図であるが，$y$ の値がマイナスの領域まで描いてしまって，このボールが地面にめり込んでしまっている！！
以後は $\bar{\ }$ をとる約束だったので， $\displaystyle y = 2 t \sin\theta - t^2 = 0 $ から $t$ の範囲は $\theta$ によって異なり，$0 \le t \le 2 \sin\theta$ であることがわかるので， 以下のようにして，地面にめりこまないように，$y$ がふたたび $0$ になる時刻までのプロットにする。
th1=pi/4
plot [t=0:2] x(th1, t*sin(th1)), y(th1, t*sin(th1)) title 'θ=π/4'
th2=pi/6
replot x(th2, t*sin(th2)), y(th2, t*sin(th2)) title 'θ=π/6'
th3=pi/3
replot x(th3, t*sin(th3)), y(th3, t*sin(th3)) title 'θ=π/3'
さて，肝心の水平到達距離を最大にする $\theta$ だが，地上から投げ上げたボールが再び地上に戻ってくるまでの時間は，（無次元化された時間で） $t = 2\sin\theta$ ということが，上述の議論でわかっているので，これを $x = 2 t \cos\theta$ に代入して，
$$ x = 4 \sin\theta \cos\theta = 2 \sin 2\theta$$これから，ただちに $x$ は $\theta=\pi/4$ のとき最大値 $2$ をとることがわかる。
あたらめて無次元化するまえの $x$ についてみると，最大到達距離は
$$x_{\max} = \bar{x}_{\max} \, y_m = \frac{v_0^2}{g}$$ケプラーの第3法則¶
以下のような内容のファイル planets.dat がカレント・ディレクトリにあるとします。
# 惑星と一部の準惑星「理科年表2019」(丸善)より。
#
# 軌道長半径 離心率 平均周期
# a(天文単位) e T (年)
水星 0.3841 0.2056 0.24085
金星 0.7233 0.0068 0.61520
地球 1.0000 0.0167 1.00002
火星 1.5237 0.0934 1.88085
木星 5.2026 0.0485 11.8620
土星 9.5549 0.0554 29.4572
天王星 19.2184 0.0463 84.0205
海王星 30.1104 0.0090 164.7701
# Ceres 2.767 0.076 4.60
冥王星 39.767 0.254 248
# Eris 67.668 0.156 5612列目の軌道長半径を x に，4列目の平均周期を y にしてプロットしてみます。
reset
set grid
set xlabel "軌道長半径 a"
set ylabel "公転周期 T"
plot "planets.dat" using 2:4 pt 7 lc 7 ps 0.7 notitle
対数グラフ¶
直線で近似できるような単純な比例関係ではなさそうです。
以下のように，set logscale xy として両対数グラフにしてみます。
set logscale xy
replot
両対数グラフにすると，データは直線に乗っているようにみえます。ということは， $\displaystyle T = K a^p$ のような関係になっていることが予想されます。
また，$a = 1$ のとき $T=1$ であることから，$K$ の値はほぼ $1$。
さらに，$a=10$ のあたりで $T$ の値は $10$ と $100$ の間であることから，$1 < p < 2$ であることが予想されます。
f(x) = K*x**p
fit f(x) "planets.dat" using 2:4 via K, p
print K, p
replot f(x)
$\displaystyle T = 1.04 \ a^{1.49}$ のような関係が得られました。
実際のケプラーの第3法則は，$T \propto a^{1.5}$ すなわち，$\displaystyle T^2 \propto a^3$ です。
データファイルから配列へ代入¶
# データの行数の取得。stats コマンドは数値データしか扱えない。
stats "planets.dat" using 2 nooutput # 1 桁目は数値ではないので，2桁目で。
N = STATS_records # データの行数
# 配列の宣言
array a[N] # 軌道長半径
array T[N] # 周期
# データファイルから配列へ代入
stats "planets.dat" using (a[$0+1] = $2, 0) nooutput
stats "planets.dat" using (T[$0+1] = $4, 0) nooutput
# 惑星名。stats は文字データは読めないらしいので，手動で代入。
array P[9] = ["水星","金星","地球","火星","木星","土星",\
"天王星","海王星","冥王星"]
# 配列の内容の表示
do for [i=1:N]{
print P[i], " ", a[i], " ", T[i]
}
簡単のために，惑星（および冥王星）を離心率 $e=0$ の円軌道として，1年の間にどれだけ公転運動するかを描いてみます。
半径 $a_i$，公転周期 $T_i$ の天体の $t=0$ での初期位置を $x$ 軸上として，時刻 $t$ における位置は
$$x = a_i \cos \frac{2\pi t}{T_i}, \quad y = a_i \sin \frac{2\pi t}{T_i}$$ですから...
reset
set parametric
set zeroaxis
set key left top
set key samplen 1
plot [t=0:T[3]] a[3]*cos(2*pi*t/T[3]), a[3]*sin(2*pi*t/T[3]) title P[3]
replot a[4]*cos(2*pi*t/T[4]), a[4]*sin(2*pi*t/T[4]) title P[4]
replot a[5]*cos(2*pi*t/T[5]), a[5]*sin(2*pi*t/T[5]) title P[5]
replot a[6]*cos(2*pi*t/T[6]), a[6]*sin(2*pi*t/T[6]) title P[6]
地球の公転軌道が円にみえるように xrange, yrange, ratio を調整して描きます。
set xrange [-2:10]
set yrange [-2:4]
set size ratio 0.5
replot
- 残りの天体の軌道も描いてみましょう。
- GIF アニメにしてみましょう。
- 一番内側の水星がちょうど一周する時間に，残りの天体はどこまで運動するかがわかるように描いてみましょう。
人口問題¶
マルサスの人口モデル¶
人口変化のしかたを明らかにし，将来の変化を予測するモデルを定式化するこ とは実社会において非常に重要な問題である。ここでは，まずマルサスの 「人口論」に書かれたアイデアを紹介する。
時刻 $t$ におけるある国の総人口を $N(t)$ とする。短い時間間隔 $dt$ における 出生数と死亡数は，ともにその時点での総人口と時間間隔に比例するから，
$$\mbox{出生数} = \alpha N dt, \quad \mbox{死亡数} = \beta N dt $$したがって，時間間隔$dt$での人口の増加$dN$は
$$dN = \alpha N dt - \beta N dt = \gamma N dt, \quad\mbox{ただし}\ \gamma \equiv \alpha - \beta $$これは微分方程式
$$\frac{dN}{dt} = \gamma N $$であり，次のように解ける。
$$N(t) = N_0 \,e^{\gamma (t - t_0)} $$ここで，$N_0 = N(t_0)$ は $t_0$ 年の人口である。また，$\gamma$ は，$t_0$ 年と $t_1$ 年の人口を使って以下のように求めることができる。
$$ \gamma = \frac{1}{t_1 - t_0} \ln\left(\frac{N(t_1)}{N_0} \right)$$以下のような内容のファイル usa.dat がカレント・ディレクトリにあるとします。
# usa.dat
#
# 年 米国の人口(単位: 百万人)
1790 3.9
1800 5.3
1810 7.2
1820 9.6
1830 12.9
1840 17.1
1850 23.2
1860 31.4
1870 38.6
1880 50.2
1890 62.9
1900 76.0
1910 92.0
1920 106.5
1930 123.2この数値データをプロットします。
# 前の設定が残っていると困るので，元に戻しておきます。
reset
set grid
set key left top
set key samplen 1
set xlabel '年'
set ylabel '人口（単位：百万人）'
plot "usa.dat" title "米国の人口" pt 7 ps 0.6
人口データから，$t_0, N_0, \gamma$ の値を求め，マルサスモデル $N(t) = N_0 \,e^{\gamma (t - t_0)}$ を重ねて描きます。
t0 = 1790
N0 = 3.9
gamma = 1./(1800-1790) * log(5.3/3.9)
Nm(t) = N0*exp(gamma*(t-t0))
replot Nm(x) title "マルサスモデル"
上の結果から，何がわかりますか？
ヴェアフルストによる修正モデル¶
マルサス・モデルは，$\gamma > 0$ の場合には人口の際限ない増加を予測する。 しかし，実際には食糧資源の供給不足，人口の過密，その他の環境的要因によ り，このような無制限な増加は続かない。
ヴェアフルストは，人口過密の要因を考慮にいれて，次のような修正を提案し た。人口は継続し続ける限り増加するが，上限 (それを$N_{\max}$としよう) があるとする。そして，人口変化は，次の各々に比例すると仮定する:
- 現在の人口 $N$
- 未使用の人口資源に対する割合 $(1 - N/N_{\rm max})$
したがって微分方程式は
$$ \frac{dN}{dt} = \gamma N \left(1 - \frac{N}{N_{\rm max}}\right) $$となり，（変数分離法によってこの式は解けて）解は
$$ N(t) = N_0 \frac{N_{\rm max}}{N_0 + (N_{\rm max}-N_0)\,e^{-\gamma (t-t_0)}} $$パラメータ $N_{\max}$ をどのように選ぶと，米国の人口変化をよく再現するか，調べてみます。
t0 = 1790
N0 = 3.9
gamma = 1./(1800-1790) * log(5.3/3.9)
Nv(t) = N0*Nmax/(N0 + (Nmax - N0)*exp(-gamma*(t-t0)))
fit Nv(x) "usa.dat" using 1:2 via Nmax
replot Nv(x) title "ヴェアフルストモデル" lc 7
# 凡例の位置を微調整して描きます。
set key default
set key at 1880,290
replot
フーリエ級数展開¶
周期$2L$の周期関数$f(x)$のフーリエ級数展開は
$$ f(x) = \frac{a_0}{2} + \sum_{n=1}^{\infty}\left( a_n \cos\frac{n\pi x}{L} + b_n \sin \frac{n\pi x}{L}\right) $$ここで，
$$ a_n = \frac{1}{L} \int_{-L}^{L} f(x)\cos\frac{n\pi x}{L}\, dx, \quad b_n = \frac{1}{L} \int_{-L}^{L} f(x)\sin\frac{n\pi x}{L}\, dx. $$矩形波のフーリエ級数展開¶
以下のように $\pi \le x < \pi$ で定義された矩形が，この領域外では周期 $2\pi$ の矩形波となっている例。
$$ f(x) = \left\{ \begin{array}{rl} -1 & \mbox{$(-\pi \le x < 0)$} \\ 1 & \mbox{$(0 \le x < \pi)$} \end{array} \right. $$gnuplot でこの周期的矩形波をどうやって表すかが，腕の見せどころ。
私としては，以下を提案する。
$$ f(x) = \mbox{sgn}(\sin x) $$ここで，$\mbox{sgn}(x)$ は，$x>0$ のとき $+1$，$x<0$ のとき $-1$ を返す関数。
reset
f(x) = sgn(sin(x))
plot [-2*pi:2*pi] [-1.5:1.5] f(x)
# 少し体裁を整える。
set samples 500
set grid
set xtics pi
set format x '%4.1P π'
replot
f1(x) = 4./pi*sin(x)
f3(x) = 4./(pi*3)*sin(3*x)
f5(x) = 4./(pi*5)*sin(5*x)
f7(x) = 4./(pi*7)*sin(7*x)
plot [-2*pi:2*pi] [-1.5:2] f(x), \
f1(x) lc 7 lw 3
plot [-2*pi:2*pi] [-1.5:2] f(x), \
f1(x) notitle dt 3, \
f1(x) + f3(x) lc 7 lw 3
plot [-2*pi:2*pi] [-1.5:2] f(x), \
f1(x) notitle dt 3, \
f1(x) + f3(x) notitle dt 3, \
f1(x) + f3(x) + f5(x) lc 7 lw 3
plot [-2*pi:2*pi] [-1.5:2] f(x), \
f1(x) notitle dt 3, \
f1(x) + f3(x) notitle dt 3, \
f1(x) + f3(x) + f5(x) notitle dt 3, \
f1(x) + f3(x) + f5(x) + f7(x) lc 7 lw 3
三角波のフーリエ級数展開¶
以下のように $\pi \le x < \pi$ で定義された三角形が，この領域外では周期 $2\pi$ の三角波となっている例。
$$ f(x) =\left\{ \begin{array}{rl} \pi + x, &\mbox{$(-\pi \le x < 0)$} \\ \pi - x, &\mbox{$(0 \le x < \pi)$} \end{array} \right. $$gnuplot でこの周期的三角波をどうやって表すかが，腕の見せどころ。
私としては，以下を提案する。
f(x) = abs(x - 2*pi*floor(x/(2*pi)) - pi)
reset
set samples 500
set grid
set xtics pi
set format x '%4.1P π'
set ytics 0.5*pi
set format y '%4.1P π'
set yrange [-0.2:4]
f(x) = abs(x - 2*pi*floor(x/(2*pi)) - pi)
plot [-3*pi:3*pi] f(x)
a0 = pi/2
f1(x) = 4./pi*cos(x)
f3(x) = 4./(9*pi)*cos(3*x)
f5(x) = 4./(25*pi)*cos(5*x)
f7(x) = 4./(49*pi)*cos(7*x)
plot [-3*pi:3*pi] f(x), a0 + f1(x) lc 7 lw 2
plot [-3*pi:3*pi] f(x), a0 + f1(x) + f3(x) lc 7 lw 2
plot [-3*pi:3*pi] f(x), a0 + f1(x) + f3(x) + f5(x) lc 7 lw 2
plot [-3*pi:3*pi] f(x), a0 + f1(x) + f3(x) + f5(x) + f7(x) lc 7 lw 2