Report 1
[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.]