MMA 无法计算 方程 t^5+t+3=0 的 5 个根的平方和

uk702

文 https://zhuanlan.zhihu.com/p/341394836 描述了用 maple 解决如下问题：

设 $t_1,t_2,...t_5$ 为方程 $t^5+t+3=0$ 的5个根，求 解以下问题：$\frac{1}{t_1^5}+\frac{1}{t_2^5}+...+\frac{1}{t_5^5$

我尝试用 Mathematica 来求解，代码如下（希望没错 ）：

1. ClearAll["Global`*"];
   
2. ans=Simplify@Solve[t^5+t+3==0];
   
3. Simplify@Fold[(#1)^(-5) + (#2)^(-5) &, t /. ans]
   

复制代码

不知什么原因，MMA 没有给出简化解来。

于是我又尝试求数值解：

1. ClearAll["Global`*"];
   
2. ans=Simplify@Solve[t^5+t+3==0] // N;
   
3. Fold[(#1)^(-5) + (#2)^(-5) &, t /. ans]
   

复制代码

MMA 给出答案为：-1129.86+1174.16 i
竟然和该文给的： $-\frac{406}{243}$ 完全不同！晕~~~

甚至，感觉手工算也不难的 $t_1^2+t_2^2+...+t_5^2$：

1. ClearAll["Global`*"];
   
2. ans = Simplify@Solve[t^5 + t + 3 == 0];
   
3. FullSimplify@Fold[#1^2 + #2^2 &, t /. ans]
   

复制代码

MMA 竟然算个没完没了！若改为数值解，

1. ClearAll["Global`*"];
   
2. ans = Simplify@Solve[t^5 + t + 3 == 0] // N;
   
3. FullSimplify@Fold[#1^2 + #2^2 &, t /. ans]
   

复制代码

MMA 给出的结果是 21.9442 + 8.73535 i 。

难道 MMA 学的韦达定理跟我们不一样？

uk702

暴汗！代码真写错了。FullSimplify@Fold[#1^2 + #2^2 &, t /. ans] 计算的是 \( (((t_1^2+t_2^2)^2+t_3^2)^2+t_4^2)^2+t_5^2 \) 的值。

要计算 \( \frac{1}{t_1^5}+\frac{1}{t_2^5}+...+\frac{1}{t_5^5} \) 的值，发现用这个代码 FullSimplify@Sum[Part[t^(-5) /. ans, i], {i, 5}]  可以。

这个代码也可以：
FullSimplify@Total[t^(-5) /. Solve[t^5 + t + 3 == 0]]

这是目前发现的最简单的写法了。

uk702

用 Fold 来计算也是可以的，但要改成：

1. ClearAll["Global`*"];
   
2. ans=Simplify@Solve[t^5+t+3==0];
   
3. ans = Prepend[ans, {t -> 0}]; (* 在 ans 的最前面添加一项 t->0 *)
   
4. FullSimplify@Fold[#1 + (#2)^(-5) &, t /. ans] (* 计算 0 + t1^(-5) + t2^(-5) + ... + t5^(-5) *)

复制代码

补充内容 (2022-12-27 08:18):
可进一步简化为: ans=Simplify@Solve[t^5+t+3==0]; FullSimplify@Fold[#1 + (#2)^(-5) &, 0, t /. ans]

wayne

可以这样，答案是 -(406/243)

1. RootSum[#^5 + # + 3 &, #^(-5) &]

复制代码

northwolves

Rootsum的计算机制是什么？ 即使是150次方，也是秒出

$RootSum[#^5 + # + 3 &, 1/#^150 &]$

$\frac{2715490287490597683785158637255586579628505626924395171690913387}{369988485035126972924700782451696644186473100389722973815184405301748249}$

mathe

递推数列呀
设$a_n=\sum_{h=1}^5 t_h^n$
于是容易得出$a(n+5)=-a(n+1)-3a(n)$
于是
\(\begin{pmatrix}a(n+5)\\a(n+4)\\a(n+3)\\a(n+2)\\a(n+1)\end{pmatrix}=
\begin{bmatrix}0&0&0&-1&-3\\
1&0&0&0&0\\
0&1&0&0&0\\
0&0&1&0&0\\
0&0&0&1&0
\end{bmatrix}
\begin{pmatrix}a(n+4)\\a(n+3)\\a(n+2)\\a(n+1)\\a(n)\end{pmatrix}
\)
然后显然a(0)=5,a(1)=a(2)=a(3)=0,a(-1)=-1/3
于是
\(\begin{pmatrix}a(n+4)\\a(n+3)\\a(n+2)\\a(n+1)\\a(n)\end{pmatrix}=
\begin{bmatrix}0&0&0&-1&-3\\
1&0&0&0&0\\
0&1&0&0&0\\
0&0&1&0&0\\
0&0&0&1&0
\end{bmatrix}^{n+1}
\begin{pmatrix}a(3)\\a(2)\\a(1)\\a(0)\\a(-1)\end{pmatrix}
\)
所以
\(a(n)=\begin{pmatrix}0&0&0&0&1\end{pmatrix}\begin{bmatrix}0&0&0&-1&-3\\
1&0&0&0&0\\
0&1&0&0&0\\
0&0&1&0&0\\
0&0&0&1&0
\end{bmatrix}^{n+1}
\begin{pmatrix}0\\0\\0\\5\\-\frac13\end{pmatrix}
\)


? f(n)=[0,0,0,0,1]*M^(n+1)*[0,0,0,5,-1/3]~
%8 = (n)->[0,0,0,0,1]*M^(n+1)*[0,0,0,5,-1/3]~
? f(-1)
%9 = -1/3
? f(-150)
%10 = 2715490287490597683785158637255586579628505626924395171690913387/36998848               5035126972924700782451696644186473100389722973815184405301748249
? f(-300)
%11 = 7373664376498616689947802511828188559166026210098992096965241773198299059348389353763026565134484527210885191229025722842141751/136891479058588375991326027382088315966463695625337436471480190078368997177499076593800206155688941388250484440597994042813512732765695774566001
? f(-1000)
%12 = 780329251088398005228108048464490410963582190572609858477362098807131513163365633040884888164208256386043441790172419796218003903978347630488120959441685971180697788629039509724914064414637384573490932683219209817093798763715361425846945448165113048373575427631800676781983912785784413197062729740388467014833103740249267956988092605300546062903960132481563949526049777380379845772926104380699253536402082670668240809773124/1322070819480806636890455259752144365965422032752148167664920368226828597346704899540778313850608061963909777696872582355950954582100618911865342725257953674027620225198320803878014774228964841274390400117588618041128947815623094438061566173054086674490506178125480344405547054397038895817465368254916136220830268563778582290228416398307887896918556404084898937609373242171846359938695516765018940588109060426089671438864102814350385648747165832010614366132173102768902855220001

wayne

我统计了一下计算耗时，可以确定 Mathematica的RootSum内置算法还不是 二分幂算法, 可能 比矩阵二分幂要快一些
计算$r^{-10^6}$，RootSum耗时2.10781 秒，MatrixPower耗时 11.7358秒，差了6倍
计算$r^{-10^7}$，RootSum耗时20.5403 秒，MatrixPower耗时 221.61秒，差了10倍
计算$r^{-10^8}$，RootSum耗时301.689 秒，MatrixPower耗时 4116.12秒，差了14倍

1. n = -10^7; Timing[RootSum[#^5 + # + 3 &, #^(n) &];]

复制代码

借用mathe的公式：

1. m={{0,0,0,-1,-3},{1,0,0,0,0},{0,1,0,0,0},{0,0,1,0,0},{0,0,0,1,0}};
   
2. Timing[First[MatrixPower[m,n-3].Transpose[{0,0,0,5,-1/3}]];]

复制代码

mathe

牛顿恒等式只处理了次数不超过变量数目的情况，次数大于变量数目时，本质上就递推数列了

mathe

精度合理估计的话，采用浮点计算会更快