Maxima解不了一個簡單方程？

01-08

我的代碼如下：
solve([x^2+y^2=1, (x-2)^2+(y-3)^2=9],[x,y]);
返回的是一個[]，難道是我的語法不對嗎？或者這是個bug？
補充：

語法應該沒錯，可以參考下這裡（譯言網 | 10分鐘教程：用Maxima解決數學問題）的例子，maxima中的賦值是「:」不是「=」

應該是bug了。

下面這個方程長得和你給的差不多，卻能給出正確的解答（可能是因為它的解是有理數）：

(%i1) solve([x^2+y^2=1,(x-2)^2+(y-3)^2=8],[x,y]); 12 5 (%o1) [[x = --, y = --], [x = 0, y = 1]] 13 13

下面這個方程和你給的有同樣的解，卻也能給出正確的解答：

(%i2) solve([x^2+y^2=1,(x-2)^2+(y-3)^2-x^2-y^2=8],[x,y]); 5/2 5/2 5/2 5/2 3 - 10 2 3 + 45 3 + 10 2 3 - 45 (%o2) [[x = - ---------, y = -----------], [x = ---------, y = - -----------]] 26 78 26 78

你給的這個方程卻偏偏解不出來：

(%i3) solve([x^2+y^2=1,(x-2)^2+(y-3)^2=9],[x,y]); (%o3) []

和它長得差不多，解也不是有理數的方程也解不出來：

(%i4) solve([x^2+y^2=1,(x-2)^2+(y-3)^2=7],[x,y]); (%o4) []

===============================

Mathematica 能給出正確的解答：

In[1]:= Solve[{x^2 + y^2 == 1, (x - 2)^2 + (y - 3)^2 == 9}, {x, y}]


               10 - 9 Sqrt[3]       3 (5 + 2 Sqrt[3])

Out[1]= {{x -&> --------------, y -&> -----------------},

                     26                    26

10 + 9 Sqrt[3] 3 (5 - 2 Sqrt[3]) &> {x -&> --------------, y -&> -----------------}} 26 26

SymPy 也能給出正確的解答：

&>&>&> from sympy import init_printing, symbols, solve, Eq &>&>&> init_printing() &>&>&> x, y = symbols("x y") &>&>&> solve([Eq(x**2 + y**2, 1), Eq((x - 2)**2 + (y - 3)**2, 9)], [x, y]) ?? ___ ___ ? ? ___ ___ ?? ??5 9?╲╱ 3 3?╲╱ 3 15? ? 9?╲╱ 3 5 3?╲╱ 3 15?? ??── + ───────, - ─────── + ──?, ?- ─────── + ──, ─────── + ──?? ??13 26 13 26? ? 26 13 13 26??

===============================

詳細分析：

在 maxima 中輸入

? solve

可以給出關於 solve 的幫助，幫助中介紹了 solve 的演算法。相關段落如下：

`solve ([&, ..., &], [&, ..., &])" solves a
system of simultaneous (linear or non-linear) polynomial equations
by calling `linsolve" or `algsys" and returns a list of the
solution lists in the variables. In the case of `linsolve" this

list would contain a single list of solutions. It takes two lists
as arguments. The first list represents the equations to be
solved; the second list is a list of the unknowns to be
determined. If the total number of variables in the equations is
equal to the number of equations, the second argument-list may be
omitted.

按這裡的說法，它實際上是調用了 algsys 函數。於是我試了

(%i6) algsys([x^2+y^2-1,(x-2)^2+(y-3)^2-9],[x,y]); (%o6) []

果然也是無解。再試試

(%i7) algsys([x^2+y^2-1,(x-2)^2+(y-3)^2-8],[x,y]); 12 5 (%o7) [[x = --, y = --], [x = 0, y = 1]] 13 13

果然給出了一對有理數解。

於是我再看 algsys 的幫助。裡邊也介紹了它的演算法。於是我決定跟著它所說的演算法試一試。

（以下的引文均出自 algsys 的幫助）

1. First the equations are factored and split into subsystems.

第一步是因式分解。這兩個方程都沒法因式分解，跳過。

2. For each subsystem &, an equation & and a variable &
are selected. The variable is chosen to have lowest nonzero
degree. Then the resultant of & and & with respect to
& is computed for each of the remaining equations & in
the subsystem &. This yields a new subsystem & in
one fewer variables, as & has been eliminated. The process
now returns to (1).

第二步用到了結式。這是解高次多項式方程組的標準的方法。於是我試了：

(%i9) resultant(x^2+y^2-1,(x-2)^2+(y-3)^2-9,x); 2 (%o9) 52 y - 60 y + 9

結果很正常。這時就只剩下 y 一個變數了，而且也只有這麼一個方程。

然後是第三步（原文比較長，我省去了無關的部分）：

3. Eventually, a subsystem consisting of a single equation is
obtained. …………
If the equation is univariate and is either linear,
quadratic, or biquadratic, then again `solve" is called if no
approximations have been introduced. …………

