Maxima を使って，関数のグラフを描くことができます。また，数値データもグラフにすることができます。

Maxima の2次元グラフ作成は draw2d() でもできますが，ここでは，plot2d() について説明します。
plot2d() でプロットできるオブジェクトは以下のとおりです。

陽関数 $y = f(x)$:
- plot2d(f(x), [x, xmin, xmax])
陰関数 $f(x, y) = 0$:
- plot2d(f(x, y) = 0, [x, xmin, xmax], [y, ymin, ymax])
媒介変数表示 $x(t), y(t)$:
- plot2d([parametric, x(t), y(t), [t, t0, t1]])
点，$x$ 座標 $y$ 座標の数値データ:
- plot2d([discrete, [[x1, y1], ..., [xn, ..., yn]]]) または
- plot2d([discrete, [x1, ..., xn], [y1, ..., yn]])

陽関数のグラフ

例として，$y = \sin x$ のグラフを $−2\pi \leq x \leq 2\pi$ の範囲で描きます。基本的な定数の一つである円周率 $\pi$ は Maxima では %pi と書きます。

In [1]:

plot2d(sin(x), [x, -2*%pi, 2*%pi])$

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

In [2]:

/* オプションの設定例 */
plot2d(sin(x), [x, -2*%pi, 2*%pi], 
  /* 表示範囲 */
  [x, -10, 10], [y, -1.1, 1.1], 
  /* グリッド */
  grid2d,   
  /* 凡例 */
  [legend, "正弦関数"], 
  /* 線の太さと色 */
  [style, [lines, 2, red]]
)$

陰関数のグラフ

例として，$x^2 + y^2 = 1$ のグラフを描きます。$y$ について解いて陽関数の形にしたり，以下で述べるような媒介変数表示にしなくても，陰関数のまま，グラフに描くことができます。

陰関数 $f(x, y) = 0$ を [x, xmin, xmax], [y, ymin, ymax] の範囲で描く書式は…

plot2d(f(x, y) = 0, [x, xmin, xmax], [y, ymin, ymax])$

In [3]:

plot2d(x**2 + y**2 - 1 = 0, [x, -1.2, 1.2], [y, -1.2, 1.2])$

縦横比・サイズの設定

縦横の比 yx_ratio を設定して円らしく見えるようにします。size も変更します。弘大 JupyterHub ではグラフを svg でインライン表示させているので，　gnuplot_svg_term_command で size 600, 600 などと設定します。

In [4]:

/* オプションの設定例 */
plot2d(x**2 + y**2 = 1, [x, -1.2, 1.2], [y, -1.2, 1.2], 
  /* 縦横比 [yx_ratio, 1] または... */
  [same_xy], 
  [gnuplot_svg_term_command, "set term svg font \",14\" size 600,600"]
)$

媒介変数表示のグラフ

半径 $1$ の円の方程式は $x^2 +y^2 = 1$ です。この円を以下のような媒介変数表示にして描きます。

$$ x=\cos t, \quad y=\sin t \quad (0 \leq t \leq 2\pi) $$

書式は…

plot2d([parametric, x(t), y(t), [t, tmin, tmax]])$

In [5]:

plot2d([parametric, cos(t), sin(t), [t, 0, 2*%pi]]
)$

In [ ]:

In [6]:

/* オプションの設定例 */
plot2d([parametric, cos(t), sin(t), [t, 0, 2*%pi]],   
  /* 縦横比 [same_xy] は縦横の表示範囲によらず円は円に見える */
  [same_xy], 
  /* 表示範囲 */
  [x, -1.5, 1.5], [y, -1.1, 1.1], 
  /* 横軸・縦軸のラベル */
  [xlabel, "x"], [ylabel, "y"]
)$

点・数値データのグラフ

以下の例では，最初（左側）の数値を $x$ 座標の値，次（右側）の数値を $y$ 座標の値として 6 個の点からなるリスト data の各点をつないだ折れ線グラフを描きます。

In [7]:

data: makelist([i, i^2], i, 0, 5);

Out[7]:

\[\tag{${\it \%o}_{13}$}\left[ \left[ 0 , 0 \right] , \left[ 1 , 1 \right] , \left[ 2 , 4 \right] , \left[ 3 , 9 \right] , \left[ 4 , 16 \right] , \left[ 5 , 25 \right] \right] \]

In [8]:

plot2d([discrete, data])$

以下のように，$x$ 座標のみの数値データ，$y$ 座標のみの数値データをプロットしてもよいです。

