Scipy库中哪个函数可以求解一元函数的根？

Rainyboy

如题，我找了半天，主要参考的是下面三个文档：
NumPy Reference
SciPy Reference Guide
NumPy User Guide
感觉整个Numpy库只是提供了数值计算的基础，最突出的特点就是扩展一系列比Python内建的列表类更好用，更接近数值计算需求的向量类，用它可以方便地声明各种各样的一维或多维数组。


SciPy吧倒是提供了很多工具，但是我找来找去就是找不到对一元函数求根的函数？类似于Matlab中的fzero.

给一段MATLAB中的fzero的用法吧，各位看看Python的这些库类中有没有类似的？没有的话估计我得自己写了……

1. p0=0.2969;
   
2. p1=-0.1260;
   
3. p2=-0.3516;
   
4. p3=0.2843;
   
5. p4=-0.1015;
   
6. %翼型内联函数
   
7. y_s = inline('p0*(k1*r+0.5)^0.5 + p1*(k1*r+0.5) + p2*(k1*r+0.5)^2 + p3*(k1*r+0.5)^3 + p4*(k1*r+0.5)^4 - k2*r','r',' p0' ,' p1', 'p2' ,'p3', 'p4','k1','k2');
   
8. k1 =1;
   
9. k2 =2;
   
10. r = fzero(y_s ,[0 0.53],[ ], p0,p1,p2,p3,p4,k1,k2);

复制代码

Rainyboy

好吧……我还是找不到，自己写了一个牛顿法求根……请大家指教！
===================代码分割线======================

# -*- coding: cp936 -*-
#牛顿法求一元函数在制定初值附近的根
#范雨 2010/12/5

def FZeros(func,bound,paralist):
    """对指定的一元函数在指定的点求导数值

        func     函数对象，其定义形式应为，（注意，该函数以一个 元组 为参数）
                 test((x,p1,p2,p3,p4,...))
        bound    求根区间及初始值(x_min,x_max,x_init)
        paralist (p1,p2,p3,p4,...)
    """
    (x_min,x_max,x_init) = bound;
    root_last = x_init;
    root_now = x_init;
    while True:
        root_now = root_last - func((root_last,) + paralist)/diff(func,root_last,paralist);
        if abs(root_last - root_now) <= 1e-5:
            break;
        if root_now >= x_min and root_now <= x_max:
            root_last = root_now;
        else:
            raise NoRootException(bound,paralist);
    return root_now;

def diff(func,x,paralist):
    """对指定的一元函数在指定的点求导数值

        func     函数对象，其定义形式应为，（注意，该函数以一个 元组 为参数）
                 test((x,p1,p2,p3,p4,...))
        x        自变量的值
        paralist (p1,p2,p3,p4,...)
    """
    dx = 1e-5;
    return (func((x + dx,) + paralist) - func((x - dx,) + paralist))/(2*dx);

class NoRootException(Exception):
   """求根过程自定义的异常"""
    def __init__(Self,bound,paralist):
        Exception.__init__(Self);
        Self.ebound = bound;
        Self.eparalist = paralist;

if __name__ == '__main__':
    def test((x,p1,p2,p3)):
        y = p1*x**2 + p2*x + p3;
        return y;
    try:
        print FZeros(test,(5,10,6),(1,-4,3));
    except NoRootException ,e:
        print "在所设置的区间" , e.ebound , "及参数表" , e.eparalist , "内没有根" ;

===================代码分割线======================

wqsong

不清楚，似乎没有。。。

Rainyboy

回复 3 # wqsong 的帖子

你是说循环试dx的大小么？

wqsong

回复 4 # Rainyboy 的帖子

嗯，调整步长，测试一下精度。。。
不知道python对浮点数的敏感度，也就是浮点数最小值。
设置一个步长在初始值和最小值之间缩放？

Rainyboy

回复 5 # wqsong 的帖子

说得倒是，确实应该改改，万一人家被导的函数本身就是1e-5量级的怎么办呢……

wqsong

Rainyboy 发表于 2010-12-5 23:31

                               
登录/注册后可看大图

回复 5 # wqsong 的帖子

说得倒是，确实应该改改，万一人家被导的函数本身就是1e-5量级的怎么办呢……

嗯嗯，哈哈，本身数量级1e-5。。。

smtmobly

optimize库中的fsolve函数可以用来对非线性方程组进行求解。它的基本调用形式如下：
fsolve(func, x0)
func(x)是计算方程组误差的函数，它的参数x是一个矢量，表示方程组的各个未知数的一组可能解，func返回将x代入方程组之后得到的误差；x0为未知数矢量的初始值。如果要对如下方程组进行求解的话：
f1(u1,u2,u3) = 0
f2(u1,u2,u3) = 0
f3(u1,u2,u3) = 0
那么func可以如下定义：

def func(x):
    u1,u2,u3 = x
    return [f1(u1,u2,u3), f2(u1,u2,u3), f3(u1,u2,u3)]
下面是一个实际的例子，求解如下方程组的解：

5*x1 + 3 = 0
4*x0*x0 - 2*sin(x1*x2) = 0
x1*x2 - 1.5 = 0
程序如下：

1. from scipy.optimize import fsolve
   
2. from math import sin,cos
   
3. 
   
4. def f(x):
   
5.     x0 = float(x[0])
   
6.     x1 = float(x[1])
   
7.     x2 = float(x[2])
   
8.     return [
   
9.         5*x1+3,
   
10.         4*x0*x0 - 2*sin(x1*x2),
    
11.         x1*x2 - 1.5
    
12.     ]
    
13. 
    
14. result = fsolve(f, [1,1,1])
    
15. 
    
16. print result
    
17. print f(result)

复制代码

输出为：

[-0.70622057 -0.6        -2.5       ]
[0.0, -9.1260332624187868e-14, 5.3290705182007514e-15]

smtmobly

不好意思，没有看到你们在这里讨论。

Rainyboy

回复 9 # smtmobly 的帖子

呵呵，谢谢！高手就是不一样啊，出手不凡，

smtmobly

http://hyry.dip.jp/pydoc
这儿有一本国人写的《用Python做数值计算》很不错，我发在论坛的很多代码是这里来的。大家可以看看，是一本Python方面很好的手册。使用与数值计算哈。

wqsong

回复 8 # smtmobly 的帖子

出手不凡。。。