這裡好長，不過也不是全部都與我們的問題有關。現在我們只剩下 y 一個變數，於是來到了這裡所說的 "univariate" 的情形。而且這時的方程（ $52y^2-60y+9$ ）是二次（quadratic）的。它說要調用 solve 函數。於是我試了：

(%i10) solve(%,y); 3/2 3/2 2 3 - 15 2 3 + 15 (%o10) [y = - -----------, y = -----------] 26 26

它給出了正確的解！

再來看第四步：

4. Finally, the solutions obtained in step (3) are substituted
into previous levels and the solution process returns to (1).

這裡沒說清楚到底把第三步的答案代入到哪個方程。於是我把 y 的第一個解代入原來的第一個方程試試：

(%i11) subst(%[1],[x^2+y^2-1,(x-2)^2+(y-3)^2-9]); 3/2 2 3/2 2 (2 3 - 15) 2 2 3 - 15 2 (%o11) [x + -------------- - 1, (x - 2) + (- ----------- - 3) - 9] 676 26

然後回到第一步。

(%i12) factor(%); 2 5/2 2 7/2 676 x - 20 3 - 343 676 x - 2704 x + 28 3 + 697 (%o12) [----------------------, -------------------------------] 676 676

還是沒法因式分解。這裡調用 factor 函數的結果只是整理了一下而已。

但要注意的是，如果 x 有有理解的話，這裡是可以分解因式的！

然後又是第二步：

(%i13) resultant(%[1],%[2],x); (%o13) 0

呃……結式為0（本來也應該為0），沒辦法繼續下去了。

於是就……解不出來了……

但如果是 $x^2+y^2=1,(x-2)^2+(y-3)^2=8$ 那個方程，由於有有理解，在回到第一步的時候，是可以因式分解出一次的式子，因此就還能解出來。至於 $x^2+y^2=1,(x-2)^2+(y-3)^2-x^2-y^2=8$ 這個方程，在回到第一步的時候，雖然也沒法因式分解，但本來就有一個一次的式子，因此也能解出來。

於是我猜想，這個bug就出在這裡。
不過也只是猜想而已，不知道對不對。

為什麼你們不報bug。。。都三年了。。最新版里還是有這個bug。。。。你們不報我來報好了Maxima -- GPL CAS based on DOE-MACSYMA / Bugs / #3208 Can"t solve this simple equation..

2016年9月19日更新:

官方覺得應該是algsys函數的問題

I have started to explore this. It interesting to compare

algsys([x^2+y^2=1,(x-2)^2+(y-3)^2=9],[x,y]);

[]

and

algsys([x^2+y^2=1, (x-2)^2+(y-3)^2=8],[x,y]);

[[x = 12/13,y = 5/13],[x = 0,y = 1]]

Tracing the higher level functions in algsys

?trace(?algsys);

?trace(?algsys0);

?trace(?algsys1);

?trace(?condensesolnl);

?trace(?distrep);

?trace(?lofactors);

?trace(?findleastvar);

?trace(?bakalevel);

?trace(?presultant);

I can see that what is attempted is equivalent to (but not using sqrtdenest)

