方程式の解を指定範囲で全て探すプログラム

以下のプログラムは、方程式を関数f(x)=0で与えたときに、指定した範囲内で解を全て求めることを意図したものである。
基本は、範囲をdivの数で分割し、それぞれの領域でScipyのbrentqを適用し、解を集めている。
しかし、関数f(x)が解の前後で符号を変えないケースがあり、このような場合の解を求めるためにNewton法も用いている。

f1d_solve_range.py

import numpy as np

def f1d_solve_range(f,range,div=1000,args=(),newtontirigger=1.E-4):
    """
Solve f(x)=0 to otain all solutions between range
range is devided by div to find solution in each subrange

When |f(x)| < 1.E-4, Newton Method is also used
because brentq can not find solution if f(x) does not change sign at a solution.
"""
    from scipy.optimize import brentq,newton

    x=np.linspace(range[0],range[1],div)
    dx=x[1]-x[0] 

    #solution 1 by bysection
    solution1=[]
    for x0  in x[:-1]:
        try:
            a=brentq(f,x0,x0+dx,args=args)
            solution1.append(a)
        except:
            continue
    solution1=np.unique(solution1)

    # solution2 by Netwon method just in case f(a) f(b) has same sign
    solution2=[] 
    for x0 in x:
        if np.abs(f(x0,*args))<newtontirigger:
            b=newton(f,x0,args=args)
            solution2.append(b)
    solution2=np.unique(solution2)
    solution2.clip(range[0],range[1])  # remove out of range

    # final solution
    solution=solution1
    for s2 in solution2:
        if np.alltrue(np.abs(solution-s2)>dx):
            print "New solution added by Newtonian!:",s2
            solution=np.append(solution,s2)
    
    return np.sort(solution)  

例１　

#example 1
    def f(x):
        return x**3 
    
    print "ex1=",f1d_solve_range(f,(-10.,10))

結果:

ex1= [ 0.]

例2　

#example 2
    def f(x):
        return x**3-x 
    print "ex2=",f1d_solve_range(f,(-10.,10)) 

結果：
ex2= [-1.  0.  1.]


例3

#example 3
    def f(x):
        return x**3-x**2 
    print "ex3=",f1d_solve_range(f,(-10.,10)) 
    x=np.linspace(-2,2,100)

結果：
ex3= New solution added by Newtonian!: 1.92717530048e-08
[  1.92717530e-08   1.00000000e+00]
解の前後で符号を変えないケースでNewton法で解が求まっている。

のグラフを参照

例4

#example 4
    def f(x):
        return x**2+1 
    print "ex4=",f1d_solve_range(f,(-10.,10)) 

結果：

ex4= []

解なし。

例5

#example 5
    def f(x,a,b,c):
        return a*x**2 + b*x +c 
    print "ex5=",f1d_solve_range(f,(-10.,10),args=(1,2,1)) 

結果:
ex5= New solution added by Newtonian!: -0.999999982865
[-0.99999998]