Return to Maxima でコンピュータ演習

Maxima でグラフ作成 plot2d 編

Maxima の2次元グラフ作成は draw2d() でもできますが，ここでは，plot2d() について説明します。

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 [3]:
/* オプションの設定例 */
plot2d(
  sin(x), [x, -2*%pi, 2*%pi], 

  /* 表示範囲 */
  [x, -10, 10], [y, -1.1, 1.1], 
  /* x 軸 y 軸の表示・非表示 true/false */
  [axes, false], 
  /* グリッド */
  grid2d, 
  
  /* 凡例 */
  [legend, "正弦関数"], 
  /* 線の太さと色 */
  [style, [lines, 2, red]]
)$

陰関数のグラフ

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

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

In [7]:
/* オプションの設定例 */
plot2d(
  x**2 + y**2 = 1, 
  [x, -1.2, 1.2], [y, -1.2, 1.2], 

  /* 縦横比 [same_xy] または... */
  [yx_ratio, 1]
)$

媒介変数表示のグラフ

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

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

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

In [11]:
/* オプションの設定例 */
plot2d(
  [parametric, cos(t), sin(t), [t, 0, 2*%pi]], 
  
  /* 縦横比 [yx_ratio, 1] または... */
  [same_xy], 
  /* 表示範囲 */
  [x, -1.1, 1.1], [y, -1.1, 1.1], 
  /* 横軸・縦軸のラベル */
  [xlabel, ""], [ylabel, ""]
)$

点・数値データのグラフ

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

