问题描述
假设我在 a
和 b
之间定义了一个函数 f(x)
.这个函数可以有很多零,也可以有很多渐近线.我需要检索这个函数的所有零.最好的方法是什么?
Suppose I have a function f(x)
defined between a
and b
. This function can have many zeros, but also many asymptotes. I need to retrieve all the zeros of this function. What is the best way to do it?
实际上,我的策略如下:
Actually, my strategy is the following:
- 我在给定的点数上评估我的函数
- 我检测是否有符号变化
- 我找到了正在改变符号的点之间的零
我验证找到的零是否真的是零,或者这是否是渐近线
- I evaluate my function on a given number of points
- I detect whether there is a change of sign
- I find the zero between the points that are changing sign
I verify if the zero found is really a zero, or if this is an asymptote
U = numpy.linspace(a, b, 100) # evaluate function at 100 different points
c = f(U)
s = numpy.sign(c)
for i in range(100-1):
if s[i] + s[i+1] == 0: # oposite signs
u = scipy.optimize.brentq(f, U[i], U[i+1])
z = f(u)
if numpy.isnan(z) or abs(z) > 1e-3:
continue
print('found zero at {}'.format(u))
这个算法似乎有效,但我看到了两个潜在的问题:
This algorithm seems to work, except I see two potential problems:
- 它不会检测不穿过 x 轴的零(例如,在像
f(x) = x**2
这样的函数中)但是,我不认为它我正在评估的函数可能会发生这种情况. - 如果离散化点太远,它们之间可能会多于一个零,并且算法可能无法找到它们.
- It will not detect a zero that doesn't cross the x axis (for example, in a function like
f(x) = x**2
) However, I don't think it can occur with the function I'm evaluating. - If the discretization points are too far, there could be more that one zero between them, and the algorithm could fail finding them.
您是否有更好的策略(仍然有效)来查找函数的所有零点?
Do you have a better strategy (still efficient) to find all the zeros of a function?
我认为这对这个问题并不重要,但对于那些好奇的人来说,我正在处理光纤中波传播的特征方程.函数看起来像(其中V
和ell
是之前定义的,ell
是一个正整数):
I don't think it's important for the question, but for those who are curious, I'm dealing with characteristic equations of wave propagation in optical fiber. The function looks like (where V
and ell
are previously defined, and ell
is an positive integer):
def f(u):
w = numpy.sqrt(V**2 - u**2)
jl = scipy.special.jn(ell, u)
jl1 = scipy.special.jnjn(ell-1, u)
kl = scipy.special.jnkn(ell, w)
kl1 = scipy.special.jnkn(ell-1, w)
return jl / (u*jl1) + kl / (w*kl1)
推荐答案
我看到的主要问题是你是否真的可以找到所有根——正如评论中已经提到的那样,这并不总是可能的.如果你确定你的函数不是完全病态的(sin(1/x)
已经提到过),下一个是你对丢失一个或几个根的容忍度.换句话说,它是关于你准备走多长时间以确保你没有错过任何一个——据我所知,没有通用的方法来为你隔离所有的根,所以你必须自己做.你展示的已经是合理的第一步.一些评论:
The main problem I see with this is if you can actually find all roots --- as have already been mentioned in comments, this is not always possible. If you are sure that your function is not completely pathological (sin(1/x)
was already mentioned), the next one is what's your tolerance to missing a root or several of them. Put differently, it's about to what length you are prepared to go to make sure you did not miss any --- to the best of my knowledge, there is no general method to isolate all the roots for you, so you'll have to do it yourself. What you show is a reasonable first step already. A couple of comments:
- Brent 的方法在这里确实是一个不错的选择.
- 首先,处理分歧.由于在您的函数中分母中有 Bessels,您可以先求解 它们的根 - 最好在例如 Abramovitch 和 Stegun (数学世界链接).这将比使用您正在使用的临时网格更好.
- 你可以做什么,一旦你找到了两个根或分歧,
x_1
和x_2
,在区间[x_1+epsilon> 中再次运行搜索, x_2-epsilon]
.继续直到找不到更多的根(布伦特的方法保证收敛到 a 根,前提是有一个). - 如果你不能列举所有的分歧,你可能需要在验证候选人确实是分歧的时候更加小心:给定
x
不要只检查f(x)
很大,检查一下,例如|f(x-epsilon/2)|>|f(x-epsilon)| 用于epsilon
的几个值(1e-8、1e-9、1e-10 之类的). - 如果你想确保你没有简单地接触零的根,寻找函数的极值,对于每个极值,
x_e
,检查f 的值(x_e)
.
- Brent's method is indeed a good choice here.
- First of all, deal with the divergencies. Since in your function you have Bessels in the denominators, you can first solve for their roots -- better look them up in e.g., Abramovitch and Stegun (Mathworld link). This will be a better than using an ad hoc grid you're using.
- What you can do, once you've found two roots or divergencies,
x_1
andx_2
, run the search again in the interval[x_1+epsilon, x_2-epsilon]
. Continue until no more roots are found (Brent's method is guaranteed to converge to a root, provided there is one). - If you cannot enumerate all the divergencies, you might want to be a little more careful in verifying a candidate is indeed a divergency: given
x
don't just check thatf(x)
is large, check that, e.g.|f(x-epsilon/2)| > |f(x-epsilon)|
for several values ofepsilon
(1e-8, 1e-9, 1e-10, something like that). - If you want to make sure you don't have roots which simply touch zero, look for the extrema of the function, and for each extremum,
x_e
, check the value off(x_e)
.
这篇关于如何使用 numpy(和 scipy)查找函数的所有零?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!