In [9]:

X: makelist(i, i, 0, 5);
Y: makelist(i**2, i, 0, 5);

Out[9]:

\[\tag{${\it \%o}_{16}$}\left[ 0 , 1 , 2 , 3 , 4 , 5 \right] \]

Out[9]:

\[\tag{${\it \%o}_{17}$}\left[ 0 , 1 , 4 , 9 , 16 , 25 \right] \]

オプションを設定して，各点を線でつながずにプロットする例。

ここでは，[style, points]を設定して線でつながずに点で描きます。プロットオプションの style には，points の他にも，lines，linespoints， dots が設定できます。

In [10]:

/* オプションの設定例 */
plot2d([discrete, X, Y],   
  /* 表示範囲 */
  [x, 0, 5], [y, 0, 28],
  /* 点と大きさ */
  [style, [points, 2]], 
  /* 色 blue, red, green, magenta, black, cyan から */
  [color, red]
)$

複数のグラフを重ねて表示

Maxima で複数のグラフを重ねて表示する例を示します。

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

$x$ の範囲は $-5 \leq x \leq 5$ です。
以下の例でわかるように，重ねて描きたい関数を [x^2 - 1, 4*x - 5] のように [ と ] で囲んだリスト形式にして plot2d コマンドを使います。

In [11]:

plot2d([x^2 - 1, 4*x - 5], [x, -5, 5])$
system("cp ~/.maxplot.svg ./maxplot09.svg")$

いくつかオプションを設定して描く例です。

In [12]:

plot2d(
  [x^2 - 1, 4*x - 5], [x, -5, 5],

  /* 表示範囲 */
  [y, -5, 10],
  /* 線の色の設定 */ 
  [color, red, black],   
  /* 線の太さの設定。最初の線の太さを 3 に */
  [style, [lines, 3], lines],    
  /* 凡例の位置。グラフの線とかぶらないように */
  [gnuplot_preamble, "set key bottom"] 
)$

数値データファイルを読み込む

Maxima では，あらかじめ作成された数値データファイルを読み込んでグラフを描くこともできます。

In [13]:

/* 読み込むデータファイルを作成しておく */
tmp: makelist([i, i^2], i, 0, 5)$
write_data(tmp, "data.txt")$

データファイルを読み込む Maxima の関数は，read_nested_list を使います。

In [14]:

dat: read_nested_list("data.txt");

Out[14]:

\[\tag{${\it \%o}_{26}$}\left[ \left[ 0 , 0 \right] , \left[ 1 , 1 \right] , \left[ 2 , 4 \right] , \left[ 3 , 9 \right] , \left[ 4 , 16 \right] , \left[ 5 , 25 \right] \right] \]

リスト dat に読み込んだ数値データを plot2d() でグラフにします。

In [15]:

plot2d([discrete, dat],
       [x, 0, 6], [y, 0, 28],        /* x, y の表示範囲を設定 */ 
       [style, [points, 1.5]], [color, blue],   /* 大きさ，色 */
       [gnuplot_preamble, "set key top left"],  /* 凡例の位置 */ 
       [legend, "データ"]                   /* 凡例表示の設定 */
       )$

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

前節の数値データを読み込んだリスト dat と理論曲線 $y = x^2$ の 2 つのグラフを重ねて表示してみます。

In [16]:

