本文介绍了如何使用 numpy(和 scipy)查找函数的所有零?的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

假设我在 ab 之间定义了一个函数 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:

  1. 我在给定的点数上评估我的函数
  2. 我检测是否有符号变化
  3. 我找到了正在改变符号的点之间的零
  4. 我验证找到的零是否真的是零,或者这是否是渐近线

  1. I evaluate my function on a given number of points
  2. I detect whether there is a change of sign
  3. I find the zero between the points that are changing sign
  4. 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:

  1. 它不会检测不穿过 x 轴的零(例如,在像 f(x) = x**2 这样的函数中)但是,我不认为它我正在评估的函数可能会发生这种情况.
  2. 如果离散化点太远,它们之间可能会多于一个零,并且算法可能无法找到它们.
  1. 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.
  2. 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?

我认为这对这个问题并不重要,但对于那些好奇的人来说,我正在处理光纤中波传播的特征方程.函数看起来像(其中Vell是之前定义的,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_1x_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 and x_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 that f(x) is large, check that, e.g. |f(x-epsilon/2)| > |f(x-epsilon)| for several values of epsilon (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 of f(x_e).

这篇关于如何使用 numpy(和 scipy)查找函数的所有零?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!

08-04 07:42
查看更多