我正在尝试使用 mpmath.findroot 的多维牛顿法以高精度数值求解方程组。这是一个示例系统:
def f(x_0, x_1, x_2, x_3, x_4, x_5, y_0, y_1, y_2, y_3, y_4, y_5, l_0, l_1,l_2, l_3, l_4, l_5, l_6, l_7, l_8, l_9, l_10, l_11, l_12, l_13, l_14):
return -x_0, -y_0 + 1, -x_5, -y_5, -x_1 - x_2, -y_1 + y_2, -x_3 - x_4, -y_3 + y_4, -(x_1 - x_5)^2 - (y_1 - y_5)^2 + 1, -(x_1 - x_4)^2 - (y_1 - y_4)^2 + 1, -(x_0 - x_5)^2 - (y_0 - y_5)^2 + 1, -(x_3 - x_4)^2 - (y_3 - y_4)^2 + 1, -(x_2 - x_3)^2 - (y_2 - y_3)^2 + 1, -(x_2 - x_5)^2 - (y_2 - y_5)^2 + 1, -(x_0 - x_5)^2 - (y_0 - y_5)^2 + 1, 2*l_10*(x_0 - x_5) + 2*l_14*(x_0 - x_5) + l_0 - .5*y_1 + .5*y_2, 2*l_9*(x_1 - x_4) + 2*l_8*(x_1 - x_5) + l_4 + .5*y_0 - .5*y_3, 2*l_12*(x_2 - x_3) + 2*l_13*(x_2 - x_5) - .5*y_0 + .5*y_4, -2*l_12*(x_2 - x_3) + 2*l_11*(x_3 - x_4) + l_6 + .5*y_1 - .5*y_5, -2*l_9*(x_1 - x_4) - 2*l_11*(x_3 - x_4) - .5*y_2 + .5*y_5, -2*l_10*(x_0 - x_5) - 2*l_14*(x_0 - x_5) - 2*l_8*(x_1 - x_5) - 2*l_13*(x_2 - x_5) + l_2 + .5*y_3 - .5*y_4, 2*l_10*(y_0 - y_5) + 2*l_14*(y_0 - y_5) + l_1 + .5*x_1 - .5*x_2, 2*l_9*(y_1 - y_4) + 2*l_8*(y_1 - y_5) + l_5 - .5*x_0 + .5*x_3, 2*l_12*(y_2 - y_3) + 2*l_13*(y_2 - y_5) + .5*x_0 - .5*x_4, -2*l_12*(y_2 - y_3) + 2*l_11*(y_3 - y_4) + l_7 - .5*x_1 + .5*x_5, -2*l_9*(y_1 - y_4) - 2*l_11*(y_3 - y_4) + .5*x_2 - .5*x_5, -2*l_10*(y_0 - y_5) - 2*l_14*(y_0 - y_5) - 2*l_8*(y_1 - y_5) - 2*l_13*(y_2 - y_5) + l_3 - .5*x_3 + .5*x_4
import sage.libs.mpmath.all as mpmath
startsol=[0, -0.345847373297000, 0.345847770952000, -0.499999005131000, 0.500000393924000, 0, 0.999999726908000, 0.938291180327000, 0.938290404618000, 0.404864576964000, 0.404866700310000, 0, 0, 0.0205563218420000, 0.00287356421947000, -0.0200611112418000, 0.00371301649518000, 0.00363938213640000, -0.00658839856658000, -0.00414258257348000, 0.0391290458552000, 0.162094656373000, 0.0813226022151000, 0.0974643704937000, 0.158202488927000, 0.0432804352887000, 0.0813226022151000]
print f(*startsol)
mpmath.mp.dps=100
ans=mpmath.findroot(f, startsol)
不幸的是,它没有给我解决方案,但会引发错误
ZeroDivisionError:矩阵在数值上是奇异的
发生这种情况是因为 findroot 正在尝试计算雅可比吗?这是否意味着系统未确定?
我使用 scipy 的 fmin_cg 找到了起点,但我想将解决方案打磨到更高的精度。作为 fmin_cg 的函数,我最小化了函数 f 的 27 个条目的平方和。
如果 mpmath.findroot 的问题无法避免,有没有更好的方法来解决这个系统的高精度问题?