/* もう少しオプションを指定して描く */
plot2d([[discrete, dat], x**2], [x, 0, 5], 
       [style, [points, 1.5], lines], 
       [color, blue, red],
       [gnuplot_preamble, "set key top left"], 
       [legend, "データ", "理論曲線 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$
を使って冥王星の軌道を描きます。

参考：冥王星 – Wikipedia

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

In [17]:

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);

Out[17]:

\[\tag{${\it \%o}_{31}$}r\left(a , e , \varphi\right):=\frac{a\,\left(1-e^2\right)}{1+e\,\cos \varphi}\]

Out[17]:

\[\tag{${\it \%o}_{32}$}x\left(a , e , \varphi\right):=r\left(a , e , \varphi\right)\,\cos \varphi\]

Out[17]:

\[\tag{${\it \%o}_{33}$}y\left(a , e , \varphi\right):=r\left(a , e , \varphi\right)\,\sin \varphi\]

媒介変数表示でプロット

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

In [18]:

aP: 39.445$
eP: 0.250$

In [19]:

plot2d([parametric, x(aP, eP, phi), y(aP, eP, phi), [phi, 0, 2*%pi]],
        [same_xy],
        [xlabel, "x"], [ylabel, "y"], [legend, "冥王星"],
        [x,-50, 50], [y, -50, 50]
      )$

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

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

参考：海王星 – Wikipedia

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

In [20]:

aN: 30.181$
eN: 0$

In [21]:

plot2d([[parametric, x(aP, eP, phi), y(aP, eP, phi), [phi, 0, 2*%pi]],
        [parametric, x(aN, eN, phi), y(aN, eN, phi), [phi, 0, 2*%pi]]],
        [same_xy],
        [xlabel, ""], [ylabel, ""], 
        [legend, "冥王星", "海王星"],
        [x,-50, 50], [y, -50, 50]
      )$

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

弘前市の月別平年気温のデータを使ってリスト HiroT を作ります。

参考：弘前　平年値（年・月ごとの値）　主な要素

In [22]:

HiroT: [
[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 [23]:

plot2d([discrete, HiroT], [style, [points, 1]])$

グラフを見ると，月別平年気温は 12 ヶ月周期の三角関数でフィットできるようにみえます。

では，以下のように関数フィットをしてみましょう。

まず，最小二乗法のパッケージ lsquares をロードします.これは lsquares_estimates を初めて使う際，最初に 1 回だけ行えばよいです。

In [24]:

load(lsquares)$

次に，次のような関数を定義して，最小二乗法でパラメータ a, b, c を求めます。

$$f(x) = a + b \sin\left(\frac{2\pi}{12}\right)+ c \cos\left(\frac{2\pi}{12}\right)$$

In [25]:

f(x):= a + b*sin(%pi*x/6)+ c*cos(%pi*x/6);

Out[25]:

\[\tag{${\it \%o}_{46}$}f\left(x\right):=a+b\,\sin \left(\frac{\pi\,x}{6}\right)+c\,\cos \left(\frac{\pi\,x}{6}\right)\]

最小二乗法を行う関数 lsquares_estimates には数値データを行列として入れてやる必要があります。

このために apply('matrix, HiroT) でリスト HiroT を行列に変換しています。

In [26]:

/* 表示は 5 桁くらいで... */
fpprintprec: 5$

lsquares_estimates(apply('matrix, HiroT), [x, y], y=f(x), [a, b, c])$
float(%);

Out[26]:

\[\tag{${\it \%o}_{49}$}\left[ \left[ a=10.533 , b=-8.3403 , c=-9.1611 \right] \right] \]

最小二乗法によって求めた a, b, c の値を f(x) に代入して fit とし，数値データと合わせて，この関数のグラフも描きます。

In [27]:

fit: ev(f(x), %);

Out[27]:

\[\tag{${\it \%o}_{50}$}-8.3403\,\sin \left(\frac{\pi\,x}{6}\right)-9.1611\,\cos \left(\frac{\pi\,x}{6}\right)+10.533\]

In [28]:

plot2d([[discrete, HiroT], fit], [x, 1, 12],  
       [style, [points, 1], lines], 
       /* 表示範囲 */ [x, 0.8, 12.2],
       [color, blue, red],
       [legend, "月別平年気温", "関数フィット"], 
       /* x 軸の目盛は 1 ごとに。凡例は center に */
       [xtics, 1], [gnuplot_preamble, "set key center"],
       [xlabel, "月"], [ylabel, "°C"]
     )$

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本の陽関数で挟まれた領域を塗りつぶすこともできます。

gnuplot, Maxima, Python グラフ作成・比較

Maxima によるグラフ作成 plot2d 編

陽関数のグラフ

陰関数のグラフ

縦横比・サイズの設定

媒介変数表示のグラフ

点・数値データのグラフ

複数のグラフを重ねて表示

数値データファイルを読み込む

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

海王星と冥王星の軌道

焦点を原点とした楕円のグラフ

媒介変数表示でプロット

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

Python によるグラフ作成 SymPy Plotting Backends 編

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

SymPy Plotting Backends で正規分布をσごとに塗りわける

matplotlib.pyplot.plot で2曲線の描画範囲を個別に設定して plot と塗りつぶし

SymPy Plotting Backends で2曲線の描画範囲を個別に設定して plot と塗りつぶし

グラフの上下左右軸に主目盛・副目盛をつける

SymPy Plotting Backends でグラフを描いて Matplotlib 流にオプション設定する

SymPy Plotting Backends を使ってみる

このページの目次