2025-04-05

判別イデアル？（判別式の類似）

1変数多項式の判別式は重解を持つための条件だが、もちろん3重解を持つための条件、あるいは重解を2つ持つための条件も考えることはできる。その場合の条件は（単一の）係数の多項式というよりは、イデアルになる（判別イデアルと呼ぶべきか？）。よく知られているように、判別式や終結式は行列式としての表示があるが、判別イデアルのいい表示があるのかはわからない。計算には singular を使った。

例として、5次式を考えてみる。斉次化は後ですることにして、x^5 の係数（元々 0 ではないと仮定）は 1 とする。他の係数には自然な重みをつけておく。なお、coeffs()は、指定した変数についての係数からなる1列の行列を出力する。イデアル型のコンストラクタideal()は行列を受け取ると、全ての成分により生成されるイデアルを返す。次のようなコードが自然であると思われる。

ring r=(0),(x,A,B,a0,a5,a4,a3,a2,a1),(lp(4),wp(5,4,3,2,1));
poly f=a0*x^5+a1*x^4+a2*x^3+a3*x^2+a4*x+a5;
poly g=reduce(f,(x-A)^3);
poly h=reduce(f,(x-A)^2*(x-B)^2);
matrix mg,mh;
mg=coeffs(g,x);
ideal P0=subst(ideal(mg),a0,1);
ideal P=eliminate(P0,A);
mh=coeffs(h,x);
ideal Q0=subst(ideal(mh),a0,1);
ideal Q=eliminate(Q0,AB);

> size(P);
3
> size(Q);
3

ここで、項順序としてa5, ... , a1の順にすると、求めたいイデアルP, Qのグレブナ基底のサイズは小さ目になる。すなわち、両方のイデアルともサイズは3個である。順序を逆にすると、サイズはそれぞれ6個、7個となる：

ring rr=(0),(a1,a2,a3,a4,a5),wp(1,2,3,4,5);
size(std(imap(r,P)));
size(std(imap(r,Q)));
setring r;
=>
6
7

ここでimap()は、環の間でイデアル等を、変数は共通としてやり取りするための関数である。斉次化したいときは、homog()を使う。ただし、この関数は斉次化のための変数を指定しないときは斉次イデアルかどうかを判定する。

> homog(P);
0

> homog(P,a0);
出力略

実験のために、多項式の解のひとつが x=0 であること、つまりa5=0を仮定すると、以下のような結果を得る。何をしているかというと、準素分解をしている（詳しくはマニュアル参照）。

> LIB "primdec.lib";
> list pd1=primdecGTZ(std(subst(P,a5,0)));
> pd1;
[1]:
   [1]:
      _[1]=108*a3^2+32*a2^3-108*a3*a2*a1-9*a2^2*a1^2+27*a3*a1^3
      _[2]=12*a4+a2^2-3*a3*a1
   [2]:
      _[1]=108*a3^2+32*a2^3-108*a3*a2*a1-9*a2^2*a1^2+27*a3*a1^3
      _[2]=12*a4+a2^2-3*a3*a1
[2]:
   [1]:
      _[1]=a3^3
      _[2]=-a3^2+3*a4*a2
      _[3]=a4*a3
      _[4]=a4^2
   [2]:
      _[1]=a3
      _[2]=a4

第1成分は、（生成的には）3重解が x=0 ではない場合、第2成分は、3重解が x=0 である場合に対応している。Qでも似たようなことができるが、省略。

最後に、折角なので計算結果を貼り付ける。