(%i4) load(sqdnst); (%o4) "/usr/share/maxima/5.38.0/share/simplification/sqdnst.mac" (%i5) eq1:y^2+x^2-1; (%o5) y^2+x^2-1 (%i6) eq2:y^2-6*y+x^2-4*x+4; (%o6) y^2-6*y+x^2-4*x+4 (%i7) r:resultant(eq1,eq2,y); (%o7) 52*x^2-40*x-11 (%i8) X:solve(r,x); (%o8) [x = -(3^(5/2)-10)/26,x = (3^(5/2)+10)/26] (%i9) Y1:solve(subst(X[1],eq2),y); (%o9) [y = -(sqrt(4077-28*3^(7/2))-78)/26, y = (sqrt(4077-28*3^(7/2))+78)/26] (%i10) Y2:solve(subst(X[2],eq2),y) ; (%o10) [y = -(sqrt(28*3^(7/2)+4077)-78)/26, y = (sqrt(28*3^(7/2)+4077)+78)/26] (%i11) subst([X[1],Y1[1]],[eq1,eq2]) $ (%i12) expand(sqrtdenest(%)); (%o12) [0,0] (%i13) subst([X[1],Y1[2]],[eq1,eq2]) $ (%i14) expand(sqrtdenest(%)); (%o14) [378/13-(4*3^(5/2))/13,0] (%i15) subst([X[2],Y2[1]],[eq1,eq2]) $ (%i16) expand(sqrtdenest(%)); (%o16) [0,0] (%i17) subst([X[2],Y2[2]],[eq1,eq2]) $ (%i18) expand(sqrtdenest(%)); (%o18) [(4*3^(5/2))/13+378/13,0]

Comparing the traces in the two cases above, I think that the failure
occurs in the back-substitution of the solutions for x and y into eq1.
I"m guessing that algsys doesn"t simplify the expression to 0 and
therefore rejects the correct solutions. More digging required.

2016年10月4日更新:

問題已經解決了,下圖是官方給我發的郵件

具體的bug看這裡https://sourceforge.net/p/maxima/bugs/3208/

受到 @AlephAlpha 的啟發，我用另一種方法解決了這個問題

如圖

我來解釋一下

i1，i2是把這兩個方程存進了e1，e2中

之前試過，在maxima5.36.1中仍存在不能解此題的bug，所以可以直接用resultant函數來進行結式運算

i3為了消去x，解出y

i5消去y，解出x

結果和我筆算結果相同，與用戶AlephAlpha（我就不再艾特了）在mathematica和SymPy中的結果也相同

所以說，遇到此類問題，直接上resultant，maxima的solve函數實在不敢恭維

p.s.真正的bug在哪兒我也不知道，solve函數的源代碼超長，看了一眼就關了...

p.p.s.不過我對它的內核很感興趣，可能之後對其做修復，然後移植到我企劃的計算器中用

p.p.p.s.真正做出來也許都猴年馬月了...不要期待太多

p.p.p.p.s.對啦，這貨對同一未知數的多次方程solve函數也解不了（如x^10+x^4=100),此時應該用allroots函數或realroots函數

語法沒錯啊，我懷疑是 bug……

用 sage 給搞定了：

sage: solve([x^2+y^2==1,(x-2)^2+(y-3)^2==9],x,y) [[x == -9/26*sqrt(3) + 5/13, y == 3/13*sqrt(3) + 15/26], [x == 9/26*sqrt(3) + 5/13, y == -3/13*sqrt(3) + 15/26]]

Maxima我沒用過，我猜測你是不是用錯了賦值等號「=」以及判斷等號「==」？

Mathematica代碼這樣寫：

Solve[{x^2 + y^2 == 1, (x - 2)^2 + (y - 3)^2 == 9}, {x, y}]

完全沒有問題，得到了正確的結果

可以解出，有圖為證。

不清楚你為毛解不出，我的事5.40版本哦。

---sorry,沒有看前面BUG已修正。

16.04.02版本bug依舊

最新版16.12.0已經解決

用這個軟體的初衷是因為可以在ssh下工作,別無它用