mathhashimotoの日記 (original) (raw)

関数のグラフを描くときに、零点を知りたいときがある。ところが、0 になっているか若干微妙な状況のときには何らかの工夫をする必要がある。（あるいはsolve()などにより解いてしまうのもひとつの考えである。）本記事では、グラフを使う場合の解決策について考察する。

対数

対数を使うのは安直かつ自然な方法だと思われる。ただし、関数の値が負の場合に対応するために
\operatorname{sign}(f(x)) \cdot \log_{10} |f(x)|
のように符号を掛けることにする。コードとしては

my_log(f,y):=block( [a,b], a:ev(f,[x=y]), b:if a<0 then -1 else 1, return(b*log(abs(a))/log(10)) );

のように実装した。ここでfは変数xについての式である。例：

f:12x^3-12x^2-3*x+4; plot2d([%,my_log(%,x)],[x,-2,2],[y,-10,10]);

3次関数のグラフとその「対数化」

グラフより、0.5 < x < 1.0あたりの極小値は負の値ではなさそうということがわかる。実際には、

poly_discriminant(f,x);
diff(f,x);
solve(%,[x]);
map(lambda([X],ev(f,X)),%),factor;
float(%);

=> -864 36x^2-24x-3 [x=-(sqrt(7)-2)/6,x=(sqrt(7)+2)/6] [(7^(3/2)+19)/9,-(7^(3/2)-19)/9] [4.16891768638357,0.05330453583865156]

より当該極小値は 0.05 程度である。もちろん実解が1つしかない3次式なので判別式の値は負である。より凝ったことをしたいのであれば、次のように色分けすることもできる。

my_log_p(f,y):=block( [a,b], a:ev(f,[x=y]), b:if a<0 then -1 else 1, return(blog((a))/log(10)) ); my_log_m(f,y):=block( [a,b], a:ev(f,[x=y]), b:if a<0 then -1 else 1, return(blog((-a))/log(10)) ); plot2d([f,my_log_p(f,x),my_log_m(f,x)],[x,-2,2],[y,-10,10]);

結果：

3次関数のグラフとその「対数化」・色分け

グラフの描画の際にエラーが出た部分は描画されず無視されることを利用している。

シグモイド関数

我々の問題意識としては、関数の絶対値が大きい範囲は注目していないので、シグモイド関数を使うのもひとつの考えだと思われる。ここでは \tanh x を使うことにする。
\tanh x=\dfrac{e^x-e^{-x}}{e^x+e^{-x}}
また、この関数を原点で「ジャンプ」するように改変した関数を使えば、関数の値の符号が変わった場所がわかりやすくなる。

シグモイド関数

そうすると、上と同じfを使うと次のようなグラフができる。

my_jump(x):=if x<0 then x-1 else x+1;

lambda([x0],tanh(ev(f,x=x0))); lambda([x0],my_jump(tanh(ev(f,x=x0)))); plot2d([f,%th(2),%th(1)],[x,-3,3],[y,-6,6]);

f のグラフ・シグモイド関数を使用

対比として、f0 = f-0.1の場合は以下のようになる。

f0:f-0.1; lambda([x0],tanh(ev(f0,x=x0))); lambda([x0],my_jump(tanh(ev(f0,x=x0)))); plot2d([f0,%th(2),%th(1),(1/2)*my_log(f0,x)-4],[x,-3,3],[y,-6,6]);

f0 のグラフ・シグモイド関数および対数を使用

シグモイド関数と対数の組み合わせ

上記の \tanh x を使う場合ではグラフが平たくなりがちなので、対数とシグモイド関数 \tfrac{y}{1+|y|} を組み合わせたものを使ってみる。
2 \cdot \sigma + \sigma \cdot \dfrac{y}{1+|y|}, \quad \sigma=\operatorname{sign}(f(x)), ~ y=\log_{10} |f(x)|
f(x)=x の場合のグラフ：

シグモイド関数と対数の組み合わせ

my_log_sigmoid(f,x0):=block( [a,b,c], a:ev(f,[x=x0]), b:if a<0 then -1 else 1, c:log(abs(a))/log(10), return(2b+bc/(1+abs(c))) );

lambda([x0],my_log_sigmoid(f0,x0)); plot2d([f0,%],[x,-3,3],[y,-6,6]);

f0の場合のグラフ：