> homog(P,a0);
_[1]=20*a0^2*a4^2-75*a0^2*a5*a3+3*a0*a3^2*a2-9*a0*a4*a2^2-a0*a4*a3*a1+45*a0*a5*a2*a1-a3^2*a1^2+3*a4*a2*a1^2-12*a5*a1^3
_[2]=500*a0^3*a5*a4-45*a0^2*a3^3+160*a0^2*a4*a3*a2-225*a0^2*a5*a2^2-60*a0^2*a4^2*a1-350*a0^2*a5*a3*a1+23*a0*a3^2*a2*a1-69*a0*a4*a2^2*a1-6*a0*a4*a3*a1^2+330*a0*a5*a2*a1^2-6*a3^2*a1^3+18*a4*a2*a1^3-72*a5*a1^4
_[3]=625*a0^2*a5^2-15*a0*a4*a3^2+40*a0*a4^2*a2+25*a0*a5*a3*a2-250*a0*a5*a4*a1-2*a3^2*a2^2+6*a4*a2^3+6*a3^3*a1-19*a4*a3*a2*a1-15*a5*a2^2*a1+9*a4^2*a1^2+40*a5*a3*a1^2
> homog(Q,a0);
_[1]=80*a0^3*a4^2-50*a0^3*a5*a3+27*a0^2*a3^2*a2-36*a0^2*a4*a2^2-54*a0^2*a4*a3*a1+30*a0^2*a5*a2*a1+4*a0*a2^4-18*a0*a3*a2^2*a1+42*a0*a4*a2*a1^2-8*a0*a5*a1^3-a2^3*a1^2+4*a3*a2*a1^3-8*a4*a1^4
_[2]=1000*a0^4*a5*a4+135*a0^3*a3^3-180*a0^3*a4*a3*a2-200*a0^3*a5*a2^2+280*a0^3*a4^2*a1-700*a0^3*a5*a3*a1+20*a0^2*a3*a2^3-9*a0^2*a3^2*a2*a1-68*a0^2*a4*a2^2*a1-132*a0^2*a4*a3*a1^2+460*a0^2*a5*a2*a1^2+12*a0*a2^4*a1-59*a0*a3*a2^2*a1^2+20*a0*a3^2*a1^3+116*a0*a4*a2*a1^3-104*a0*a5*a1^4-3*a2^3*a1^3+12*a3*a2*a1^4-24*a4*a1^5
_[3]=625*a0^5*a5^2+135*a0^4*a4*a3^2+240*a0^4*a4^2*a2-600*a0^4*a5*a3*a2-250*a0^4*a5*a4*a1+108*a0^3*a3^2*a2^2-124*a0^3*a4*a2^3-54*a0^3*a3^3*a1-234*a0^3*a4*a3*a2*a1+360*a0^3*a5*a2^2*a1-71*a0^3*a4^2*a1^2+290*a0^3*a5*a3*a1^2+16*a0^2*a2^5-80*a0^2*a3*a2^3*a1+9*a0^2*a3^2*a2*a1^2+183*a0^2*a4*a2^2*a1^2+62*a0^2*a4*a3*a1^3-250*a0^2*a5*a2*a1^3-8*a0*a2^4*a1^2+36*a0*a3*a2^2*a1^3-8*a0*a3^2*a1^4-70*a0*a4*a2*a1^4+40*a0*a5*a1^5+a2^3*a1^4-4*a3*a2*a1^5+8*a4*a1^6

2025-03-09

singular における適切な項順序、パラメーター付きの場合

項順序の設定

言うまでもなく、singular において（グレブナ基底において）、項順序を適切に定めることが重要であり、項順序によっては計算がいつまでも終わらない*1ということが起こり得る。ただし、groebner()により経験的に計算が早く終わると考えられる項順序でグレブナ基底を求めてくれるので、std()でうまくいかない場合は試してみるといいだろう。また、option(prot)により計算の進捗状況を表示させられるので、その情報も参考になる。

さて、考えているイデアルがパラメーター（定数だとみなしたい不定元）を含むようなときは、項順序の「直和」を考えるというまくいくときがある。例：

ring r=0,(x,y,a0,a1,a2,b0,b1,b2,c0,c1,c2),(dp(2),dp(9));
ideal I=
a0*x^2+a1*x*y+a2*y^2,
b0*x^2+b1*x*y+b2*y^2,
c0*x^2+c1*x*y+c2*y^2;

これは3つの2変数斉次2次式が共通解を持つという条件を求めるためのイデアルのつもりである。こうすると、最初の2変数についての順序と最後の9変数（パラメーター）についての順序を、前者を（辞書式順序で）優先する形で順序を入れたことになる。グレブナ基底を求めると、

