需要提高 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）这个 class 函数可以通过适当大阶的多项式近似为任意精度，并且（b）有有效的算法来找到（所有）多项式的根.幸运的是，Numpy 提供了执行多项式逼近和求多项式根的函数。

让我们考虑一个具体的函数

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

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

以下代码使用 polyfit 和 poly1d 通过阶数 f_poly 的多项式函数在感兴趣的范围 (-10<x<10) 上近似 func 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 returns 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包提供了一个函数，可以直接给出一个函数在一个区间内的所有根。它也基于执行多项式逼近的想法，但是，它使用更复杂（但更有效）的方法。它会自动选择近似的最佳多项式阶数（无需用户输入），多项式根实际上等于真实根（无需通过 fsolve 对其进行细化）。

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

相同的根

import pychebfun

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

Answer 2

在两个固定点之间（即 df/dx=0），您有一个或零个根。在您的情况下，可以通过分析计算两个固定点：

[-c + log(1/(b - sqrt(b*(b - 4)) - 2)) + log(2),
 -c + log(1/(b + sqrt(b*(b - 4)) - 2)) + log(2)]

因此您需要在三个区间中找到一个零。使用 Sympy 可以让您免于手动进行计算。它的 sy.nsolve() 允许在一个区间内稳健地找到一个零：

import sympy as sy

a, b, c, x = sy.symbols("a, b, c, x", real=True)

# The function:
f = -(x+a) + b / (1 + sy.exp(-(x + c)))
df = f.diff(x)  # calculate f' = df/dx
xxs = sy.solve(df, x)  # Solving for f' = 0 gives two solutions

# numerical values:
pp = {a: 5, b: 10, c: .5}  # values for a, b, c
fpp = f.subs(pp)
xxs_pp = [xpr.subs(pp).evalf() for xpr in xxs]  # numerical stationary points
xxs_pp.sort()  # in ascending order

# resulting intervals:
xx_low = [-1e9,      xxs_pp[0], xxs_pp[1]]
xx_hig = [xxs_pp[0], xxs_pp[1],       1e9]

# calculate roots for each interval:
xx0 = []
for xl_, xh_ in zip(xx_low, xx_hig):
    try:
        x0 = sy.nsolve(fpp, (xl_, xh_), solver="bisect")  # calculate zero
    except ValueError:  # no solution found
        continue
    xx0.append(x0)

print("The zeros are:")
print(xx0)

sy.plot(fpp)  # plot function

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

need to improve accuracy in fsolve to find multiples roots

python

mathematical-optimization

scipy