f0 のグラフ（シグモイド関数と対数の合成）

x 軸に接している場合

別の例として、グラフが x 軸に接している場合も考えてみる*1。
g=x^6 (16 x+9)^2

g:x^6*(16*x+9)^2;

lambda([x0],-2+tanh(ev(g,x=x0))); lambda([x0],-5+my_jump(tanh(ev(g-1E-7,x=x0)))); plot2d([g,%th(2),%th(1)],[x,-3,3],[y,-7,2]);

g のグラフ・シグモイド関数を使用

もちろん式からも明らかではあるが、グラフから極小値は限りなく 0 に近いらしいということが読み取れる。また、上述のようにシグモイド関数と対数を組み合わせた場合は、次のようになる。

plot2d([g, -3+my_log_sigmoid(g,x), -5+my_log_sigmoid(g-1E-7,x), -2], [x,-3,3],[y,-7,2]);

g のグラフ（シグモイド関数と対数の合成）

2変数の場合

2変数関数は1変数の場合より難しくなる。実際のところ、implicit()により陰関数グラフを描画しようとすると、どうしてもところどころ正確ではないグラフになってしまうようである。次の例を考える。

h:256y^8+288y^7-480x^2y^6+81y^6-1026x^2y^5 +513x^4y^4-486x^2y^4+405x^4y^3-270x^6y^2 +729x^4y^2+243x^6y+81x^8;

まずは素直にimplicit()を使って陰関数グラフを描画してみる。

tt:1; mag:10^0; draw( dimensions=[1200,500], gr2d( title="h", ip_grid=[200,200], implicit(magh,x,-tt,tt,y,-tt,tt) ), gr2d( title="h+epsilon", color=red, ip_grid=[200,200], implicit(mag(h)+1E-6,x,-tt,tt,y,-tt,tt) ), gr2d( title="h-epsilon", color=green, ip_grid=[200,200], implicit(mag*(h)-1E-6,x,-tt,tt,y,-tt,tt) ), columns=3 );

h のグラフ・implicit() 使用; mag=1

原点周辺の様子を知るための補助として、mag=10^5として描画すると、次のようになる。

h のグラフ・implicit() 使用; mag=10^5

さて、plot2d()で等高線を描画することもできるので、その機能も試してみる。

tt:0.5; plot2d([contour,h],[x,-tt,tt],[y,-tt,tt],[grid,400,400]);

等高線 (1)

上述の対数・シグモイド関数を使った場合は以下のようになる。

my_log_sigmoid_xy(f,x0,y0):=block( [a,b,c], a:ev(f,[x=x0,y=y0]), b:if a<0 then -1 else 1, c:log(abs(a))/log(10), return(2b+bc/(1+abs(c))) );

tt:0.5; plot2d([contour,my_log_sigmoid_xy(h,x,y)], [x,-tt,tt],[y,-tt,tt],[grid,200,200]);

等高線 (2)

現時点では違いがよくわからないのだが、draw3d()にも等高線を描画する機能がある（マニュアルも参照のこと）。

tt:1/2; draw3d( explicit(h,x,-tt,tt,y,-tt,tt), contour_levels = {0,1E-8,1E-6,1E-4,1E-2,0.1,0.2,0.3,0.4}, contour = map, surface_hide = true) $

等高線 (3)

パラメータcontour_lebelsにより等高線の高さを指定できる。なお中括弧は maxima では集合を意味する。またcontour=bothとすればhの3次元グラフも同時に表示される。

補足：hの最低次部分は 81 {{y}^{2}} {{( {{y}^{2}}-3 {{x}^{2}}) }^{2}}なので、 y=0も原点における接線である。ところが、この接線に対応した曲線の成分は係数が実数にはならないので、グラフには表れない*2。実際、その既約成分は2つあって、共に原点で非特異であって、その片方を級数表示すると
y=\dfrac{\sqrt{-3}-1}{6} x^2+\dfrac{11 \sqrt{-3}-12}{243} x^4+\dfrac{545 \sqrt{-3}-612}{19683} x^6+\cdots

(sqrt(-3)-1)/(6) * x^2 + (11*sqrt(-3)-12)/(243) * x^4

