[TOC]

注意，我为二分法，迭代法与牛顿法，分别写了一个函数，它们都可以不带任何附加参数地处理高维情况。如果我下面展示的脚本可以得到正确的配置来运行，那么运行下面的这个脚本将会得到它们全部的结果：

from math import exp
from binary_search import BS
from p32_3_iteration import iter
from newton import newton_solve
from collections import namedtuple

def main():
    f = lambda x: exp(x) + 10 * x - 2
    mode = 1
    if mode == 0:
        inte = namedtuple('inteval', 'left right')
        inte.left = 0
        inte.right = 1
        BS(f, inte)
    elif mode == 1:
        g = lambda x:(2-exp(x))/10
        x0 = 0
        iter(x0, g, f)
    elif mode == 2:
        x0 = 0
        fp = lambda x:exp(x) + 10
        newton_solve(f, x0, fp=fp)
if __name__ == '__main__':
    import sys
    sys.exit(int(main() or 0))

所有代码基于 Python 3.5.2 ，外部依赖是 numpy 。

二分法

我定义了一个函数 BS，接收一个函数 f ，一个存放着左端点和右端点的 namedtuple ，和一个精确度 $\epsilon$ 作为参数。注意到最大迭代次数，也就是 global_step被我规定成了 1000 ，因为测试中不会遇到很大的迭代次数， 1000 是个比较适中的数值。我在这个函数里加入了一个异常检测，用来处理双根的情况。其他的异常情况，比如，迭代次数用尽，函数值接近于 0 ，也都有考虑。

代码如下：

from collections import namedtuple
import math
f = lambda x: x**3+4*x**2-10
def BS(f, inte, eps=1e-5):   

    global_step = 0
    tol = 1e-8
    while global_step < 1000:
        new_guess = (inte.right + inte.left)/2
        fn = f(new_guess)
        fl = f(inte.left)
        fr = f(inte.right)
        if math.fabs(fn) < tol:
            print('find a solosolo solution {}, function value:{}'.format(new_guess, fn))
            break        
        if inte.right-inte.left < eps:
            print('find solution with epsilon within 1e-5 {}, function value:{}, in timestep:{}'.format(new_guess, fn, global_step))
            break
        if fn*fl < -0:
            inte.right = new_guess
        elif fr*fn < -0:
            inte.left = new_guess
        else:
            raise ValueError('double root detected!')
        global_step += 1
    # else:
    #     print('no proper solution found')
    
def main():
    eps = 1e-5    
    inte = namedtuple('inteval', 'left right')
    inte.left = 1
    inte.right = 2
    BS(f, inte, eps)

if __name__ == '__main__':
    import sys
    sys.exit(int(main() or 0))

单元测试的输出为：

find solution with epsilon within 1e-5 1.3652305603027344, function value:9.030992742964372e-06, in timestep:17

迭代法

iter 这个函数接受四个参数，第一个参数 x0 为初始迭代点，第二个参数是迭代式，第三个参数是要求根的原函数，第四个默认参数是精确度。注意到为了兼容高维情况，所有的字面值，函数的参数及返回值都要求是 numpy.ndarray 的对象。最后跳出迭代循环的条件规定为相邻两个迭代值之差的 $L_2\ norm$ 小于给定的精确度。

import numpy as np

def iter(x0, f, ori_f, eps=1e-5):
    global_step = 0
    x0 = np.array(x0).reshape((-1, 1))
    prev_x = x0
    while global_step < 1000:
        new_x = f(x0)
        if np.linalg.norm(new_x-x0, ord=None) < eps:
            print('in step {}, solution value {} found, funtion value {}'.format(global_step, new_x, ori_f(new_x)))
            break
        x0 = new_x
        global_step += 1
    else:
        print('no proper root found')

def main():
    f1 = lambda x: np.sqrt(10/(4+x))  
    f2 = lambda x: 0.5*np.sqrt(10-x**3)
    ori_f = lambda x:x**3+4*x**2-10
    print('first iter')
    iter(1.5, f1, ori_f)
    print('second iter')
    iter(1.5, f2, ori_f)

if __name__ == '__main__':
    import sys
    sys.exit(int(main() or 0))

单元测试的输出为：

first iter
in step 5, solution value 1.3652305756734338 found, funtion value 9.28481537343373e-06
second iter
in step 15, solution value 1.3652332557424998 found, funtion value 5.354194795970102e-05

牛顿法

函数接受五个参数，f 为所求函数，x0为初始点， fp为导数的函数，注意，若导数没有显式给出，在函数内部将会用 $\frac{f(x+\epsilon)-f(x-\epsilon)}{2\epsilon}$ 计算导数的数值。eps为精确度。注意我使用了 np.linalg.tensorinv 来计算矩阵的逆。

import numpy as np

def newton_solve(f, x0, fp=None, eps=1e-5, N=1000):
    tol = 1e-8
    if not fp:
        fp = lambda x:(f(x+tol)-f(x-tol))/(2*tol)
    x = np.array(x0, ndmin=1)
    for i in range(N):
        fx = f(x)
        fpx = fp(x)
        if np.linalg.norm(fpx) < tol:
            raise ValueError("step {} derivative smaller than 1e-8".format())
        delta = np.dot(np.linalg.tensorinv(fpx, ), fx)
        print(delta)
        xnew = x - delta     
        if np.linalg.norm(delta) < eps:
            print("time step {} solution converge, x={}, f(x)={}".format(i, xnew, f(xnew)))
            return
        else:
            x = xnew
    print("no solution found in 1000 timesteps")
def main():
    f = lambda x: x**3+4*x**2-10
    fp = lambda x:3*x**2+8*x
    newton_solve(f, 1.5, fp)

if __name__ == '__main__':
    import sys
    sys.exit(int(main() or 0))

单元测试的输出为：

time step 3 solution converge, x=[ 1.36523001], f(x)=[ 0.]