1

我正在尝试使用 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 的问题无法避免，有没有更好的方法来解决这个系统的高精度问题？

4

1 回答 1

1

是的，系统是欠定的，它的雅可比矩阵的秩为 24（最多）而不是 27。您可以使用 SymPy 进行检查：

import sympy as sp
vars = sp.var('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')
F = sp.Matrix([-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 $
J = F.jacobian(vars)
print(J.rank())

更仔细地查看雅可比矩阵可以看出冗余方程是什么。例如，这些是编号为 0-3 和 10 的方程式：

-x_0      = 0 
-y_0 + 1  = 0
-x_5      = 0
-y_5      = 0 
-(x_0 - x_5)^2 - (y_0 - y_5)^2 + 1  = 0

显然，最后一个是前四个的结果。这种冗余需要消失。

使用符号雅可比行列式

如果 mpmath.findroot 没有给出计算雅可比行列式的函数，它将使用数值微分，从而引入额外的错误。如果可用，最好提供雅可比函数。这是一个完整的小代码示例。

import sympy as sp
import mpmath as mp
vars = sp.var('x, y, z')
F = [x**2 + y**3 + z**4 - 6, x + y + z - 2, 3*x**2 + y - 5]
J = sp.Matrix(F).jacobian(vars)
f = lambda x0,y0,z0 : [Fc.subs(list(zip(vars, [x0,y0,z0]))) for Fc in F]
Jac = lambda x0,y0,z0 : J.subs(list(zip(vars, [x0,y0,z0]))).tolist()
start = [1, 1, 1]
print(mp.findroot(f, start, J = Jac))
于 2016-04-25T21:48:53.820 回答