用python找到超越方程的根答案

【问题标题】：Find root of a transcendental equation with python用python找到超越方程的根
【发布时间】：2017-08-20 05:21:09
【问题描述】：

我必须解决以下超越方程

cos(x)/x=c

对于给定的常数 c。

例如我在 Mathematica 中做了一个简短的代码，其中我为常数 c 生成了一个随机值列表

const = Table[RandomReal[{0, 5}], {i, 1, 10}]

(*{1.67826, 0.616656, 0.290878, 1.10592, 0.0645222, 0.333932, 3.59584, \
2.70337, 3.91535, 2.78268}*)

比我定义的函数

f[x_, i_] := Cos[x]/x - const[[i]]

并开始寻找根源：

Table[FindRoot[f[x, i] == 0, {x, 0.1}][[1, 2]], {i, 1, Length[const]}]
(*{0.517757, 0.947103, 1.21086, 0.694679, 1.47545, 1.16956, 0.26816, \
0.347764, 0.247615, 0.338922}*)

现在我很想在 python 中编写类似的程序（可能使用 numpy？），但我真的找不到任何好的现有答案来解决这样的问题。有人可以帮忙吗？

【问题讨论】：

使用随机样本解决这个问题是否重要？如果没有，请查看Newton's method
@Michail，不，绝对不是！我只使用随机生成器来提供一些合成数据。在实际情况下，我会有一些实验数据。
你想如何处理多个根？找到所有这些重要还是要找到接近 0 的？
@ChrisMueller 我只对第一个正根感兴趣。没有了。

标签： python numpy numeric transcendental-equation

【解决方案1】：

您可以使用sympy：

>>> from sympy import cos, Symbol, nsolve
>>> x = Symbol('x')
>>> consts = [random.random() for _ in range(10)]
>>> [nsolve(cos(x)/x - c, x, 1) for c in consts]
[mpf('0.89659506789294669'),
 mpf('0.96201114853313738'),
 mpf('0.74186728791161379'),
 mpf('1.1720944924353926'),
 mpf('0.92953351945607071'),
 mpf('0.96626530553984035'),
 mpf('1.4270719610604761'),
 mpf('0.85968954499458035'),
 mpf('0.86682911058530746'),
 mpf('0.91591678333479274')]

【讨论】：

你知道这如何处理多个根吗？例如，如果您设置 c=0.1，则该等式似乎有 ~10 个根。
nsolve 只解决一个根，在幕后它使用mpmath.findroot() 所以你也许可以直接使用那个较低级别的库。

【解决方案2】：

我过去实现这一目标的一种方法是使用scipy.optimize.minimize 来找到平方函数的最小值。

from scipy.optimize import minimize
from numpy import cos

def opt_fun(x, c):
    return (cos(x)/x - c)**2

const = 1.2
res = minimize(lambda x: opt_fun(x, const), x0=0.001)

# Check if the optimization was successful
print(res.success)
# >> True

# Extract the root from the minimization result
print(res.x[0])
# >> 0.65889256782472172

这绝不是万无一失的，但它可以快速准确。例如，如果有多个根，minimize 会从您选择的初始点找到“下坡方向”的根，这就是我在上面选择一个小的正值的原因。

另一个需要注意的问题是数量级差异很大的数字，这对于最小化问题总是如此。在您的等式中，随着c 变得非常大，第一个正根变得非常小。如果您最终试图在这种情况下找到根源，您可能需要将x 缩放到接近 1 才能获得准确的结果 (an example here)。

【讨论】：

scipy.optimize 还提供了各种root finding functions，尽管我怀疑，对于这个例子，你提出的基于优化的方法表现得更差。
@Stelios：我在下面使用了这样的求根算法。 c=1.2 示例的值确实几乎相同。
通过寻找函数平方的最小值来找到方程的根并不是真正值得推荐的，因为它在根和原始函数的零值最小值之间引入了歧义，并且最小值类似于可以确定根的精度的平方根。最好通过其他答案提到的方法直接找到根源。

【解决方案3】：

对于这种简单的单变量函数，您可以使用 Chebfun 的 python 实现轻松找到感兴趣区间内的所有根。我知道有两个，Chebpy 和 pychebfun，它们都很出色。

例如，使用 Chebpy，可以执行以下操作以在区间 [0.5, 12] 中找到 cos(x)/x - 0.05 的根：

from chebpy import chebfun

x = chebfun('x', [0.5, 12])
c = 0.05
f = np.cos(x)/x - c

rts = f.roots()
print(rts)

[ 1.4959 4.9632 7.4711 11.6152]

【讨论】：

这看起来是一个非常优雅的解决方案，但我在 Windows 上安装 ChebPy 时遇到了很多麻烦......你甚至无法想象。
你有在windows 10, python 3.6.0上安装它的经验吗？
@skrat 从cheby的安装说明来看，windows可能不支持。也许试试 pychebfun，我已经用 python 3.5 在 windows 机器上成功安装了它。

【解决方案4】：

或者，您可以使用root：

import numpy as np
from scipy.optimize import root

def func_cos(x, c):
    return np.cos(x) / x - c

crange = range(1, 11)

res = [root(func_cos, 0.5, args=(ci, )).x[0] for ci in crange]

那么res如下：

[0.73908513321516056,
 0.45018361129487355,
 0.31675082877122118,
 0.24267468064089021,
 0.19616428118784215,
 0.16441893826043114,
 0.14143076140757282,
 0.12403961812459068,
 0.11043425911223313,
 0.099505342687387879]

如果您对使用root 求解方程组感兴趣，可以查看this answer。

【讨论】：