逸珺 · 2020年09月14日

数学之美:牛顿-拉夫逊迭代法原理及其实现

首发:嵌入式客栈
作者:逸珺

导读

前面刚转了一篇文章提到了牛顿-拉夫逊(拉弗森)(Newton-Raphson method)方法,感觉这个数学方法很有必要相对深入写一篇文章来总结分享印证一下自己的理解。这是写本文的由来,如果发现文章中有错误之处,请留言交流讨论。

什么是牛顿-拉夫逊方法?

牛顿其人:Isaac Newton(1642年12月25日– 1727年3月20日) 是一位英国数学家,物理学家,天文学家,神学家和作家,被公认为有史以来最有影响力的科学家之一,并且是科学革命的关键人物。他的书《自然哲学的数学原理》于1687年首次出版,奠定了古典力学的基础。牛顿还为光学做出了开创性的贡献,并与戈特弗里德·威廉·莱布尼兹(Gottfried Wilhelm Leibniz)发展了无穷微积分的学科。

640.jpg

牛顿

拉弗森Joseph Raphson 生卒不详,其最著名的著作是1690年出版的《通用分析方程》。它包含一种方法,现在称其为牛顿-拉夫森方法,用于近似方程式的求根。艾萨克·牛顿(Isaac Newton)在1671年写的《通量法》中开发了一个非常相似的公式,但是这项工作要到1736年才出版,这是拉夫森分析之后近50年。但是,该方法的Raphson版本比Newton方法更简单,因此通常被认为是更好的方法。

640-1.jpg

所以,牛顿迭代法(简写)就是一种近似求解实数域与复数域求解方程的数学方法。那么这个方法是具体是什么原理呢?

牛顿迭代如何迭代?

直接看数学公式描述如何迭代不直观,先来看动图就很容易理解牛顿迭代法为什么叫迭代法以及怎样迭代的:

牛顿迭代法是原理是根据一个初始点在该点做切线,切线与X轴相交得出下一个迭代点的坐标,再在处做切线,依次类推,直到求得满足精度的近似解为止。

640.gif

由前面描述知道,牛顿迭代法是用来近似求解方程的,这里有两个点需要说明:

  • 为啥要近似求解?很多方程可能无法直接求取其解
  • 迭代法非常适合计算机编程实现,实际上计算机编程对于牛顿迭代法广为应用

来看看,数学上如何描述的?

其中为函数在处的一阶导数,也就是该点的切线。

来简单推一推上面公式的由来,直线函数方程为:

知道一个直线的一个坐标点以及斜率则该直线的方程就很容易可以得知:

那么该直线与轴的交点,就是也即等式的解:

啥时候停止迭代呢?

  1. 计算出
  2. 给出一个初始假定根值,利用上面迭代式子进行迭代
  3. 计算绝对相对迭代近似误差
  4. 将绝对相对近似误差与预定的相对误差容限进行比较。如果,则迭代步骤2,否则停止算法。另外,检查迭代次数是否已超过允许的最大迭代次数。如果是这样,则需要终止算法并退出。另一个终止条件是:

如何编码呢?

由于牛顿迭代法主要目的是解方程,当然也有可能用于某一个数学函数求极值,所以无法写出通用的代码,这里仅仅给出一个编代码的思路。相信掌握了思路,对于各种实际应用应该能很快的写出符合实际应用的代码。

假定一函数为

其波形图如下:

640.png

其一阶导数为:

那么对于该函数的根:

从图上大致可以知道有两个根,如果直接解方程,则很难求出其根,可以编个代码试试:

#include <stdio.h>#include <math.h>#include <stdlib.h>/*假定待求根函数如下*/#define    F(x)    (2*(x)*(x)-10*cos(x)+(x)-80)/*其一阶导数为*/#define   DF(x)    (4*(x)+10*sin(x)+1)float newton_rooting(float x0,float precision,float min_deltax,int max_iterations){     float xn,xn1,fn,fn1,dfn;     float deltax;     int step = 0;     xn  = x0;     xn1 = x0;     do{       xn  = xn1;       fn  = F(xn);       dfn = DF(xn);       /*判0*/       if( fabs(dfn) <1e-6 )       {            if( fabs(fn)>precision )                return NAN;            else                return fn;       }       xn1 = xn - fn/dfn;       fn1 = F(xn1);       deltax = fabs(xn1-xn);              step++;       if( step>max_iterations )       {           if( fabs(fn1)<precision )               return xn;           else               return NAN;       }     }     while( fabs(fn1)>precision || deltax>min_deltax );     return xn1;}void main(){     float root_guess = 23.0f;     float precision = 0.00001f;     float min_deltax = 0.001f;     float root;     int step = 7;     root = newton_rooting( root_guess,precision,min_deltax,step );     printf("根为: %f,函数值为:%f\n", root,F(root));     root_guess = -23;     root = newton_rooting( root_guess,precision,min_deltax,step );     printf("根为: %f,函数值为:%f\n", root,F(root));}

结果:

根为: 6.457232, 函数值为:0.000004根为: -6.894969,函数值为:-0.000008

函数值已经很接近于0了,如果还需要更为精确的值,则可以选择在此基础上进一步求解,比如利用二分法逼近。

需要注意些啥?

  • 求斜率可能为0,如为0时,则可能找到了函数的极值,比如:

640-1.png

  • 如果选择的初始猜测根的接近方程f(x)=0中函数f(x)的拐点 ,Newton-Raphson方法可能开始偏离根。然后,它可能会又收敛回到根。例如:

640-2.pngimage-20200913235435866

  • 如果选择的初值不合适,可能会跳掉一些根,比如:

640-3.png

所以实际应用时,需要知道自己待求解模型的大致情况,在合理的加以调整。

有哪些应用?

  • 比如知道某系统的传递函数,求解传函的参数,可以将上述方法推而广之,求解多维变量方程组,求导就变成求偏导了
  • 又比如设计一电路测量某物质的阻抗
  • ....

总结一下

牛顿迭代法在解决实际问题时,利用迭代求方程近似根的数学原理,在工程中有着很好的实用价值。比如求一个趋势的极值,传递函数参数辨识等都有广泛的实际应用。本文抛砖引玉,有可能文章也有很多错误疏漏的地方,如有不同看法或者发现错误,欢迎留言交流指正。

推荐阅读

推荐阅读
关注数
2882
内容数
266
分享一些在嵌入式应用开发方面的浅见,广交朋友
目录
极术微信服务号
关注极术微信号
实时接收点赞提醒和评论通知
安谋科技学堂公众号
关注安谋科技学堂
实时获取安谋科技及 Arm 教学资源
安谋科技招聘公众号
关注安谋科技招聘
实时获取安谋科技中国职位信息