> ideal J=std(I);
> size(J);
21
> factorize(J[1]);
[1]:
   _[1]=1
   _[2]=a2*b1*c0-a1*b2*c0-a2*b0*c1+a0*b2*c1+a1*b0*c2-a0*b1*c2
   _[3]=y
[2]:
   1,1,2

というような結果を得る。つまり、自然にグレブナ基底の1つ目では x が（結果的に）消去されている。なお、係数行列の行列式が因子として現れている。

消去について補足

このイデアルでは射影直線の斉次座標 x,y で考えているので、変数 x,y の消去とは、通常の意味の消去ではなく、成分に分解（準素イデアル分解）したときに x=y=0 に含まれるようなものを全て除くことに他ならない。従って、次のようにイデアルの商*2をとるとよい。

intvec v;
ideal K0,K1,K2,K3,K4;
K0=J;
v=v,size(K0);v;K1=quotient(K0,ideal(x,y));
v=v,size(K1);v;K2=quotient(K1,ideal(x,y));
v=v,size(K2);v;K3=quotient(K2,ideal(x,y));
v=v,size(K3);v;K4=quotient(K3,ideal(x,y));
size(std(reduce(K3,std(K2))));
size(std(reduce(K4,std(K3))));

=>

0,21
0,21,22
0,21,22,21
0,21,22,21,21
18
0

K0, K1, ... と増大列になり、K4=0 (mod K3) となった時点で増大が止まることから、K3 で「消去」が終了した。結果は以下の通りで、8個の条件の式が得られた。

> ideal L=std(reduce(K3,ideal(x,y)));
// ** _ is no standard basis
> L;
L[1]=a2*b1*c0-a1*b2*c0-a2*b0*c1+a0*b2*c1+a1*b0*c2-a0*b1*c2
L[2]=b2^2*c0^2-b1*b2*c0*c1+b0*b2*c1^2+b1^2*c0*c2-2*b0*b2*c0*c2-b0*b1*c1*c2+b0^2*c2^2
L[3]=a2*b2*c0^2-a1*b2*c0*c1+a0*b2*c1^2-a2*b0*c0*c2+a1*b1*c0*c2-a0*b2*c0*c2-a0*b1*c1*c2+a0*b0*c2^2
L[4]=a2^2*c0^2-a1*a2*c0*c1+a0*a2*c1^2+a1^2*c0*c2-2*a0*a2*c0*c2-a0*a1*c1*c2+a0^2*c2^2
L[5]=a2*b0*b2*c0-a0*b2^2*c0-a1*b0*b2*c1+a0*b1*b2*c1-a2*b0^2*c2+a1*b0*b1*c2-a0*b1^2*c2+a0*b0*b2*c2
L[6]=a2^2*b0*c0-a0*a2*b2*c0-a1*a2*b0*c1+a0*a2*b1*c1+a1^2*b0*c2-a0*a2*b0*c2-a0*a1*b1*c2+a0^2*b2*c2
L[7]=a2^2*b0^2-a1*a2*b0*b1+a0*a2*b1^2+a1^2*b0*b2-2*a0*a2*b0*b2-a0*a1*b1*b2+a0^2*b2^2
L[8]=a1*b0*b2^2*c0-a0*b1*b2^2*c0+a2*b0^2*b2*c1-a1*b0*b1*b2*c1+a0*b1^2*b2*c1-a0*b0*b2^2*c1-a2*b0^2*b1*c2+a1*b0*b1^2*c2-a0*b1^3*c2-a1*b0^2*b2*c2+2*a0*b0*b1*b2*c2

コメント?

ここで、項順序の違いが明確に計算時間に影響するような気の利いた例を挙げたいところだが、それは今後の課題とする。なお、どの項順序だと計算が早くなるかは、まだよくわかっていない。上述の例は少なくとも理論的にはいい順序といえると思うが、計算の早さはまた別の話になる可能性もある。quotient()も計算に非常に時間がかかることがよくある。

*1:理論上、有限時間で終わるかはともかく

*2:環の商ではない

2025-03-07

maxima で計算中にメモリが足りなくなる問題