In [13]:
data: makelist([i, i^2], i, 0, 5);
Out[13]:
\[\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 [14]:
plot2d([discrete, data])$

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

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

In [16]:
/* オプションの設定例 */
plot2d(
  [discrete, data], 
  
  /* 表示範囲 */
  [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 [18]:
plot2d([x^2 - 1, 4*x - 5], [x, -5, 5])$

いくつかオプションを設定して描く例です。

In [20]:
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"] 
)$

上のグラフをみると，2 つのグラフは 1 点で接しているように見えます。

実際に $y = 4x – 5$ が接線であることを示してみましょう。

$x^2 – 1 = 4 x – 5$ の解が1つであることを示せば十分ですね。

In [22]:
sol: solve(x^2 - 1 = 4*x - 5, x);
Out[22]:
\[\tag{${\it \%o}_{22}$}\left[ x=2 \right] \]

$y = f(x)$ の $x = x_1$ における接線の方程式は，
$$ y – f(x_1) = f'(x_1) \,(x – x_1)$$ ですから…

In [23]:
f(x):= x^2 - 1;
define(df(x), diff(f(x), x));

x1: rhs(sol[1]);
y = df(x1) * (x - x1) + f(x1), ratsimp;
Out[23]:
\[\tag{${\it \%o}_{23}$}f\left(x\right):=x^2-1\]
Out[23]:
\[\tag{${\it \%o}_{24}$}{\it df}\left(x\right):=2\,x\]
Out[23]:
\[\tag{${\it \%o}_{25}$}2\]
Out[23]:
\[\tag{${\it \%o}_{26}$}y=4\,x-5\]

練習 3-1

$y = f(x) = x^2 -1$ の $x = 1$ における接線，および $ x = 3$ における接線の方程式を求め，これら3本の線をグラフに描け。

In [24]:
f(x):= x^2 - 1;
define(df(x), diff(f(x), x));
Out[24]:
\[\tag{${\it \%o}_{27}$}f\left(x\right):=x^2-1\]
Out[24]:
\[\tag{${\it \%o}_{28}$}{\it df}\left(x\right):=2\,x\]

$y = f(x)$ の $x = x_1$ における接線の方程式は，
$$ y – f(x_1) = f'(x_1) \,(x – x_1)$$ ですから…

In [ ]:
In [ ]:
数値データファイルを読み込む
Maxima では，あらかじめ作成された数値データファイルを読み込んでグラフを描くこともできます。
以下のような内容のファイル test.dat がカレントディレクトリにあるとします。（ない場合は，弘大 JupyterHub の Home で「新規」から「テキストファイル」を選択し，以下の数値をコピペし，「test.dat」にリネームし，保存してください。）
0   0
1   1
2   4
3   9
4   16
5   25


データファイルを読み込む Maxima の関数は，read_nested_list を使います。
In [25]:
dat: read_nested_list("test.dat");

Out[25]:
\[\tag{${\it \%o}_{29}$}\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 [26]:
plot2d([discrete, dat],
       [x, 0, 6], [y, 0, 28],        /* x, y の表示範囲を設定 */ 
       [style, [points, 1]], [color, red],      /* 大きさ，色 */
       [gnuplot_preamble, "set key top left"],  /* 凡例の位置 */ 
       [legend, "データ"]                   /* 凡例表示の設定 */
       )$

数値データと理論曲線を重ねて表示
前節の数値データファイル test.dat と理論曲線 $y = x^2$ の 2 つのグラフを重ねて表示してみます。
In [28]:
/* もう少しオプションを指定して描く */
plot2d([[discrete, dat], x**2], [x, 0, 5], 
       [style, [points, 1.5], lines], 
       [gnuplot_preamble, "set key top left"], 
       [legend, "データ", "理論曲線 y=x^2"]
       )$

練習 3-2
楕円のグラフ（楕円中心を原点にして）
媒介変数表示のグラフを複数重ねて表示する例です。
原点を中心とし，長半径 $a$，単半径 $b$ の楕円の式は
$$\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1$$
媒介変数表示では，$x=a \cos t, \ y=b \sin t, \ (0\leq t\leq 2\pi)$と書けます。
 $a$，$b$ に適当な値を入れて複数の楕円のグラフを重ねて描きなさい。
In [30]:
/* ヒント：媒介変数 t 表示の x, y を関数として定義 */
x(a, t):= a*cos(t)$
y(b, t):= b*sin(t)$

In [ ]:
2つの楕円のグラフを重ねて描く場合は，
[parametric, x(2,t), y(1,t), [t, 0, 2*%pi]]
`

の部分を，例えば $a = 2, b = 1$ の楕円と $a = 1, b = 1.5$ の楕円として
[[parametric, x(2,t), y(1,t), [t, 0, 2*%pi]], 
 [parametric, x(1,t), y(1.5,t), [t, 0, 2*%pi]]
]

のようにします。
In [ ]:
海王星と冥王星の軌道
焦点を原点とした楕円のグラフ
太陽からの万有引力を受けて運動する惑星は，太陽（二体問題では質量中心）を焦点の一つとした楕円運動を描いて運動します。
 焦点の一つを原点とし，長半径 $a$，離心率 $e$ の楕円の方程式は，極座標 $r$，$\varphi$ を使って以下のように表すことができます。
$$ r= \frac{a (1−e^2)}{ 1 + e\cos\varphi}$$
さて，かつては第9惑星，現在では準惑星の一つである冥王星も楕円軌道を描きます。冥王星の軌道長半径 $𝑎_P=39.445 \mbox{au}$，離心率 $𝑒_P=0.250$
 を使って冥王星の軌道を描きます。
まず，極座標表示の楕円の式を関数として定義します。
In [31]:
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[31]:
\[\tag{${\it \%o}_{36}$}r\left(a , e , \varphi\right):=\frac{a\,\left(1-e^2\right)}{1+e\,\cos \varphi}\]
Out[31]:
\[\tag{${\it \%o}_{37}$}x\left(a , e , \varphi\right):=r\left(a , e , \varphi\right)\,\cos \varphi\]
Out[31]:
\[\tag{${\it \%o}_{38}$}y\left(a , e , \varphi\right):=r\left(a , e , \varphi\right)\,\sin \varphi\]
冥王星の軌道長半径と離心率：
In [32]:
aP: 39.445$
eP: 0.250$

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

では，この冥王星の軌道と，その内側をまわる海王星の軌道を重ねて描いてみましょう。
海王星の軌道長半径は $a_N=30.181 \mbox{au}$，離心率は $e_N=0.0097$ と小さいので簡単のために $e_N=0$ として扱います。
実際の2つの天体の軌道は同じ平面上にありませんが，太陽からの距離のみをみるために，ここでは同一平面上に描きます。
In [35]:
aN: 30.181$
eN: 0$

In [36]:
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]]
       ],
        [yx_ratio, 1],
        [xlabel, "x"], [ylabel, "y"], [legend, "冥王星", "海王星"],
        [x,-50, 50], [y, -50, 50]
      )$

さて，この冥王星と海王星の軌道のグラフをいると，軌道が交差しているように見えます。このことを確かめます。
表示される領域を制限して，交差しているあたりを拡大してみます。
In [38]:
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]]
       ],
        [yx_ratio, 1],
        [xlabel, "x"], [ylabel, "y"], [legend, "冥王星", "海王星"],
        [x, 26, 28], [y, 12, 14]
      )$

海王星の軌道は半径 $a_N$ の円としていますから，以下の式を満たす角度 $\varphi$
 を求めればよいことになります。
$$ r(a_P, e_P, \varphi) = a_N$$
$0 < \varphi < \pi/2$
 のどこかで交差しているようですから，このあたりで数値的に解を求める方法は以下の通りです。
In [40]:
phi1: find_root(r(aP, eP, t) = aN, t, 0, %pi/2);

Out[40]:
\[\tag{${\it \%o}_{49}$}0.4485997043415258\]
$ -\pi/2 < \varphi < 0$ の範囲でも交差していますね。
 対称性から答えは phi1 に負号をつけたものになるのは明らかですが，念のために確かめます。
In [41]:
phi2: find_root(r(aP, eP, t) = aN, t, -%pi/2, 0);

Out[41]:
\[\tag{${\it \%o}_{50}$}-0.4485997043415258\]
求めた角度はラジアン単位ですから，度に直してみます。
In [42]:
fpprintprec: 4$ /* 表示するケタ数を 4 に */ 
float(phi1*180/%pi);
% * 2;

Out[42]:
\[\tag{${\it \%o}_{52}$}25.7\]
Out[42]:
\[\tag{${\it \%o}_{53}$}51.41\]
つまり一周 360 度のうち，この角度にあたる分は，海王星のほうが冥王星より太陽から遠い位置にあることがわかります。
原点と 2 つの軌道の交点を結ぶ直線を描く方法は以下の通りです。
In [43]:
plot2d([[discrete, [[0, 0], [aN*cos(phi1), aN*sin(phi1)]]],
        [discrete, [[0, 0], [aN*cos(phi2), aN*sin(phi2)]]]
       ],
        [yx_ratio, 1], [x,-50, 50], [y, -50, 50], 
        [xlabel, "x"], [ylabel, "y"], [legend, "", ""],
        [color, black, black]
      )$

練習 3-3
原点と2つの軌道の交点を結ぶ直線を追加して，冥王星と海王星の軌道を描け。
In [ ]:
月別平年気温の関数フィット
弘前市の月別平年気温のデータを使ってリスト HiroT を作ります。
In [45]:
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]
];

Out[45]:
\[\tag{${\it \%o}_{56}$}\left[ \left[ 1 , -1.5 \right] , \left[ 2 , -1.0 \right] , \left[ 3 , 2.3 \right] , \left[ 4 , 8.6 \right] , \left[ 5 , 14.3 \right] , \left[ 6 , 18.3 \right] , \left[ 7 , 22.3 \right] , \left[ 8 , 23.5 \right] , \left[ 9 , 19.4 \right] , \left[ 10 , 12.9 \right] , \left[ 11 , 6.5 \right] , \left[ 12 , 0.8 \right] \right] \]
このデータをグラフに描いてみます。
In [46]:
plot2d([discrete, HiroT], [style, [points, 1]])$

In [47]:
/* 弘大 JupyterHub では 
   set_draw_defaults(file_name="~/.maxplot",terminal='svg)$
   されているので。*/
system("cp ~/.maxplot.svg ./maxplot16.svg")$

グラフを見ると，月別平年気温は 12 ヶ月周期の三角関数でフィットできるようにみえます。
では，以下のように関数フィットをしてみましょう。
まず，最小二乗法のパッケージ lsquares をロードします.これは lsquares_estimates を 初めて使う際，最初に 1 回だけ行えばよいです。
In [48]:
load(lsquares)$

次に，次のような関数を定義して，うまくフィットするように定数 a0, a1, b1 を求めます。
In [49]:
f(x) := a0/2 + a1*cos(2*%pi*x/12) + b1*sin(2*%pi*x/12);

Out[49]:
\[\tag{${\it \%o}_{60}$}f\left(x\right):=\frac{a_{0}}{2}+a_{1}\,\cos \left(\frac{2\,\pi\,x}{12}\right)+b_{1}\,\sin \left(\frac{2\,\pi\,x}{12}\right)\]
参考：フーリエ級数展開
一般に周期 $2L$ の周期関数は，以下のようにフーリエ級数展開できるのであった。したがって，ここでやっていることは，12ヶ月分の月平年気温という数値データに対して，フーリエ級数展開みたいなことを数値的に行っていると考えることができる。
$$ f(x) = \frac{a_0}{2} + \sum_{n=1}^{\infty} \bigl( a_n \cos \left(n \frac{\pi}{L} {x}\right) + b_n \sin \left(n\frac{\pi}{L} {x}\right) \bigr) $$
最小二乗法を行う関数 lsquares_estimates には数値データを行列として入れてやる必要があります。
このためにリスト HiroT を行列 M に変換します。（リストと行列では表示が異なります。行列は行列らしく表示されています。）
In [50]:
M: apply('matrix, HiroT);

Out[50]:
\[\tag{${\it \%o}_{61}$}\begin{pmatrix}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 \\ \end{pmatrix}\]
In [51]:
/* 有効数字 3 桁ですから... */
fpprintprec: 3$

lsquares_estimates(M, [x, y], y=f(x), [a0, a1, b1])$
float(%);

Out[51]:
\[\tag{${\it \%o}_{64}$}\left[ \left[ a_{0}=21.1 , a_{1}=-9.16 , b_{1}=-8.34 \right] \right] \]
最小二乗法によって求めた a0, a1, b1 の値を fit に代入し，この関数のグラフを描きます。
In [52]:
fit: ev(f(x), %);

Out[52]:
\[\tag{${\it \%o}_{65}$}-8.34\,\sin \left(\frac{\pi\,x}{6}\right)-9.16\,\cos \left(\frac{\pi\,x}{6}\right)+10.5\]
In [53]:
plot2d(fit, 
       [x, 0.5, 12.5], [y, -5, 28],
       [style,lines], [color, red],
       [gnuplot_preamble, "set xtics 1"],
       [xlabel, "月"], 
       [ylabel, "月別平年気温"], 
       [title, "弘前市の月別平年気温"]
     )$

In [54]:
/* 弘大 JupyterHub では 
   set_draw_defaults(file_name="~/.maxplot",terminal='svg)$
   されているので。*/
system("cp ~/.maxplot.svg ./maxplot17.svg")$

練習 3-4
  1. フィットさせる関数を以下のような g(x) にし，最小二乗法で求めてみよ。
  2. 月別平年気温の数値データとフィット曲線のグラフを重ねて描け。
  3. 月平年気温の平均値を求めてみよ。（g(x) の定数部分 a0/2 が平均値。）
In [55]:
g(x):= a0/2 + a1*cos(2*%pi*x/12) + b1*sin(2*%pi*x/12) 
            + a2*cos(4*%pi*x/12) + b2*sin(4*%pi*x/12)$

In [ ]:
月平年気温のデータを入れたリスト HiroT を使って平均値を求める例。
In [56]:
/* HiroT[i][2] が i 月の気温だから... */
sum(HiroT[i][2], i, 1, length(HiroT))/length(HiroT);

Out[56]:
\[\tag{${\it \%o}_{69}$}10.5\]
月平年気温のデータを入れた行列 M を使って平均値を求める例。
In [57]:
/* 行列 M の2列目の気温を取り出す */
col(M, 2);

Out[57]:
\[\tag{${\it \%o}_{70}$}\begin{pmatrix}-1.5 \\ -1.0 \\ 2.3 \\ 8.6 \\ 14.3 \\ 18.3 \\ 22.3 \\ 23.5 \\ 19.4 \\ 12.9 \\ 6.5 \\ 0.8 \\ \end{pmatrix}\]
In [58]:
/* 行列 M の i 行 2列目が i 月の気温だから... */
sum(M[i, 2], i, 1, length(M))/length(M);

Out[58]:
\[\tag{${\it \%o}_{71}$}10.5\]