numpy（とscipy）を使って関数のすべてのゼロを見つける方法は？

Question

f(x)がaとbの間に定義されているとします。この関数は多くのゼロを持つことができますが、多くの漸近線も持ちます。私はすべてこの関数のゼロを取得する必要があります。それを行う最善の方法は何ですか？

実は、私の戦略は以下の通りです：

1. 私は私がポイント間のゼロを見つける看板
2. の変更があるかどうかを検出ポイント
3. の与えられた数に私の機能を評価

4. たゼロが本当にゼロである場合、または私は検証に署名変化していること、これは漸近線

   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)) 
   

0であれば

このアルゴリズムが動作しているようですが、私は2つの潜在的な問題を参照除い：

1. をそれは（例えば、f(x) = x**2のような機能に）x軸と交差しないゼロでは検出されませんが、私はそれが評価している機能で起こるとは思わない。
2. 離散化ポイントがあまりにも遠すぎると、それらの間に1より多くのゼロが存在する可能性があり、アルゴリズムがそれらを見つけることができない可能性があります。

関数のすべてのゼロを見つけるには、より効率的な方法がありますか？

---

私はそれが質問のために重要だとは思いませんが、興味がある人のために、私は、光ファイバにおける電波伝搬の特性方程式を扱っています。機能は次のように（どこVとellが以前に定義されており、ellは正の整数）になります。

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) 

Answer 1

私はこれを見る主な問題は、すでにコメントで述べてきたように、これは常に可能ではない---です。あなたの機能が完全に病理学的ではないことが確かな場合（sin(1/x)は既に述べられています）、次は、ルートまたはそのいくつかを逃すことに対するあなたの寛容です。別の言い方をすると、それはあなたが何も見逃していないことを確かめるためにどのくらいの長さを覚悟しているかについてです。私の知る限りでは、あなたのためにすべての根を分離する一般的な方法はないので、自分でやってください。あなたが示していることは、合理的な第一歩です。いくつかのコメント：

• ブレントの方法は確かに良い選択です。
• まず、相違点を扱います。あなたの関数では分母にベッセルがあるので、最初にの根をとして解くことができます。例えば、AbramovitchとStegun（Mathworld link）を見てください。これは、使用しているアドホックグリッドを使用するよりも優れています。
• 2つのルーツまたはダイバージェントを見つけたら、x_1とx_2の間で、[x_1+epsilon, x_2-epsilon]の間隔で検索を再度実行します。それ以上の根が見つからなくなるまで続行する（ブレントの方法は、にのルートに収束することが保証されている。
• すべての相違点を列挙できない場合は、候補者が確かに発散性であることを確認するために少し注意したいことがあります。xを指定すると、f(x)が大きいことを確認するだけではありません。 epsilon（1e-8,1e-9,1e-10、そのようなもの）のいくつかの値については、|f(x-epsilon/2)| > |f(x-epsilon)|です。
• 単にゼロに触れる根を持たないようにするには、関数の極値を調べ、各極値の場合はx_eの値をf(x_e)にチェックします。

Answer 2

なぜあなたはnumpyに限られていますか？ scipyのダウンロードは、正確に何をしたいんパッケージがあります：数値プログラミングは難しいですので、それをしない:)

---

をとにかく、場合：私は学んだ

http://docs.scipy.org/doc/scipy/reference/optimize.nonlin.html

1レッスンあなたは自分でアルゴリズムを構築するには死んでいる、scipyのドキュメントページにリンクしている（永遠にロードする、btw）は、アルゴリズムのリストを提供する。私が以前に使った方法の1つは、問題に必要な程度に関数を離散化することです。 （つまり、\ delta xを調整して、問題の特性サイズよりもはるかに小さくなるようにします。）これにより、関数の特徴（符号の変更など）を探すことができます。そして、幼稚園以降の線分の導関数を簡単に計算することができます。したがって、離散化された関数には明確な一次導関数があります。特性のサイズよりも小さくなるようにdxを調整したので、問題に重要な機能の機能を見逃すことはありません。

「特性のサイズ」が意味するものを知りたい場合は、長さまたは1 /長さの単位で関数のパラメータを探します。すなわち、ある関数f（x）について、xは長さの単位を持ち、fは単位を持たないと仮定する。次に、xが増えていくものを探します。たとえば、cos（\ pi x）を離散化する場合、xを乗算するパラメータ（xの長さの単位がある場合）は、1 /長さの単位を持つ必要があります。したがって、cos（\ pi x）の特徴的なサイズは1/\ piです。離散化をこれよりずっと小さくすると、問題は発生しません。確かに、このトリックは常に機能するとは限りませんので、いくつかの修正を行う必要があります。あなたが実際にすべてのルーツを見つけることができるかどうか