maxima で割と膨大な記憶容量が必要な計算をしていると、"Heap exhausted, game over" というエラーメッセージが出る場合がある*1。参考：
stackoverflow.com
メッセージの通り、メモリが足りなくなったのでこれ以上計算を続けられないということである。3つの解決策を以下に提示する。

(1) プログラムの見直し

上記リンク先にもあるように、そもそも記憶領域を無駄遣いしているかどうかを確認するのは当然重要であろう。

(2) メモリの割り当て量を増やす

設定を変更して maxima に割り当てられているメモリの量を増やす。上記リンク先で引用されている記事も参照：
http://arcoramexico.net/wxmaxima-16.12.2/info/wxmaxima.html#My-sbcl_002dbased-maxima-runs-out-of-memory
wxMaxima を使う場合は、設定の「追加パラメーター」の欄に、指定された書式（比較的新しいバージョンなら説明が書いてある）により割り当て容量を設定できる。自分の環境では、追加パラメーター

-X "--dynamic-space-size 4000"

により4000MBと設定した場合は以下の結果を得た。

:lisp (sb-ext:dynamic-space-size)
=>
4194304000

(3) ガベージコレクション

garbage_collect()によりガベージコレクションができないか試みてくれる。マニュアルによると、maxima は自動的にメモリの解放を行ってくれるようだが、自分の経験ではこのコマンドを使うことで状況が改善したこともあった。また、kill()により、もう使わない大容量のオブジェクトを消すことも検討した方がいいかもしれない。

*1:全く詳しくないので、maxima や lisp や scbl の区別というか切り分けはしていない。

2025-01-25

emacs 便利なキーバインド

メモ

複数行をN個のスペースでインデントする（選択してから）：
C-u N C-x Tab
indentation - Indent several lines in Emacs - Stack Overflow
ここで、C-u N は N 回繰り返すコマンドで、C-x Tab は選択範囲をインデントするコマンドで、この2つの合成をしている。ただし、C-x Tab 単独で実行すると、矢印キーでインデント幅を設定することになる。

www.youtube.com

2024-08-05

maxima でグラフを描く・関数の零点（未完成）

対数
シグモイド関数
シグモイド関数と対数の組み合わせ
x 軸に接している場合
2変数の場合

関数のグラフを描くときに、零点を知りたいときがある。ところが、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:12*x^3-12*x^2-3*x+4;
plot2d([%,my_log(%,x)],[x,-2,2],[y,-10,10]);

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

	poly_discriminant(f,x);
	diff(f,x);
	solve(%,[x]);
	map(lambda([X],ev(f,X)),%),factor;
	float(%);
=>
	-864
	36*x^2-24*x-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(b*log((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(b*log((-a))/log(10))
);
plot2d([f,my_log_p(f,x),my_log_m(f,x)],[x,-2,2],[y,-10,10]);

結果：

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

シグモイド関数

我々の問題意識としては、関数の絶対値が大きい範囲は注目していないので、シグモイド関数を使うのもひとつの考えだと思われる。ここでは $\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]);

対比として、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]);

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

上記の $\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(2*b+b*c/(1+abs(c)))
);

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

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

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

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

2変数の場合

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

h: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;

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

tt:1;
mag:10^0;
draw(
    dimensions=[1200,500],
    gr2d(
        title="h",
        ip_grid=[200,200],
        implicit(mag*h,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
);

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

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

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

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

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(2*b+b*c/(1+abs(c)))
);

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

現時点では違いがよくわからないのだが、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) $

パラメータ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
 + (545*sqrt(-3)-612)/(19683) * x^6;

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

*1:この式は後述の2変数関数を y 軸に制限して得られる

*2:「表れない」といっても、implicit() を使ったグラフにもその断片のようなものが観察される

2024-08-03

maxima でグラフを描く・draw2d 編

基本
計算例：三つ葉グラフ
補足：消去について
極座標表示
点オブジェクト
点からの射影
はさみグラフ
補足：零の閾値

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

複数のグラフを描画したいときは、単に順番に引数にしていけばよい。ここで、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次式が求める定義方程式である。この式は有名な式であって、以下の表記を使うことも多い。
$(x^2+y^2)^2+y^3-3 x^2 y$
この曲線は原点で3つの接線 $y=0$ , $y=\pm \sqrt{3} x$ を持っている*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)
);

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

