gnuplot によるグラフ作成基本編

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 と書きます。

In [1]:
reset
In [2]:
plot [-2*pi:2*pi] sin(x)

いくつかオプションを設定する例。

In [3]:
# オプションの設定例
# グリッド
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) $$

In [4]:
reset
In [5]:
set parametric
plot [0:2*pi] cos(t), sin(t)

Out[5]:
	dummy variable is t for curves, u/v for surfaces

いくつかオプションを設定する例。

In [6]:
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

Out[6]:
	dummy variable is t for curves, u/v for surfaces

点・数値データのグラフ

以下の例では,配列 X に $x$ 座標の値,配列 Y に $y$ 座標の値を入れて,6個の点のグラフを描きます。

In [7]:
reset
In [8]:
array X[6]
array Y[6]
do for [i=1:6]{
  X[i] = i-1
  Y[i] = (i-1)**2
}
In [9]:
# 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" 
In [10]:
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成分が書かれたファイルを読み込んでベクトルを描く例。

In [11]:
# ファイルを作成する
# 始点のx座標,y座標,x成分,y成分
set print "v.txt"
  print 0, 0, 0.5, 1
set output
In [12]:
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 をつけてみる。

In [13]:
plot [-2:2][-2:2] "+" using (0):(0):(0.5):(1) w vec filled head notitle

複数のベクトルやベクトル場を描く例。斜方投射の速度ベクトルを描いてみます。まずは,速度ベクトルの始点と成分が書かれたファイルをつくり,それを読み込んでベクトルを描く例。

In [14]:
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
In [15]:
# ファイルの内容確認
!cat vec.txt
Out[15]:
0.0 0.0 0.5 0.5
1.0 0.75 0.5 0.25
2.0 1.0 0.5 0.0
3.0 0.75 0.5 -0.25
4.0 0.0 0.5 -0.5
In [16]:
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

In [17]:
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)
}
In [18]:
# 配列の内容を確認
print X, Y, Vx, Vy
Out[18]:
[0.0,1.0,2.0,3.0,4.0]
[0.0,0.75,1.0,0.75,0.0]
[0.5,0.5,0.5,0.5,0.5]
[0.5,0.25,0.0,-0.25,-0.5]

In [19]:
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

さらには以下のようなベクトル場のグラフも描くことができます。参考までに。

iwaki

複数のグラフを重ねて表示

gnuplot で複数のグラフを重ねて表示する例を示します。

$y = x^2 – 1$ と $y = 4 x – 5$ の 2 つのグラフを重ねて描きます。

$x$ の範囲は $-5 \leq x \leq 5$ です。

In [20]:
reset
In [21]:
plot [-5:5] x**2 - 1, 4*x - 5

いくつかオプションを設定して描く例です。線の太さの設定は

linewidth (lw) 数値
In [22]:
# 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 では,あらかじめ作成された数値データファイルを読み込んでグラフを描くこともできます。

In [23]:
reset
In [24]:
# 読み込むファイルを作成
set print "data.txt"
  do for [i=0:5]{
    print i, i**2
  }
set print
In [25]:
# ファイル data.txt の内容を表示
! cat data.txt
Out[25]:
0 0
1 1
2 4
3 9
4 16
5 25
In [26]:
# 凡例の表示位置
set key left top

plot [0:6][0:28] "data.txt" pt 7 lc "red" title "データ"

数値データと理論曲線を重ねて表示

前節の数値データファイル data.txt と理論曲線 $y = x^2$ の 2 つのグラフを重ねて表示してみます。

In [27]:
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$
を使って冥王星の軌道を描きます。

まず,極座標表示の楕円の式を関数として定義します。

In [28]:
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)

媒介変数表示でプロット

冥王星の軌道長半径と離心率を入れて,まずは媒介変数表示でプロットしてみます。

In [29]:
reset
In [30]:
aP = 39.445
eP = 0.250
In [31]:
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 "冥王星"

Out[31]:
	dummy variable is t for curves, u/v for surfaces

では,この冥王星の軌道と,その内側をまわる海王星の軌道を重ねて描いてみましょう。

海王星の軌道長半径は $a_N=30.181 \mbox{au}$,離心率は $e_N=0.0097$ と小さいので簡単のために $e_N=0$ として扱います。

実際の2つの天体の軌道は同じ平面上にありませんが,太陽からの距離のみをみるために,ここでは同一平面上に描きます。

In [32]:
aN = 30.181
eN = 0
In [33]:
replot x(aN, eN, phi), y(aN, eN, phi) title "海王星"
set output "gplot12.svg"
replot
set output

参考:極座標表示でプロット

極座標表示でプロットする場合は,set polar とします。(解除は unset polar です。)

In [34]:
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 "海王星"

Out[34]:
	dummy variable is t for curves

月別平年気温の関数フィット

弘前市の月別平年気温のデータを使って配列 HiroT を作ります。

In [35]:
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]
In [36]:
# 読み込むファイルを作成
set print "hiro.txt"
  do for [i=1:12]{
    print i, HiroT[i]
  }
set print
In [37]:
# ファイルの内容表示
! cat hiro.txt
Out[37]:
1 -1.5
2 -1.0
3 2.3
4 8.6
5 14.3
6 18.3
7 22.3
8 23.5
9 19.4
10 12.9
11 6.5
12 0.8
In [38]:
reset
In [39]:
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)$$

In [40]:
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
Out[40]:
iter      chisq       delta/lim  lambda   a             b             c            
   0 2.2382958688e+03   0.00e+00  8.16e-01    1.000000e+00   1.000000e+00   1.000000e+00
   1 1.9206962103e+01  -1.16e+07  8.16e-02    1.003158e+01  -7.406230e+00  -8.144960e+00
   2 4.7566283262e+00  -3.04e+05  8.16e-03    1.053305e+01  -8.339219e+00  -9.159939e+00
   3 4.7566133159e+00  -3.16e-01  8.16e-04    1.053333e+01  -8.340255e+00  -9.161067e+00
iter      chisq       delta/lim  lambda   a             b             c            

After 3 iterations the fit converged.
final sum of squares of residuals : 4.75661
rel. change during last iteration : -3.15567e-06

degrees of freedom    (FIT_NDF)                        : 9
rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 0.726989
variance of residuals (reduced chisquare) = WSSR/ndf   : 0.528513

Final set of parameters            Asymptotic Standard Error
=======================            ==========================
a               = 10.5333          +/- 0.2099       (1.992%)
b               = -8.34026         +/- 0.2968       (3.559%)
c               = -9.16107         +/- 0.2968       (3.24%)

correlation matrix of the fit parameters:
                a      b      c      
a               1.000 
b              -0.000  1.000 
c               0.000  0.000  1.000 
In [41]:
print a, b, c
Out[41]:
10.5333333317856 -8.34025525998032 -9.16106711406752

最小二乗法によって求めた $a, b, c$ の値をもった $f(x)$ のグラフを重ねてプロットしてみます。

In [42]:
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) ...
In [43]:
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)"