我正在尝试用多维牛顿方法数值求解一组方程组,以获得高精度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: matrix is numerically singular
这是因为findroot试图计算雅可比?这是否意味着这个系统是不确定的?在
我用scipy的fmin_cg找到了起点,但我想把解决方案打磨到更高的精度。作为fmin_cg的函数,我将函数f的27个条目的平方和最小化
如果mpmath.findroot公司无法避免,有没有更好的方法来解决这个系统的高精度?在
是的,系统是欠定的,它的雅可比矩阵的秩是24(最多)而不是27。您可以使用SymPy检查:
仔细观察雅可比矩阵可以看出冗余方程是什么。例如,这些方程式编号为0-10:
^{pr2}$显然,最后一个是前四个结果的结果。这种冗余需要消除。在
使用符号Jacobian
如果mpmath.findroot公司没有给定一个计算雅可比的函数,它将使用数值微分,引入额外的误差。如果可用,最好提供一个雅可比函数。下面是一个完整性的小代码示例。在
相关问题 更多 >
编程相关推荐