緑のグラフは「はさみグラフ」と（勝手に）名付けることにした。注意：実はこのグラフに「点」を付け加える必要がある（後述）。

補足：消去について

上記計算でなぜ余分な因子が出てくるかというと、eliminate()は終結式を使って変数を1つずつ消去していくが、終結式（が 0 になること）はあくまでも必要条件でしかない。例えば
$ax-1=bx-1=0$
という $x$ についての連立方程式を考えると、これが（共通）解をもつための必要十分条件は
$a-b=0, ~ a \neq 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,%);

式としては
$-\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
);

右のグラフは極座標表示である。このように、複数のグラフを並べて表示することもできる。その場合はdraw()を使う。引数は個別のグラフ（scene オブジェクトとよばれる）やオプションである。この例ではcolumns=2としているが、つまり1行に2個並べるという設定である。個別のグラフはdraw2d()の代わりにgr2d()により生成する。つまりdraw2d()とはdraw()とgr2d()を組み合わせたものである。xrange,yrangeは描画の範囲、dimensionsは出力結果（全体）の大きさを指定している。元の陰関数による描画（灰色）と比べてみると、原点以外ではほぼ同じであることがわかる。また、上述の改良版の出力結果とは概ね一致している。ついでだが、描画の範囲を変更することで原点周辺を拡大した場合は次のようになる。

上述の赤いグラフの拡大：

補足：このようなグラフを描画するときに、「点」としてグラフの断片が出てくる可能性を検討しなければならない。例えば、
$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/(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のグラフは完成したことになる。

点オブジェクト

上述のように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
);

なおflatten()によりリストのネスト構造を「平坦化」しているが、しなくても自動的に適切に解釈してくれる仕組みになっているようだ。

点からの射影

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

なおこの代数曲線の特異点は原点のみである。

補足：零の閾値

次のような関数を使うことで、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:この曲線の原点での特異性を $D_4$ 型曲線特異点という

*2:代数閉体上ならヒルベルト零点定理が成立する

*3:番号は環境依存かもしれない

2024-07-16

maxima 小ネタ集（進化するコンテンツ）

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

lisp のリストとしての式
代数的数の分母の有理化
行列の生成
三角関数の式の単純化
リスト操作：apply() と map()
kill()
経過時間表示
サブリスト
現在時刻
定義した関数の内容を確認

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+2*x*y+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+2*x*y+x^2
(%o2)	abs(y+x)
(%o3)	abs
(%o4)	y+x

代数的数の分母の有理化

フラグ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の場合に限っては単にフラグの名前だけでもよい。

行列の生成

成分が何らかの法則で決まっているのような行列は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]

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

三角関数の式の単純化

trigsimp()により三角関数の式を単純化できる。また、trigexpand()により加法定理・倍角公式等を適用して式を展開できる。

trigsimp( 2*cos(t)^2+sin(t)^2 );
=>
cos(t)^2+1

trigexpand( cos(2*t)+sin(t)^2 );
=>
cos(t)^2

参考：
Maxima/三角関数・双曲線関数 - Wikibooks
maximaにおける三角関数 - 備忘録

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

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

(%i1)	apply("+",makelist(i*x^i,i,0,10));
(%o1)	10*x^10+9*x^9+8*x^8+7*x^7+6*x^6+5*x^5+4*x^4+3*x^3+2*x^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行にまとめられる。

kill()

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

経過時間表示

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

サブリスト

sublist()はリストから条件によりサブリストを作る。また、sublist_indices()はサブリストの各要素に対応したインデックスのリストを作る。

	sublist(makelist(i,i,100),primep);
	sublist_indices(%,lambda([x],x=47));
(%o1)	[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97]
(%o2)	[15]

なおprimepにより素数かどうか判定している。

現在時刻

timedate()により現在の日時を文字列として取得する。計算結果をテキストファイルとして出力するときなどに有用であると思われる。他にもいくつか時間関係の関数がある（マニュアル参照）。

定義した関数の内容を確認

dispfun()に（ユーザー定義の）関数名を入れると（複数も可）、関数の定義を返す。