maxima で2次元グラフを描くときはplot2d()を使うのが簡単な方法だが、draw2d()を使うとより複雑なことができる(本来の gnuplot の機能に近い)。基本的には、描画するオブジェクトを順番に記述していく。またオプションも必要に応じて指定する。draw2d()が受け取る引数としては、複数の引数や、リスト形式(複数のリストも可)も受け取ってくれる。
参考:
Introduction to draw (Maxima 5.47.0 Manual)
Maxima でグラフ作成 draw2d 編 - 相対論の理解とその周辺
基本
陽関数、陰関数、パラメトリック 曲線を描画するときは、それぞれexplicit(), implicit(), parametric()を使う。引数は前半が実際の関数、後半が各変数が動く範囲を指定する。なお、極座標 表示を使うpolar()もある(後述)。また、点を表示したいときはpoints()を使う。
cv1:[sin(3*t)*cos(t),sin(3*t)*sin(t)];
draw2d(
nticks=200,
ip_grid=[100,100],
proportional_axes=xy,
explicit(sin(x),x,-%pi,%pi),
implicit(x^2+y^2=1,x,-1.5,1.5,y,-1.5,1.5),
parametric(cv1[1],cv1[2],t,0,%pi)
); 3つのグラフ 複数のグラフを描画したいときは、単に順番に引数にしていけばよい。ここで、nticksはparametric()の滑らかさ(全体の区間 を何分割するか)を、ip_gridはimplicit()のグリッドの細かさを指定している。またproportional_axes=xyは x 軸と y 軸のスケールを一致させるオプションである。
計算例:三つ葉 グラフ
実は、上述のグラフの中で「三つ葉 」グラフは三角関数 を使わずに代数曲線(多項式 )として表すことができる。折角なので定義方程式を求めてみる。(おそらくこの式というかグラフは名前がついているのではないかという気もするが、今のところ不明である。)
cv1_eq:[sin(3*t)*cos(t)-x,sin(3*t)*sin(t)-y];
trigexpand(%),expand;
ev(%,[sin(t)=T,cos(t)=S]);
append(%,[T^2+S^2-1]);
eliminate(%,[T,S])[1],factor;
makelist(inpart(inpart(%,i),1),i,1,3);
(%o1) [cos(t)*sin(3*t)-x,sin(t)*sin(3*t)-y]
(%o2) [-x-cos(t)*sin(t)^3+3*cos(t)^3*sin(t),-y-sin(t)^4+3*cos(t)^2*sin(t)^2]
(%o3) [-x-S*T^3+3*S^3*T,-y-T^4+3*S^2*T^2]
(%o4) [-x-S*T^3+3*S^3*T,-y-T^4+3*S^2*T^2,T^2+S^2-1]
(%o5) x^8*(y^4+y^3+2*x^2*y^2-3*x^2*y+x^4)^2*(y^8-y^7-14*x^2*y^6+y^6+x^2*y^5+81*x^4*y^4-6*x^2*y^4+38*x^4*y^3-224*x^6*y^2+9*x^4*y^2-96*x^6*y+256*x^8)^2
(%o6) [x,y^4+y^3+2*x^2*y^2-3*x^2*y+x^4,y^8-y^7-14*x^2*y^6+y^6+x^2*y^5+81*x^4*y^4-6*x^2*y^4+38*x^4*y^3-224*x^6*y^2+9*x^4*y^2-96*x^6*y+256*x^8] 消去の過程で余分な因子が出てきてしまったが、最後のリストの2番目の4次式が求める定義方程式である。この式は有名な式であって、以下の表記を使うことも多い。
この曲線は原点で3つの接線 , を持っている*1 。折角なので、余分な因子(リストの3番目の8次式)で決まる曲線と併せて描画してみる。
[h1,h2]:[%[2],%[3]];
h3:y^12-12*x^2*y^10+6*x^2*y^9+y^9+54*x^4*y^8-54*x^4*y^7-9*x^2*y^7-100*x^6*y^6+9*x^4*y^6+162*x^6*y^5+27*x^4*y^5+33*x^8*y^4-54*x^6*y^4-138*x^8*y^3-27*x^6*y^3+72*x^10*y^2+81*x^8*y^2-72*x^10*y+16*x^12;
draw2d(
ip_grid=[200,200],
proportional_axes=xy,
line_width=2,
implicit(h1=0,x,-1.5,1.5,y,-1.5,1.5),
color=red,
implicit(h2=0,x,-2,2,y,-2,2),
color=green,
line_width=1,
implicit(h3=0,x,-2,2,y,-3,2)
); 3つの代数曲線 ここでh1,h2に同時に代入を行っている。colorは当然色を指定している。line_widthは線の太さである。なおh3は上記計算をしているときに計算間違いで偶然出てきた式だが、ついでに描画してみた。原点に「特異点 」があるので、その周辺での描画が変になる可能性があるが、仕方がない。(何となく映画の羊たちの沈黙 の檻が逆光になっているあのシーンのように見えなくもない。)実は、この赤と緑のグラフも本来は原点を通っているはずなので、この出力はおかしいということになる。これは maxima (gnuplot ) の内部でほぼ 0 の値を厳密な 0 とみなしてしまっていることにより起きている。
今後の課題: 本来はこの 0 とみなす閾値 を変更することで、より正確なグラフの描画を実現するべきところだが、うまくいかなかった。つまり、gnuplot だとset zero <value>として当該閾値 の変更をする。draw2d()の場合はオプションuser_preamble="<command>"により gnuplot にコマンドを渡すことで、当該閾値 を変更できるはずであるが、結果は変わらなかった。とりあえず応急処置としては以下のような対策がある。どうやら、当該閾値 は 1e-6 と設定されているようである。したがって、描画する関数に 1e-6 あるいは 9.99e-7 などを加えれば、「正の方向」からの当該閾値 の効果を打ち消すことができると思われる。結果は以下の通り(ついでにグリッドも細かくしてみた)。
draw2d(
ip_grid=[300,300],
proportional_axes=xy,
line_width=2,
implicit(h1+1E-6,x,-1.5,1.5,y,-1.5,1.5),
color=red,
implicit(h2+1E-6,x,-2,2,y,-2,2),
color=green,
line_width=1,
implicit(h3+1E-6,x,-2,2,y,-3,2)
); 改良版・3つの代数曲線 緑のグラフは「はさみグラフ」と(勝手に)名付けることにした。注意: 実はこのグラフに「点」を付け加える必要がある(後述)。
補足:消去について
上記計算でなぜ余分な因子が出てくるかというと、eliminate()は終結 式を使って変数を1つずつ消去していくが、終結 式(が 0 になること)はあくまでも必要条件でしかない。例えば
という についての連立方程式 を考えると、これが(共通)解をもつための必要十分条件 は
である。この1つ目の条件が終結 式である。eliminate()は引数として与えられた変数のリストの順番通りに変数を消去していくようなので、そこで結果は変わる可能性がある。実際、S,Tの順で消去を行うと、次の結果になる。
eliminate(%,[S,T]),factor;
=>
[43046721*x^4*(y^4+y^3+2*x^2*y^2-3*x^2*y+x^4)^2
*(256*y^8+288*y^7-480*x^2*y^6+81*y^6-1026*x^2*y^5+513*x^4*y^4-486*x^2*y^4+405*x^4*y^3-270*x^6*y^2+729*x^4*y^2+243*x^6*y+81*x^8)^2] なお、本来は変数消去はグレブナ基底というか項順序を使って行う。そうすれば上記のような余分な項が出てくることはない(はず)。maxima のパッケージとしては "grobner" がある。
load(grobner); としてパッケージを読み込む。例えば次のようにすると、変数消去が実行される。
poly_elimination_ideal(%,2,[S,T,x,y]);
=>
[y^4+y^3+2*x^2*y^2-3*x^2*y+x^4] 引数の2は変数リストの2番目までの変数を消去するという意味である。
三つ葉 グラフは、最初の三角関数 を使った表示においてはpolar()により描画することもできる。
draw2d(
nticks=200,
polar(sin(3*t),t,0,%pi)
); もちろん、この形の方がコードとしては簡単で見やすい(出力結果は同じなので省略)。さて、それでは上述のh2のグラフについて考える。まずは式を極座標 に変換する。
ev(h2,[x=r*cos(t),y=r*sin(t)]);
trigsimp(%),factor;
eq_r:%/r^6;
divide(%,r^2);
=>
[576*sin(t)^8-1872*sin(t)^6+2289*sin(t)^4-1248*sin(t)^2+256,
r*(132*sin(t)^7-363*sin(t)^5+326*sin(t)^3-96*sin(t))
+16*sin(t)^6-24*sin(t)^4+9*sin(t)^2] 最後の行は苦し紛れだが、結果として出てきたrの2次式をrについて整理するために割り算を使った。ともかく、rがtについて解けたことになるので、それを使ってグラフを表示する。
solve(eq_r,[r]);
[rt1,rt2]:map(rhs,%); 式としては
のようになっている(2次方程式 なので平方根 のプラスマイナスの違いがある)。実はpolar()はrが虚数 の範囲でも実部を取って描画を行うので、(おそらく自動化もできるとは思うが)そこは手作業で除く必要がある。
solve(4-5*sin(t)^2,[t]);
float(%);
t0:ev(t,%[2]);
=>
solve: using arc-trig functions to get a solution.
Some solutions will be lost.
[t=-asin(2/sqrt(5)),t=asin(2/sqrt(5))]
[t=-1.10714871779409,t=1.10714871779409]
1.10714871779409 というわけで、t0=1.10...として-t0 < t < t0の範囲で描画を行えばよい。
draw(
gr2d(
xrange=[-1,1],
yrange=[-2,1.5],
nticks=200,
proportional_axes=xy,
color=grey,
line_width=5,
ip_grid=[200,200],
implicit(h2=0,x,-2,2,y,-2,2),
color=blue,
line_width=1,
polar(rt1,t,-t0,t0),
polar(rt2,t,-t0,t0)
),
gr2d(
xaxis=true,
dimensions=[1500,500],
xrange=[-1.2,1.2],
yrange=[-2,2],
explicit(rt1,t,-t0,t0),
explicit(rt2,t,-t0,t0)
),
columns=2
); h2 のグラフ 右のグラフは極座標 表示である。このように、複数のグラフを並べて表示することもできる。その場合はdraw()を使う。引数は個別のグラフ(scene オブジェクトとよばれる)やオプションである。この例ではcolumns=2としているが、つまり1行に2個並べるという設定である。個別のグラフはdraw2d()の代わりにgr2d()により生成する。つまりdraw2d()とはdraw()とgr2d()を組み合わせたものである。xrange,yrangeは描画の範囲、dimensionsは出力結果(全体)の大きさを指定している。元の陰関数による描画(灰色)と比べてみると、原点以外ではほぼ同じであることがわかる。また、上述の改良版の出力結果とは概ね一致している。ついでだが、描画の範囲を変更することで原点周辺を拡大した場合は次のようになる。
h2 のグラフ(拡大) 上述の赤いグラフの拡大:
h2 のグラフ(上述の赤いグラフの拡大) 補足: このようなグラフを描画するときに、「点」としてグラフの断片が出てくる可能性を検討しなければならない。例えば、
の「グラフ」は原点のみである*2 。いつものように(?)陰関数定理を使うことにより、このような点は特異点 でなければならないことがわかる(ただし十分条件 ではない)。実際に特異点 を maxima に見つけてもらうことにすると、
[h2,diff(h2,x),diff(h2,y)];
solve(%,[x,y]);
sing_pts:map(lambda([X],ev([x,y],X)),%);
=>
[[0,0],[-5/(3*2^(3/2)),5/3],[5/(3*2^(3/2)),5/3]] となり、原点を除く2点をグラフに追加しなければならないことがわかる。
draw2d(
xrange=[-1,1],
yrange=[-2,2],
nticks=200,
proportional_axes=xy,
polar(rt1,t,-t0,t0),
polar(rt2,t,-t0,t0),
color=red,
point_type=2,
point_size=2,
points(makelist(sing_pts[i],i,2,3))
); h2 のグラフ(完成) これでh2のグラフは完成したことになる。
点オブジェクト
上述のようにpoints()により点オブジェクトを生成できる。引数の形式は3通りある(マニュアル参照)。オプションとしては、point_typeにより点を表す記号を、point_sizeにより記号の大きさを指定できる。記号は以下の通りである*3 。
pts:flatten(
makelist([point_type=i,points([[i,1]])],i,0,15)
);
draw2d(
xrange=[-1,16],
yrange=[0.8,1.2],
point_size=3,
pts
); point_type なおflatten()によりリストのネスト構造を「平坦化」しているが、しなくても自動的に適切に解釈してくれる仕組みになっているようだ。
点からの射影
上述の計算では、三角関数 を使って曲線をパラーメータ表示したが、大体同じやり方だが純代数的な手法を用いることもできる。つまり、(例えば)原点をとって、原点を通る直線の傾きをパラメータとする。言い換えると、グラフの各点 に対して、比 を対応させるということである。このような対応を点からの射影と呼ぶ。傾きをsとして、
ev(h2,[y=s*x],factor);
h_s:%/x^6;
poly_discriminant(%,x),factor;
=>
-3*(s-2)*s^4*(s+2)*(s^2-8)^2*(s^2-3)^2 となる。h_sはxの2次式なので、この式をxについて解くためには判別式の平方根 をとればよい。偶然だが、判別式は2乗の因子を除くと2次式、つまり になる。射影直線の2点分岐2重被覆はまた射影直線になるので、その座標をパラメータ
とすればよい。
u^2=-(s+2)/(s-2);
s_u:solve(%,s);
ev(h_s,%,factor);
solve(%,[x]);
x_u:factor(%[1]); 2つの解は、uを-uに置き換えることで移り合う。その片方をx_uとした。実際、
ev(y=s*x,[s_u],factor);
y_u:ev(%,[x_u],factor);
ev(h2,[x_u,y_u],factor);
=>
0 となり、グラフのパラメータ付けを与えていることが確認できる。これを使ってグラフを描画してみると、
plot2d([parametric,rhs(x_u),rhs(y_u),[u,-100,100]]); 原点からの射影によるパラメータ付け となる。傾きを使っているので仕方ないが、グラフ描画の細かさにムラがある。また、一周させるには(このままでは)uの範囲を限りなく広くしなければならない。なお、このような式はいつでもあるわけではなく、代数曲線が射影直線と双有理同値(幾何種数が 0 になる)のときに限る。
はさみグラフ
折角なので、上記のh3についても正確なグラフを描いてみる。(怒っているクリーパーの顔グラフと呼んでもいい気もしてきた。なおh3+1E-3=0とすると、怒っている顔のグラフになる。)
ev(h3,[y=s*x],factor);
h_s:%/x^9;
poly_discriminant(h_s,x),factor;
=>
-27*s^8*(s^2-3)^8*(s^6-6*s^4+9*s^2+4)^3 h_sはxの3次式である。詳細は省略するが、判別式は常に0 以下の値をとるので、実解は 1 つだけあることがわかる。
solve(h_s,[x])[3],factor;
x_s:rhs(%);
draw(
dimensions=[1500,500],
gr2d(
xrange=[-1.5,1.5],yrange=[-2.5,1.5],
nticks=800,
parametric(x_s,s*x_s,s,-30,10)
),
gr2d(
ip_grid=[500,500],
implicit(h3+1E-6,x,-1.5,1.5,y,-2.5,1.5)
),
columns=2
); h3 のグラフ・比較 なおこの代数曲線の特異点 は原点のみである。
補足:零の閾値
次のような関数を使うことで、implicit()の 0 の閾値 問題の対策をすることもできなくもないが、いいやり方なのかはよくわからない。少なくとも処理時間は増えた気がする。
my_jump(f,x0,y0):=block(
[a,b],
a:ev(f,[x=x0,y=y0]),
b:if a>0 then 1 else -1,
return(a+b)
);
draw2d(
ip_grid=[300,300],
implicit(my_jump(h3,x,y),x,-3,3,y,-3,3)
); 要は、値が正なら 1 を、負なら -1 を加えることで 0 の前後で値をジャンプさせている。出力結果は問題なさそうである。