のようになっている（複素共役をとるともう一方の級数表示になる）。なお、このような議論は、特異点におけるブローアップを考えると見通しよく理解される。

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(3t)*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(3t)-x,sin(t)sin(3t)-y] (%o2) [-x-cos(t)sin(t)^3+3cos(t)^3sin(t),-y-sin(t)^4+3cos(t)^2sin(t)^2] (%o3) [-x-ST^3+3S^3T,-y-T^4+3S^2T^2] (%o4) [-x-ST^3+3S^3T,-y-T^4+3S^2T^2,T^2+S^2-1] (%o5) x^8(y^4+y^3+2x^2y^2-3x^2y+x^4)^2*(y^8-y^7-14x^2y^6+y^6+x^2y^5+81x^4y^4-6x^2y^4+38x^4y^3-224x^6y^2+9x^4y^2-96x^6y+256x^8)^2 (%o6) [x,y^4+y^3+2x^2y^2-3x^2y+x^4,y^8-y^7-14x^2y^6+y^6+x^2y^5+81x^4y^4-6x^2y^4+38x^4y^3-224x^6y^2+9x^4y^2-96x^6y+256x^8]

消去の過程で余分な因子が出てきてしまったが、最後のリストの2番目の4次式が求める定義方程式である。この式は有名な式であって、以下の表記を使うことも多い。
(x^2+y^2)^2+y^3-3 x^2 y
この曲線は原点で3つの接線 y=0, y=\pm \sqrt{3} x を持っている*1。折角なので、余分な因子（リストの3番目の8次式）で決まる曲線と併せて描画してみる。

h3:y^12-12x^2y^10+6x^2y^9+y^9+54x^4y^8-54x^4y^7-9x^2y^7-100x^6y^6+9x^4y^6+162x^6y^5+27x^4y^5+33x^8y^4-54x^6y^4-138x^8y^3-27x^6y^3+72x^10y^2+81x^8y^2-72x^10y+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 になること）はあくまでも必要条件でしかない。例えば
ax-1=bx-1=0
という x についての連立方程式を考えると、これが（共通）解をもつための必要十分条件は
a-b=0, ~ a \neq 0
である。この1つ目の条件が終結式である。eliminate()は引数として与えられた変数のリストの順番通りに変数を消去していくようなので、そこで結果は変わる可能性がある。実際、S,Tの順で消去を行うと、次の結果になる。

eliminate(%,[S,T]),factor; => [43046721x^4(y^4+y^3+2x^2y^2-3x^2y+x^4)^2 (256y^8+288y^7-480x^2y^6+81y^6-1026x^2y^5+513x^4y^4-486x^2y^4+405x^4y^3-270x^6y^2+729x^4y^2+243x^6y+81*x^8)^2]

なお、本来は変数消去はグレブナ基底というか項順序を使って行う。そうすれば上記のような余分な項が出てくることはない（はず）。maxima のパッケージとしては "grobner" がある。

load(grobner);

としてパッケージを読み込む。例えば次のようにすると、変数消去が実行される。

poly_elimination_ideal(%,2,[S,T,x,y]); => [y^4+y^3+2x^2y^2-3x^2y+x^4]

引数の2は変数リストの2番目までの変数を消去するという意味である。

極座標表示

三つ葉グラフは、最初の三角関数を使った表示においてはpolar()により描画することもできる。

draw2d( nticks=200, polar(sin(3*t),t,0,%pi) );

もちろん、この形の方がコードとしては簡単で見やすい（出力結果は同じなので省略）。さて、それでは上述のh2のグラフについて考える。まずは式を極座標に変換する。

ev(h2,[x=rcos(t),y=rsin(t)]); trigsimp(%),factor; eq_r:%/r^6; divide(%,r^2); => [576sin(t)^8-1872sin(t)^6+2289sin(t)^4-1248sin(t)^2+256, r*(132sin(t)^7-363sin(t)^5+326sin(t)^3-96sin(t)) +16sin(t)^6-24sin(t)^4+9*sin(t)^2]

最後の行は苦し紛れだが、結果として出てきたrの2次式をrについて整理するために割り算を使った。ともかく、rがtについて解けたことになるので、それを使ってグラフを表示する。

solve(eq_r,[r]); [rt1,rt2]:map(rhs,%);

