《用Python学习数值分析--解方程》

2017-08-24

  对于二次方程,我们可以直接使用求根公式,这种方式称为“解析式”。对于更高次的方程,使用求根公式就困难多了,虽然难,但仍然是可行的;四次方程依然可行。但是,根据”伽罗瓦定理“告诉我们,五次及五次以上的方程没有解析解,你再也找不到公式来找到根了。相对于“解析式”,对于没有公式的方程,只能采用数值解法。本篇讨论单个方程数值求解的各种方法,用Python 的scipy算法库做简单介绍,深入的理论部分请参考《Numerical Analysis》教材。
  数值计算时,一个非常重大的问题是:什么时候选用什么算法?数值法需要选择起始点,计算区间[a, b]。怎么选呢?其实很简单,就是猜测。很不靠谱吧。
  常用方法有:

  1. 包围法(bracketing methods)
    • 二分法(Bisection Method)
    • 盈不足法(False position method)
  2. 迭代法
    • 牛顿法(Newton's Method)(类牛顿法有 Halley's method,Householder's method,)
    • 割线法(Secant method)
  3. 混合法
    • Brent's method 

  以方程f1 :1*x^5 + 2*x^4 + 3*x^3 + 4*x^2 + 5*x^1 + 6*x^0 = 0 为例。函数曲线图像如下所示。

 

一,包围法

  包围法,也翻译为“交叉法”,我觉得“交叉”这个词不够准确。这一类方法依赖于“介值定理”。当区间[a, b]足够小时,比如说单浮点数、双浮点数即将不能准确的表示a和b时,我们认为继续算下去没有什么意义了,可以停止了,那么,a或者b就可以被当作根了。当然,在实际计算时,截止计算时,区间[a, b]的精度也许会低,以减少计算量。包围法是非常朴素和原始的计算方法,都在两千年前被人们使用了。如下函数f2 = x * cos(x)

 

1,二分法
  这个应该相当直观了。如上面的函数f2,解的数量是无限多的。假设需要求解[-10, 10]之间的解。
  scipy.optimize.bisect(f, a, b, args=(), xtol=2e-12, rtol=8.8817841970012523e-16, maxiter=100, full_output=False, disp=True)
  该算法只需要三个参数,函数方程f, 猜测的寻根区间[a, b],算法复杂度为O(logN)

以scipy.optimize.bisect(f1, -10, 10, full_output=True) 调用时,
(-1.4917979881408883,       converged: True
           flag: 'converged'
 function_calls: 46
     iterations: 44
           root: -1.4917979881408883)

对于f2,计算出来以下解,很明显,这并不是我们想要的。
(0.0,       converged: True
           flag: 'converged'
 function_calls: 3
     iterations: 1
           root: 0.0)

2,盈不足法
  和二分法、割线法类似,该算法速度更快,不会像割线法一样会发散。不少学者认为这个方法来自于中国的盈不足术(例子请参考中文wiki)。这个方法能够处理的问题复杂多了,广义的盈不足术则指通过双假设法将其他数学问题转化为盈亏问题。能够解决非线性问题。
  在scipy中,并没有提供该算法,而是提供了它的一个变体:Ridders' method。

scipy.optimize.ridder(f1, -10, 10, xtol=2e-30, rtol=2e-15, maxiter=1000, full_output=True)
(-1.4917979881399008,       converged: True
           flag: 'converged'
 function_calls: 16
     iterations: 7
           root: -1.4917979881399008)
可以看到,三个重要的默认参数 xtol=_xtol, rtol=_rtol, maxiter=_iter,取值分别是:
_iter = 100                              # 计算迭代的次数
_xtol = 2e-12                          #绝对误差,必须是非负
_rtol = 4*finfo(float).eps    #相对误差   finfo(float).eps=8.881784197e-16
但是,单精度最小值是 1.175494351e-38,双精度最小值是2.22507e-308,所以,我们选择容忍的误差还是很大的。两个算法中相对误差rtol能够取的最小值是  finfo(float).eps,是double的epsilon。

  从上面的曲线图可以看到,根肯定在[-20, 20]之间,可能在[-5, 5]之间,但是,面对真正的问题时,你看不到函数曲线,如果你猜测的区间是[10, 20],可以看到你是绝对找不到根的,算法计算永远不会收敛。这也是数值算法的最大问题。相对于找到全局解,找到局部解容易的多。

二,迭代法
相比较于”bracketing methods“,最大的区别在于迭代法需要靠猜。感觉这事儿不靠谱吧。对,本来用解析法解决不了的问题,现在想要解决,肯定是需要付出某些代价的。
1,割线法
  这个算法可以看作是牛顿法的有限差分近似。但是,它被提出的时间要比牛顿法早3000年。牛顿法是二阶的,需要导数,3000年前的古人肯定是不知道微积分的概念的。割线法就相当于牛顿法的离散近似。scipy中并没有单独的函数来实现,而是内置于scipy.optimize.newton(unc, x0, fprime=None, args=(), tol=1.48e-08, maxiter=50, fprime2=None) 中了。如果提供了导数fprime,就使用牛顿法,否,则使用割线法。
2,牛顿法(也称 Newton–Raphson method)
  这个方法假设函数有连续的导数,如果初始点距离根很远,那么算法不会收敛。但,当收敛时,是二次收敛的,比二分法要快。
难点有:

  1. 计算函数的导数
  2. 在根附件收敛失败
  3. 在”驻点“处程序终止,无法完成计算
  4. 根大于1时收敛速度慢 

root = scipy.optimize.newton(f2, -10, f2_p, maxiter=3000)
print(root)

-10.9955742876

Halley's method要求二阶可导,Householder's method 要求最高d+1阶可导。 
所以,我们可以看到,或者在之后的解方程组算法中亦可见到,牛顿法是核心内容。是一定要熟练掌握的。

三,混合法
1,Brent's method
scipy 中 scipy.optimize.brent(func, args=(), brack=None, tol=1.48e-08, full_output=0, maxiter=500)

root = scipy.optimize.brent(f2,  full_output=True, maxiter=5000)
(-0.86033358893022038, -0.5610963381910451, 10, 11)

这个函数并不计算方程的解,而是给出最有点。现在还不知道这个算法有什么用途。

总结:
用数值方法解单变量高次方程(一般次数也不能太高,否则误差太大),是数值计算中最简单的问题。现实中的问题建模为单变量方程的情况较少,极少有系统是受单因素影响。多变量系统才是我们关注的重点,而且,也只是多变量线性系统,多变量非线性系统我们一般很难解,也没有什么意义。我们做不到的事情真的还有很多。

ref:

  1. http://zh.numberempire.com/graphingcalculator.php
  2. https://en.wikipedia.org/wiki/Root-finding_algorithm
  3. https://en.wikipedia.org/wiki/Newton%27s_method
  4. https://en.wikipedia.org/wiki/Householder%27s_method
  5. https://en.wikipedia.org/wiki/Secant_method
  6. https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.bisect.html
  7. https://en.wikipedia.org/wiki/Brent%27s_method
如果有任何意见,欢迎留言讨论。


[ 主页 ]
COMMENTS
POST A COMMENT

(optional)



(optional)