使 Maxima 代码运行得更快收敛半径

问题描述

我试图找到与应用于 Mandelbrot 集的牛顿方法相关的某些泰勒级数的收敛半径。为此，我编写了一些 Maxima 代码；计算出的 R 至少对于第一个注释案例（1 / inf 而不是 1 / 2）是不正确的，所以我试图打印一些值以查看它在数字上的表现。然而，Maxima 花费了太多时间（第一个值 1 秒，第二个值 17 秒，我放弃等待第三个值）。在第三个评论案例中的表现甚至更慢。我怎样才能加快速度？

f(c,z) := z^2 + c;
/* h(c) := f(c,f(c,0)); /* simpler case is much faster */ */
h(c) := ratsimp(f(c,0))) / f(c,0)); /* polynomial with roots at period 3 components */
/* h(c) := ratsimp(f(c,0)))) / f(c,0))); /* next case is much slower */ */
hh(c) := ratsimp(diff(h(c),c));
N(z) := ratsimp(z - h(z) / hh(z)); /* Newton's method */
M(z) := ratsimp(subst(c = c + z,N(c)) - c); /* perturbed */
for ct in allroots(h(c) = 0) do
(
  print(ct),/* coefficient of Taylor series */
  a(n) := at(diff((subst(lhs(ct) = bfloat(rhs(ct)),M(z))),z,n) / factorial(n),z = 0),R : 1 / limit(abs(a(n+1)/a(n)),n,infinity),/* radius of convergence */
  print(R),/* prints 1/inf,but this is incorrect at least for the first commented case */
  print(bfloat(1 / abs(a(11)/a(10)))),/* prints after 1 second */
  print(bfloat(1 / abs(a(101)/a(100)))),/* prints after 17 seconds */
  print(bfloat(1 / abs(a(1001)/a(1000)))) /* gave up waiting */
);

EDIT 纠正 infinity 与 inf 的混淆，考虑第一个注释案例 h(c) = f(c,0)) = c^2 + c。然后 M(z) = (z^2 - h(c)) / (2 * z + 2 * c + 1) 为根 c = 0 给出 M(z) = z^2 / (1 - (-2 * z))，它是几何级数 z^2 (1 + (-2) * z + (-2)^2 * z^2 + ... + (-2)^n z^n + ...) 的总和。现在系数 a(n) = (-2)^n 所以 1 / limit(abs(a(n + 1)/a(n)),inf) 应该是 1/2，但 Maxima 报告除以零：

f(c,z) := z^2 + c;
h(c) := f(c,0));
hh(c) := ratsimp(diff(h(c),c));
N(z) := ratsimp(z - h(z) / hh(z));
M(z) := ratsimp(subst(c = c + z,N(c)) - c);
for ct in solve(h(c) = 0) do
(
  print(ct),a(n) := at(diff((subst(lhs(ct) = bfloat(rhs(ct)),invR : limit(abs(a(n+1)/a(n)),inf),print(abs(a(10+1)/a(10))),print(abs(a(100+1)/a(100))),print(abs(a(1000+1)/a(1000))),print(invR)
);

哪个打印

...
c = - 1 
`rat' replaced -1.0B0 by -1/1 = -1.0B0
`rat' replaced -1.0B0 by -1/1 = -1.0B0
2.0b0 
2.0b0 
2.0b0 
0 
c = 0 
`rat' replaced 1.0B0 by 1/1 = 1.0B0
`rat' replaced 1.0B0 by 1/1 = 1.0B0
2.0b0 
2.0b0 
2.0b0 
0 

我想以数字方式评估其他示例的高阶项以查看其行为，因为似乎没有正确计算限制。但是我的 Maxima 代码太慢而无法实现这一点，虽然 SymPy 在周期 3（三次）情况下运行正常（请参阅链接问题），但对于周期 4（度数 6）它也需要太长时间。

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）