式としては
-\frac{132 {{\sin{(t)}}^{7}}+\sqrt{4-5 {{\sin{(t)}}^{2}}} \left( 4 \cdot {{3}^{\frac{5}{2}}} {{\sin{(t)}}^{6}}-59 \sqrt{3} {{\sin{(t)}}^{4}}+8 \cdot {{3}^{\frac{3}{2}}} {{\sin{(t)}}^{2}}\right) -363 {{\sin{(t)}}^{5}}+326 {{\sin{(t)}}^{3}}-96 \sin{(t)}}{1152 {{\sin{(t)}}^{8}}-3744 {{\sin{(t)}}^{6}}+4578 {{\sin{(t)}}^{4}}-2496 {{\sin{(t)}}^{2}}+512}
のようになっている（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 のグラフ（上述の赤いグラフの拡大）

補足：このようなグラフを描画するときに、「点」としてグラフの断片が出てくる可能性を検討しなければならない。例えば、
x^2+y^2=0
の「グラフ」は原点のみである*2。いつものように（？）陰関数定理を使うことにより、このような点は特異点でなければならないことがわかる（ただし十分条件ではない）。実際に特異点を maxima に見つけてもらうことにすると、

[h2,diff(h2,x),diff(h2,y)]; solve(%,[x,y]); sing_pts:map(lambda([X],ev([x,y],X)),%); => [[0,0],[-5/(32^(3/2)),5/3],[5/(32^(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()によりリストのネスト構造を「平坦化」しているが、しなくても自動的に適切に解釈してくれる仕組みになっているようだ。

点からの射影

上述の計算では、三角関数を使って曲線をパラーメータ表示したが、大体同じやり方だが純代数的な手法を用いることもできる。つまり、（例えば）原点をとって、原点を通る直線の傾きをパラメータとする。言い換えると、グラフの各点 (x,y) に対して、比 (x:y) を対応させるということである。このような対応を点からの射影と呼ぶ。傾きをsとして、

ev(h2,[y=sx],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次式、つまり(s+2)(s-2)になる。射影直線の2点分岐2重被覆はまた射影直線になるので、その座標をパラメータ
u=\sqrt{-\dfrac{s+2}{s-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=sx],factor); h_s:%/x^9; poly_discriminant(h_s,x),factor; => -27s^8*(s^2-3)^8*(s^6-6s^4+9s^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 の前後で値をジャンプさせている。出力結果は問題なさそうである。

本記事は内容を追加していく予定である。出力結果は適当に編集している場合がある。

(1) lisp のリストとしての式

冪や絶対値の中身を取り出すときは、maxima における式は lisp のリストとして実装されていることを利用する。inpart()関数によりリストの要素を取り出す。例：

x^2+2*x*y+y^2;
f:factor(%);
inpart(f,0);
inpart(f,1);
inpart(f,2);

(%o1) y^2+2xy+x^2 (%o2) (y+x)^2 (%o3) "^" (%o4) y+x (%o5) 2

あるいは

x^2+2*x*y+y^2;
g:sqrt(factor(%));
inpart(g,0);
inpart(g,1);

(%o1) y^2+2xy+x^2 (%o2) abs(y+x) (%o3) abs (%o4) y+x

(2) 代数的数の分母の有理化

フラグalgebraicがtrueであるときratsimp()により代数的数の分母を有理化することができる。なお、次のようにするとフラグの変更は1行限りで、それ以降はフラグは元の値に戻る。

(1+sqrt(2))^-1; ratsimp(%), algebraic:true; algebraic; 1 (%o1) ----------- sqrt(2) + 1 (%o2) sqrt(2) - 1 (%o3) false

ただし、これはインタラクティブモードのときのみ有効なやり方である。バッチファイルの場合は以下のようにすればよい。

1+2^(1/4)+2^(3/4); ev( ratsimp( %^-1 ) , algebraic:true); 3/4 1/4 (%o1) 2 + 2 + 1 1/4 sqrt(2) - 2 - 1 (%o2) - ------------------ 3

なお、algebraic=trueとしても同じである。trueの場合に限っては単にフラグの名前だけでもよい。

(3) 行列の生成

成分が何らかの法則で決まっているのような行列はgenmatrix()により定義できる。なおリストについてはmakelist()という関数がある。例：

lambda([i,j],i^j); genmatrix(%,3,3); => [1, 1, 1] [2, 4, 8] [3, 9, 27]

ここでlambda()はラムダ式である。第1引数は変数のリストである。genmatrix()については、第1引数はラムダ式など成分の法則、第2、第3引数は行列のサイズを指定している。もうひとつの例：

lambda([i,j],if i=j then 1 else 0); genmatrix(%,2,3); => [1, 0, 0] [0, 1, 0]

いわゆるクロネッカーのデルタを使っている。

(5) リスト操作：apply() と map()

例えばリストの和を取りたいときは

(%i1) apply("+",makelist(ix^i,i,0,10)); (%o1) 10x^10+9x^9+8x^8+7x^7+6x^6+5x^5+4x^4+3x^3+2x^2+x

のようにapply()を使うとよい。他の例：

min(1,2,3);
v:[1,2,3];
apply(min,v);

(%o1) 1 (%o2) [1,2,3] (%o3) 1

わざわざmin(v[1],v[2],v[3])とする必要はない。別の例：

apply(matrix,[[1,2],[3,4]]);

(%o1) matrix( [1, 2], [3, 4] )

一方で、ある関数に対して、リストの各要素にその関数を作用させて、その結果のリストを作りたいときはmap()を使う。例：

pp(x):=x+1;
map(pp,[1,2,3]);
map(lambda([x],x+1),[4,5,6]);

(%o1) pp(x):=x+1 (%o2) [2,3,4] (%o3) [5,6,7]

最後の行はラムダ式を使っている。こうすると1行にまとめられる。

(6) kill()

kill(all)とすると、ユーザー定義の変数等が全てリセットされる。なお、"all" の制限として、kill(labels)とすれば入出力の番号（ラベル）がリセットされる。なお入出力の番号自体は変数linenumに格納されている。

(7) 出力関連

フラグshowtimeをtrueにすると、コマンドを実行する度に、経過時間が表示される。

本記事では、主にパラメータ付けされた3次元空間中の曲面を扱う。ところで maxima が実際に使っているのは gnuplot なので、凝ったグラフ等を出力したい場合は、gnuplot を勉強する方がいいと思われる。例えば、parametric_surface()により曲面を表示することもできるが、本記事では扱わない。

例1

まず、曲面というよりは曲線だが、以下の簡単な例から始める。この式は螺旋曲線を定義する。

p1:[cos(t),sin(t),t]; plot3d(p1,[t,0,32%pi],[s,0,1],[grid,100,1]);

螺旋曲線

どうやら、3次元空間内の曲線をそのまま描画することはできないようで、上記コードのように、「ダミー」変数sを導入している。gridは変数t, sがそれぞれ動く区間を何分割するかを明示的に指定している。参考：
忘却の微分方程式(28) 3次元プロット、MathematicaとMaxima | デバイスビジネス開拓団

例2

次に、これもよくある例だが、球面を描画する。uは緯度、vは経度である。z 軸方向を地軸としている。なお、画像の通り、maxima では座標軸は右手系でz 軸正の向きが上方向、x 軸正の向きとy 軸正の向きがそれぞれ手前方向と奥方向である。マウスでドラッグすると視点を変えることもできる。また、単なる感想だが、このようなグラフィックが無料のソフトで気軽に表示できるというのは、大昔の自分なら夢のようなことである。

/* u: longitude ; v: latitude */ p2:[cos(u)*cos(v),cos(u)*sin(v),sin(u)]; plot3d(p2,[u,-%pi/2,%pi/2],[v,-%pi,%pi],same_xyz);

球面

オプションsame_xyzは x, y, z の各軸のスケールを一致させる。折角なので面積要素 dS も考えると、
(dS)^2 = (dx \wedge dy)^2 + (dy \wedge dz)^2 +(dz \wedge dx)^2
の uv 平面への引き戻しは以下のように求められる。

Xu:diff(p2,u); Xv:diff(p2,v); XX:matrix(Xu,Xv)$ XX.transpose(XX); trigsimp(%); => [1, 0], [0, cos(u)^2]

従って dS = \cos(u) \, du \wedge dv を得る。これの積分値は単位球の表面積 4 \pi である。なお、trigsimp()は三角関数の式をなるべく簡単にする関数である（ようである）。

例3

主半径R=3、従半径r=1のトーラス*1を描画する。主半径方向の角はu、従半径方向の角はvとしている。

p3:[(R+rcos(v))cos(u),(R+rcos(v))sin(u),rsin(v)]; ev(p3,[R=3,r=1]); plot3d(%,[u,0,2%pi],[v,0,2*%pi],same_xyz);

トーラス

この例でも dS を求めると、

Xu:diff(p3,u); Xv:diff(p3,v); XX:matrix(Xu,Xv)$ XX.transpose(XX); factor(trigsimp(%)); => [(r*cos(v)+R)^2, 0], [0, r^2]

より、 dS=(R+r \cos(v)) r \, du \wedge dv を得る。従ってトーラスの表面積は (2 \pi R) \times (2 \pi r) に等しい（長方形の面積風の式にしてみた）。

*1:「トーラス」が厳密に何を意味するかは状況によると思われる。今の場合は、3次元ユークリッド空間の2次元部分多様体としての構造を考えている。

lisp のリストの長さ

多項式の項数を取得するためには、lisp のlength()を使う。（正確にはこの関数の返り値から1を引く。）maxima から lisp の関数を参照する場合は?length()のように?を付ければよい。なお、例えば?length(x)とすると、リストではないというエラーが出るので、atom()により場合分けして対処する必要がある。自分の勉強のため、もう少し応用ができるようなやり方も考えてみる。あらかじめ .lisp ファイルを作っておいて、

(defun len_lisp (x) (length x))

のように書いておく。そうすると、load()で読み込めば?len_lisp()として参照できる。あるいは、手動で

:lisp (defun $foo (x) (length x))

とすれば、foo()として参照できる。lisp 側で$を付けると maxima 側では?は付けない。ただし、:lispはbatchload()で読み込む maxima バッチファイル内では使えないことに注意。

計算例

最後の節に、自前の関数series_trunc()のコードを書いておいた。変数のリストと、各変数の冪の範囲を指定する。例：

(1+x+x^2)(1+y+y^2); series_trunc([x,y],[[1,2],[0,1]],%); => x^2y+x*y+x^2+x

xの冪指数1から2、yの冪指数0から1の範囲に該当する項だけを集めている。

公式を確かめてみる

それでは早速、テータ関数の公式を適当に考えて、実験してみる。とりあえず |n| \leq 5 の範囲で打ち切ることにして、

N:5; makelist(q^(i^2)y^i,i,-N,N)$ t00:SumL(%); makelist((-1)^iq^(i^2)*y^i,i,-N,N)$ t01:SumL(%); makelist(q^((i+1/2)^2)*y^i,i,-N,N)$ t10:SumL(%);

とする。ここでqは \exp(\pi i \tau)、yは \exp(2 \pi i z) を表す。またSumL()はリストの和をとる自前の関数である。そうすると、公式
\vartheta_{00}(\tau,0)^4-\vartheta_{01}(\tau,0)^4-\vartheta_{10}(\tau,0)^4=0
を確かめるためには

series_trunc([q],[[0,N^2]],t00^4-t10^4-t01^4); ev(%,y=1);

とすればよい。結果は 0 となる。なお、代入y=1の前の式ではqはどの項でも奇数乗として現れる。なお、4乗の計算は時間がかかるので、「2乗の2乗」として計算した方が時間を節約できるだろう。さて、別の公式
2 \vartheta_{00}(2\tau,z)^2
= (\vartheta_{00}(4 \tau) + \vartheta_{10}(4 \tau)) \vartheta_{00}(\tau,z) + (\vartheta_{00}(4 \tau) - \vartheta_{10}(4 \tau)) \vartheta_{01}(\tau,z)
を確かめるために

2*ev(t00,[q=q^2])^2 -ev(t00+t10,[q=q^4,y=1])*t00 -ev(t00-t10,[q=q^4,y=1])*t01$ series_trunc([q],[[0,N^2]],%);

とすると、やはり結果は 0 となる。なお、maxima ではこのような代入をすると、例えば|q|*q^4のように絶対値が出てくることがあるので、それを考慮しなければならない。また、上記計算においてはqの冪指数は常に非負であることに注意する。そうでなければ有限の範囲で打ち切るという操作と掛け算は整合しない場合がある。（環論的にいうなら、打ち切る範囲がイデアルになっていると都合がいいということ。）

コード

あまり高等なコードではないが、一応コードは以下の通り。もちろん最適化も考えていない。rat_mono_pow()は

rat_mono_pow(x,x^(2/3)*y^2); => 2/3

のように、単項式について指定された変数の冪指数を求める関数である。

len_lisp(x):=?length(x);

rat_mono_pow(var,f):=block( [f0,n,k,t,e], if numberp(f) then return(0), /* to apply 'partition()', 'inpart(f,0)' should be "" / if (not atom(f)) and inpart(f,0)="" then( f0:partition(f,var)[2] )else( f0:f ), if (not atom(f0)) and inpart(f0,0)="" then( n:len_lisp(f0)-1, t:inpart(f0,1) )else( n:-1, t:f0 ), k:1, /* the k-th factor of 'f0' / e:0, / exponent / do( if (not atom(t)) and inpart(t,0)=abs then( t:inpart(t,1) ), if atom(t) then( e:e+(if t=var then 1 else 0) )else( / should be of "^" type such as abs(x)^3 */ e:e+inpart(t,2) ), if n=-1 or k=n then return(), k:k+1, t:inpart(f0,k) ), return(e) );

series_trunc(vars,ranges,f):=block( [m,f0,g,flag_mono,flag_range,t,N,k,e], m:length(vars), f0:expand(f), g:0, /* resulting series (polynomial) / if atom(f0) or inpart(f0,0)#"+" then( flag_mono:true, t:f0 )else( flag_mono:false, t:inpart(f0,1), / 1st term / N:len_lisp(f0)-1 / # of terms / ), k:1, / the k-th term of 'f0' / do( flag_range : for i:1 thru m do( e:rat_mono_pow(vars[i],t), if e < ranges[i][1] or ranges[i][2] < e then( return() ), if not atom(t) then( if inpart(t,0)="" then( t:partition(t,vars[i])[1]vars[i]^e )else( / type "^" or "abs" */ t:vars[i]^e ) ) ), if flag_range=done then g:g+t, if flag_mono=true or k=N then return(), k:k+1, t:inpart(f0,k) ), return(g) );

