gnuplot を使って,関数のグラフを描くことができます。また,数値データもグラフにすることができます。
gnuplot の2次元グラフ作成でプロットできる対象は以下のとおりです。
- 陽関数 $y = f(x)$:
plot [xmin:xmax] f(x)
- 陰関数 $f(x, y) = 0$ を直接 plot する機能はないようです。
- 媒介変数表示 $x(t), y(t)$:
set parametric; plot [tmin:tmax] x(t), y(t)
- 点,$x$ 座標 $y$ 座標の数値データ:
- ファイルから読み込んで
plot "filename"
- 配列のプロット
plot X using 1:(Y[$1])
など
- ファイルから読み込んで
- ベクトル
- ファイルから読み込んで
plot "filename" w vec
- 配列のプロット
plot X u (X[$1]):(Y[$1]):(Vx[$1]):(Vy[$1]) w vec
など
- ファイルから読み込んで
また,2本の陽関数で挟まれた領域を塗りつぶすこともできます。
陽関数のグラフ
例として,$y = \sin x$ のグラフを $−2\pi \leq x \leq 2\pi$ の範囲で描きます。基本的な定数の一つである円周率 $\pi$ は gnuplot では pi
と書きます。
reset
plot [-2*pi:2*pi] sin(x)
いくつかオプションを設定する例。
# オプションの設定例
# グリッド
set grid
# x 軸 y 軸の表示
set zeroaxis
# 線の太さ:lw (linewidth),線の色:lc (linecolor),凡例: title
# グラフの x 描画範囲と x 軸の表示範囲を別々に設定する例
# sample と書くのがキモ
plot [-10:10][-1.1:1.1] sample \
[-2*pi:2*pi] sin(x) lw 2 lc "red" title "正弦関数"
媒介変数表示のグラフ
半径 $1$ の円の方程式は $x^2 +y^2 = 1$ です。この円を以下のような媒介変数表示にして描きます。
$$ x=\cos t, \quad y=\sin t \quad (0 \leq t \leq 2\pi) $$
reset
set parametric
plot [0:2*pi] cos(t), sin(t)
いくつかオプションを設定する例。
reset
set parametric
# グリッド
set grid
# x 軸 y 軸の表示
set zeroaxis
# 縦横比を 1:1 に
set size square
# 表示範囲
set xrange [-1.1:1.1]
set yrange [-1.1:1.1]
plot [0:2*pi] cos(t), sin(t) notitle
点・数値データのグラフ
以下の例では,配列 X
に $x$ 座標の値,配列 Y
に $y$ 座標の値を入れて,6個の点のグラフを描きます。
reset
array X[6]
array Y[6]
do for [i=1:6]{
X[i] = i-1
Y[i] = (i-1)**2
}
# plot X, Y としたいところですが,それは別の意味。
plot X using 1:(Y[$1])
# plot Y using (X[$1]):2 でも可
いくつかオプションを指定してプロットする例。ここでは, w lp (with linespoints)
で各点を線でつないて描きます。
pointtype (pt) 0~15
with (w) points (p)
lines (l)
linspoints (lp)
linecolor (lc) 0~8 くらい,または "red" 等
plot Y using (X[$1]):2 w lp pt 7 lc 6 notitle
ベクトルを描く
原点を始点とし,成分が $v_x = 0.5, v_y = 1$ のベクトルを描く例。with vector (w vec)
をつけます。
まず,始点のx座標,y座標,x成分,y成分が書かれたファイルを読み込んでベクトルを描く例。
# ファイルを作成する
# 始点のx座標,y座標,x成分,y成分
set print "v.txt"
print 0, 0, 0.5, 1
set output
reset
# x 軸 y 軸の表示
set zeroaxis
# 縦横比を 1:1 に
set size square
# ファイルを読み込んでベクトルを描く
plot [-2:2][-2:2] "v.txt" w vec notitle
special file-name "+"
を使うと,データファイルを作らずにベクトル場を描くことができる。またベクトルの矢がしょぼいので,filled head
をつけてみる。
plot [-2:2][-2:2] "+" using (0):(0):(0.5):(1) w vec filled head notitle
複数のベクトルやベクトル場を描く例。斜方投射の速度ベクトルを描いてみます。まずは,速度ベクトルの始点と成分が書かれたファイルをつくり,それを読み込んでベクトルを描く例。
x(t) = 2*t
y(t) = 2*t - t**2
vx(t) = 0.5
vy(t) = 0.5 - 0.5*t
# 速度ベクトルの始点と成分が書かれたファイルを作成
set print "vec.txt"
do for [i=1:5]{
t = 0.5*(i-1)
print x(t), y(t), vx(t), vy(t)
}
set print
# ファイルの内容確認
!cat vec.txt
reset
set xrange [-0.5:4.5]
set yrange [-0.5:2]
set size ratio 0.5
set grid
plot "vec.txt" w vec filled head notitle
x(t) = 2*t
y(t) = 2*t - t**2
vx(t) = 0.5
vy(t) = 0.5 - 0.5*t
# 速度ベクトルの始点と成分の配列を作成
array X[5]
array Y[5]
array Vx[5]
array Vy[5]
do for [i=1:5]{
t = 0.5*(i-1)
X[i] = x(t)
Y[i] = y(t)
Vx[i] = vx(t)
Vy[i] = vy(t)
}
# 配列の内容を確認
print X, Y, Vx, Vy
reset
set xrange [-0.5:4.5]
set yrange [-0.5:2]
set size ratio 0.5
set grid
plot X using (X[$1]):(Y[$1]):(Vx[$1]):(Vy[$1]) \
w vec filled head notitle
さらには以下のようなベクトル場のグラフも描くことができます。参考までに。
複数のグラフを重ねて表示
gnuplot で複数のグラフを重ねて表示する例を示します。
$y = x^2 – 1$ と $y = 4 x – 5$ の 2 つのグラフを重ねて描きます。
$x$ の範囲は $-5 \leq x \leq 5$ です。
reset
plot [-5:5] x**2 - 1, 4*x - 5
いくつかオプションを設定して描く例です。線の太さの設定は
linewidth (lw) 数値
# x 軸 y 軸の表示
set zeroaxis
# グリッドの表示
set grid
# x 軸の目盛は 1 ごとに
set xtics 1
# y 軸の目盛は 1 ごとに
set ytics 1
# y の表示範囲を [-2:6] に。
# 1行が長くなる場合は \ で改行
plot [-5:5][-2:6] x**2 - 1 lw 2 title "x^2 - 1", \
4*x - 5 title "4 x - 5"
数値データファイルを読み込む
Maxima では,あらかじめ作成された数値データファイルを読み込んでグラフを描くこともできます。
reset
# 読み込むファイルを作成
set print "data.txt"
do for [i=0:5]{
print i, i**2
}
set print
# ファイル data.txt の内容を表示
! cat data.txt
# 凡例の表示位置
set key left top
plot [0:6][0:28] "data.txt" pt 7 lc "red" title "データ"
数値データと理論曲線を重ねて表示
前節の数値データファイル data.txt
と理論曲線 $y = x^2$ の 2 つのグラフを重ねて表示してみます。
plot [0:6][0:28] "data.txt" pt 7 lc "red" title "データ", \
x**2 lw 2 lc "blue" title "理論曲線 y=x^2"
海王星と冥王星の軌道
焦点を原点とした楕円のグラフ
太陽からの万有引力を受けて運動する惑星は,太陽(二体問題では質量中心)を焦点の一つとした楕円運動を描いて運動します。
焦点の一つを原点とし,長半径 $a$,離心率 $e$ の楕円の方程式は,極座標 $r$,$\varphi$ を使って以下のように表すことができます。
$$ r= \frac{a (1−e^2)}{ 1 + e\cos\varphi}$$
さて,かつては第9惑星,現在では準惑星の一つである冥王星も楕円軌道を描きます。冥王星の軌道長半径 $𝑎_P=39.445 \mbox{au}$,離心率 $𝑒_P=0.250$
を使って冥王星の軌道を描きます。
まず,極座標表示の楕円の式を関数として定義します。
r(a, e, phi) = a*(1-e**2)/(1 + e*cos(phi))
x(a, e, phi) = r(a, e, phi)*cos(phi)
y(a, e, phi) = r(a, e, phi)*sin(phi)
媒介変数表示でプロット
冥王星の軌道長半径と離心率を入れて,まずは媒介変数表示でプロットしてみます。
reset
aP = 39.445
eP = 0.250
set zeroaxis
set xrange [-50:50]
set yrange [-50:50]
set size square
set parametric
plot [phi=0:2*pi] x(aP, eP, phi), y(aP, eP, phi) title "冥王星"
では,この冥王星の軌道と,その内側をまわる海王星の軌道を重ねて描いてみましょう。
海王星の軌道長半径は $a_N=30.181 \mbox{au}$,離心率は $e_N=0.0097$ と小さいので簡単のために $e_N=0$ として扱います。
実際の2つの天体の軌道は同じ平面上にありませんが,太陽からの距離のみをみるために,ここでは同一平面上に描きます。
aN = 30.181
eN = 0
replot x(aN, eN, phi), y(aN, eN, phi) title "海王星"
set output "gplot12.svg"
replot
set output
参考:極座標表示でプロット
極座標表示でプロットする場合は,set polar
とします。(解除は unset polar
です。)
reset
set zeroaxis
set xrange [-50:50]
set yrange [-50:50]
set size square
set polar
plot [phi=0:2*pi] r(aP, eP, phi) title "冥王星", \
r(aN, eN, phi) title "海王星"
月別平年気温の関数フィット
弘前市の月別平年気温のデータを使って配列 HiroT
を作ります。
array HiroT[12] = [-1.5,-1.0,2.3,8.6,14.3,18.3,\
22.3,23.5,19.4,12.9,6.5,0.8]
# 読み込むファイルを作成
set print "hiro.txt"
do for [i=1:12]{
print i, HiroT[i]
}
set print
# ファイルの内容表示
! cat hiro.txt
reset
set pointsize 0.6
set grid
set xtics 1
plot [1:12] "hiro.txt" using 1:2 pt 7 title '月別平年気温'
グラフをみると,月別平年基本は 12ヶ月周期の正弦関数あるいは余弦関数のようにみえます。では,以下のような関数でフィットしてみます。
$$f(x) = a + b \sin\left(\frac{2\pi x}{12}\right)
+ c \cos\left(\frac{2\pi x}{12}\right)$$
f(x) = a + b*sin(pi*x/6) + c*cos(pi*x/6)
fit f(x) "hiro.txt" using 1:2 via a, b, c
print a, b, c
最小二乗法によって求めた $a, b, c$ の値をもった $f(x)$ のグラフを重ねてプロットしてみます。
set key sample 1
replot f(x) title "関数フィット"
領域の塗りつぶし
$y = f(x)$, $y = g(x)$, $x_1 \leq x \leq x_2$ の領域を塗りつぶす例。書式は…
plot [x1:x2] f(x) with filledcurves (w filledc) y=g(x) fillcolor (fc) ...
f(x) = 0.6*x + 0.4*cos(x)
reset
set zeroaxis
plot [0:2.5][-0.2:1.5] sample \
[0.5:2] f(x) w filledc y=0 fc "yellow" notitle, \
[0.25:2.25] f(x) lw 2 lc "blue" title "f(x)"