python - 需要提高 fsolve 的准确性以找到倍数根

Question

我正在使用此代码获取非线性函数的零点。最肯定的是，该函数应该有 1 个或 3 个零

import numpy as np
import matplotlib.pylab as plt
from scipy.optimize import fsolve

[a, b, c] = [5, 10, 0]

def func(x):
    return -(x+a) + b / (1 + np.exp(-(x + c)))

x = np.linspace(-10, 10, 1000)

print(fsolve(func, [-10, 0,  10]))
plt.plot(x, func(x))
plt.show()

在这种情况下，代码可以毫无问题地给出 3 个预期的根。但是，当 c = -1.5 时，代码会遗漏一个根，而当 c = -3 时，它会找到一个不存在的根。

我想计算许多不同参数组合的根，因此手动更改种子不是一个实用的解决方案。

我感谢任何解决方案、技巧或建议。

Answer 1

您需要的是一种自动获得函数根的良好初始估计值的方法。这通常是一项艰巨的任务，然而，对于单变量、连续函数，它相当简单。这个想法是要注意 (a) 此类函数可以通过适当大阶的多项式近似为任意精度，并且 (b) 存在用于查找(所有)多项式根的有效算法。幸运的是，Numpy 提供了执行多项式逼近和求多项式根的函数。

让我们考虑一个特定的函数

[a, b, c] = [5, 10, -1.5]

def func(x):
    return -(x+a) + b / (1 + np.exp(-(x + c)))

以下代码使用polyfit和 poly1d近似func在感兴趣的范围内 ( -10<x<10 ) 通过多项式函数 f_poly订单10 .

x_range = np.linspace(-10,10,100)
y_range = func(x_range)

pfit = np.polyfit(x_range,y_range,10)

f_poly = np.poly1d(pfit)

如下图所示，f_poly确实是 func 的一个很好的近似值.通过增加阶数可以获得更高的精度。然而，在多项式逼近中追求极高的准确性是没有意义的，因为我们正在寻找根的近似估计，这些根稍后将被 fsolve 细化。

多项式近似的根可以简单地获得为

roots = np.roots(pfit)
roots

array([-10.4551+1.4893j, -10.4551-1.4893j, 11.0027+0.j , 8.6679+2.482j , 8.6679-2.482j , -5.7568+3.2928j, -5.7568-3.2928j, -4.9269+0.j , 4.7486+0.j , 2.9158+0.j ])

正如预期的那样，Numpy 返回 10 个复根。然而，我们只对区间 [-10,10] 内的实根感兴趣。 .这些可以提取如下:

x0 = roots[np.where(np.logical_and(np.logical_and(roots.imag==0, roots.real>-10), roots.real<10))].real
x0

array([-4.9269, 4.7486, 2.9158])

数组 x0可以作为fsolve的初始化:

fsolve(func, x0)

array([-4.9848, 4.5462, 2.7192])

备注:pychebfun package提供了一个函数，直接给出一个区间内函数的所有根。它也基于执行多项式逼近的想法，但是，它使用更复杂(但更有效)的方法。它会自动选择近似的最佳多项式阶数(无需用户输入)，多项式根实际上等于真实根(无需通过 fsolve 对其进行细化)。

这个简单的代码给出了与 fsolve 相同的根

import pychebfun

f_cheb = pychebfun.Chebfun.from_function(func, domain = (-10,10))
f_cheb.roots()