ヤコビのテータ関数（の平行移動）は
\theta_{ab}(\tau,z) = \sum_n \exp(2 \pi i ((n+a)^2 \tau/2+(n+a) (n+z)))
により定まる。ここで a,b \in \mathbb{R} である。maxima でテータ関数を計算する。

MM:100; fpprec:50;

theta_ab(a,b,tau,z):=block( [s], s:0, for i:-MM thru MM do( s:s+bfloat(exp(2*%pi*%i*((i+a)^2tau/2+(i+a)(z+b)))) ), return(expand(s)) );

ここで fpprec は小数の精度を定めている。bfloat() により、この精度で近似値を算出する。また、項をどこまで足すかということで、とりあえず -100 \leq n \leq 100 の範囲としているが、本来は動的に範囲を決めるべきところではある。ただし以下に見るように、この程度で十分な精度が出るようである（もちろん引数の値による）。なお、この程度なら一瞬で計算は終わる*1。例として、公式
\theta_{00}^4(\tau,0) = \theta_{1/2,0}^4(\tau,0) + \theta_{0,1/2}^4(\tau,0)
を確認してみる。

tau0:2*%i; t00:theta_ab(0,0,tau0,0); t10:theta_ab(1/2,0,tau0,0); t01:theta_ab(0,1/2,tau0,0); t00^4-t10^4-t01^4;

