我使用 scipy.optimize.root
和 hybr
方法(最好的方法?)来查找数字函数的根
我在每次迭代时打印残差
三角洲 d 117.960112417
三角洲 d 117.960112417
三角洲 d 117.960112417
三角洲 d 117.960048733
三角洲 d 117.960112427
三角洲 d 117.960112121
增量 d 1.46141491664
增量 d 0.0322651167588
增量 d 0.000363688881595
三角洲 d 4.05494689256e-08
如何通过增加步长来加速求根,尤其是在第一次迭代之间? 我不知道该算法究竟是如何工作的,但看起来很奇怪,前 3 个结果相同,而后 3 个结果也完全相同。
Reading the doc ,我尝试修改 eps
因子,但没有成功
编辑:@sasha,这是一个非常基本的功能来说明问题
def f(X1,X2):
print ' X1 , diff , norm ' , X1 , X2 - X1 , np.linalg.norm(X2 - X1)
return X2 - X1
Xa = np.array([1000,1000,1000,1000])
Xb = np.array([2000,2000,2000,2000])
SOL = scipy.optimize.root(f,Xa,(Xb,))
结果如下 无论 X 的长度如何,我们在开始时都有 3 次相同的迭代
X1 , diff , norm [1000 1000 1000 1000] [1000 1000 1000 1000] 2000.0
X1 , diff , norm [ 1000. 1000. 1000. 1000.] [ 1000. 1000. 1000. 1000.] 2000.0
X1 , diff , norm [ 1000. 1000. 1000. 1000.] [ 1000. 1000. 1000. 1000.] 2000.0
X1 , diff , norm [ 1000.0000149 1000. 1000. 1000. ] [ 999.9999851 1000. 1000. 1000. ] 1999.99999255
X1 , diff , norm [ 1000. 1000.0000149 1000. 1000. ] [ 1000. 999.9999851 1000. 1000. ] 1999.99999255
X1 , diff , norm [ 1000. 1000. 1000.0000149 1000. ] [ 1000. 1000. 999.9999851 1000. ] 1999.99999255
X1 , diff , norm [ 1000. 1000. 1000. 1000.0000149] [ 1000. 1000. 1000. 999.9999851] 1999.99999255
X1 , diff , norm [ 2000. 2000. 2000. 2000.] [-0. -0. -0. -0.] 4.36239133705e-09
X1 , diff , norm [ 2000. 2000. 2000. 2000.] [ 0. 0. 0. 0.] 0.0
最佳答案
首先,我认为您将迭代与函数调用混淆了,这两者并不完全相同。因为您没有为求解器提供 Jacobean 函数,它必须估计 Jacobean(或者可能只是它的一部分)本身。 Jacobean 基本上是导数的多维等价物。它表示当您稍微改变输入时目标函数的输出如何变化。
大多数数值求解器通过在非常接近当前猜测的某个点评估目标函数并检查输出变化多少来以数值方式估计 Jacobean。我的猜测是您看到的前几个调用是评估目标函数,然后估计 Jacobean。您看到任何实际变化的第一个调用发生在它估计 Jacobean 之后,然后用它来计算根的下一个猜测。
如果您想自己检查一下,请尝试为求解器提供一个回调函数。它将在每次迭代中调用此函数,您可以查看它在每次迭代中的位置。我认为您会发现它仅在几次迭代中收敛,但每次迭代都会多次调用该函数。
当然,您可以通过为求解器提供一个 Jacobean 函数来避免所有这些工作,该函数可以调用它来计算某个点的 Jacobean。如果你这样做,它就不需要多次调用来估计它。
文档包含有关如何添加回调和提供 Jacobean 函数的信息。如果需要,我可以添加示例。
http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.optimize.root.html
关于python - `scipy.optimize.root` 更快的根查找,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/37661586/