結果は

-4.0091470651382934692108047007272228182227953201576b-51

まさに、今の精度で 0 とみなせる。他には、以下のように

(t00^4+t10^4)/(t00^2*t10^2); => 6.0b0

テータ定数の比の \tau=2 i についての特殊値に一致している。別の例：

tau0:rectform(2*%i/(2*%i+1)); t00:theta_ab(0,0,tau0,0); t10:theta_ab(1/2,0,tau0,0); t01:theta_ab(0,1/2,tau0,0); rectform( (t00^4+t10^4)/(t00^2*t10^2) );

ここで rectform() は複素数を x+yi の形に変換する関数である。これも結果は

6.0000000000000000000000000000000000000000000000001b0-3.2073176521106347753686437605817782545782362561261b-50*%i

で、同じ値である。最後に

tau0:%i; s00:theta_ab(0,0,tau0,0); s10:theta_ab(1/2,0,tau0,0); s01:theta_ab(0,1/2,tau0,0); s00/s10; bfloat(2^(1/4)); => 1.1892071150027210667174999705604759152929720924638b0 1.1892071150027210667174999705604759152929720924638b0

となり、完全に一致している。

*1:自分は昔から maxima を使っているが、やはりパソコンの進化によりかなり計算が速くなっていると感じる。

2倍公式