你的浏览器版本过低,可能导致网站不能正常访问!
为了你能正常使用网站功能,请使用这些浏览器。

牛顿-拉夫逊(拉弗森)(Newton-Raphson method)方法

[复制链接]
gaosmile 发布时间:2020-9-15 10:14

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

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

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

0 p  E4 I4 s; \
微信图片_20200915100704.jpg 牛顿
4 d# V+ K1 S/ g* R1 Y, C2 Y

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

微信图片_20200915100708.jpg
" Y. M8 q3 c5 \

所以,牛顿迭代法(简写)就是一种近似求解实数域与复数域求解方程的数学方法。那么这个方法是具体是什么原理呢?: F6 f  V, {, L" j2 a% s3 s/ Y

牛顿迭代如何迭代?

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

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

# |5 c6 x2 `4 {6 p  _

微信图片_20200915100711.gif

由前面描述知道,牛顿迭代法是用来近似求解方程的,这里有两个点需要说明:
3 }. _$ K0 w# V! W

  • 为啥要近似求解?很多方程可能无法直接求取其解
  • 迭代法非常适合计算机编程实现,实际上计算机编程对于牛顿迭代法广为应用. ^  Y! j( q( k$ E, H% F6 C# _

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

1600135792(1).png

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

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

1600135824(1).png

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

1600135844.png

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

1600135864(1).png

啥时候停止迭代呢?
  • 计算出
  • 给出一个初始假定根值,利用上面迭代式子进行迭代 1600135885(1).png
  • 计算绝对相对迭代近似误差 1600135936(1).png

  • 将绝对相对近似误差与预定的相对误差容限进行比较。如果\epsilon_s" data-formula-type="inline-equation" style="max-width: 100%; box-sizing: border-box !important;">,则迭代步骤2,否则停止算法。另外,检查迭代次数是否已超过允许的最大迭代次数。如果是这样,则需要终止算法并退出。另一个终止条件是: 1600135953(1).png
    / a5 N! {: W8 E! Y5 i" }* W

4 ?0 a1 _  @3 J3 F' {9 E

6 j. f- n. z+ H1 ~! e8 Z如何编码呢?

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

假定一函数为

1600135978(1).png

其波形图如下:

微信图片_20200915100715.png

其一阶导数为:

1600135996(1).png

那么对于该函数的根:

1600136018(1).png

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

#include <stdio.h>" h5 ~) }& G+ \/ ?' m- k
#include <math.h>
! l% d: F7 j  r+ {- a- Y3 P3 e/ k( ]#include <stdlib.h>
# @- c9 ]4 |; S, x8 H+ Q
# K. h  C; T; I* h) P2 G/*假定待求根函数如下*/
6 ?* z) y) ?; k9 F: g- t* s#define    F(x)    (2*(x)*(x)-10*cos(x)+(x)-80)
: ]/ s* K5 N, b! t0 F% t7 T! k" t
/*其一阶导数为*/
4 P) s0 }: G$ o- b#define   DF(x)    (4*(x)+10*sin(x)+1). ~# N% u6 s7 X& d3 b4 A
4 {$ Q: Q/ L4 f+ ~; `, Q0 v
float newton_rooting(float x0,float precision,float min_deltax,int max_iterations)5 k. e7 P; k9 j4 {3 Y$ J! }
{
; r& C% p; O; z- w- M; Y& h' j     float xn,xn1,fn,fn1,dfn;
) q5 s! m- X9 s( _% ~     float deltax;* M. R' h4 f* I* I0 A
     int step = 0;
9 g' G! H) R6 w6 i+ _     xn  = x0;
" }; B) H  R9 W+ _     xn1 = x0;
7 \$ w: y) J) V, `9 C6 O" t$ Y     do{
' ]! H! p, I, x8 D! Z. \       xn  = xn1;
, h" a% R7 D8 A9 z       fn  = F(xn);+ U6 ]' H% u8 Q7 |
       dfn = DF(xn);
% D+ a- Z- j0 H& i% W9 f( C       /*判0*/
/ p- V% {$ H8 q: s/ H, L2 o       if( fabs(dfn) <1e-6 )
7 O" L9 w- N5 T1 d       {
! S% g/ \( k2 o) E( q" r            if( fabs(fn)>precision )
" v" ^2 t0 C* d5 g+ t                return NAN;
7 w! p( H9 K9 }# F) `            else
/ K. M* y. @2 y                return fn;0 @$ H! H3 N7 k) F# v
       }
: i- ~2 P$ [' f5 p7 t
* p3 M# B+ Q' V, m$ T       xn1 = xn - fn/dfn;
* T# q* I$ s0 E3 R, _       fn1 = F(xn1);( d8 z1 f! Q5 s+ D2 I! F& W- B
       deltax = fabs(xn1-xn);: H. K# L: @) E
      
6 S" Y0 c5 _0 ^' h( c7 T3 S       step++;
1 M0 `6 z" \( k1 U' c9 u! }- T7 Z       if( step>max_iterations )+ E2 H8 a% D( G) p5 o; [
       {
2 }& A/ h0 A$ T) x9 q+ K; g           if( fabs(fn1)<precision ): Y( \: J0 r" e, k# n
               return xn;* q+ k7 \, Y7 ]: h4 j
           else
) ~$ R: D& a  }$ {7 W# X8 U2 S/ _               return NAN;
$ E. }9 t7 f8 Y       }: {  X. h* i5 I  M: o5 ]  ]
     }
# A( C7 h0 W7 [  D* F/ a     while( fabs(fn1)>precision || deltax>min_deltax );
, l4 R+ y' J9 u
7 u' `1 n9 W0 v* V     return xn1;- k7 S2 s. s7 |. C8 s, [. ~0 i
}% S3 b1 D% W$ n" R! `
, G. F- S" Z8 _/ ]7 |; C
void main()
9 B, R( l' r) r+ j; J/ Q{) ~: e' @0 O) C# G3 j1 T/ }
     float root_guess = 23.0f;( n( x% x: E, U) u* ]0 @
     float precision = 0.00001f;
' ^& ]" Q4 D6 q8 Y& z; Z     float min_deltax = 0.001f;
0 b$ I* o8 T( F# _- ^& Q     float root;( N" f8 V1 R( u% ~) ]( h% O
     int step = 7;
) _: s1 w) Y; P" X! ?. E) ?9 u7 Z( f  ], B. k
     root = newton_rooting( root_guess,precision,min_deltax,step );) \4 l0 t9 L, e! S4 w% y/ S1 P
     printf("根为: %f,函数值为:%f\n", root,F(root));
* P9 v  e+ Z& u8 ^/ }( `) I. ]
2 z$ w+ W1 o: s7 Y# ]     root_guess = -23;# I9 l) K+ L9 n/ w; r/ {: T
     root = newton_rooting( root_guess,precision,min_deltax,step );& D+ E* r) ^/ l" ~/ Q
     printf("根为: %f,函数值为:%f\n", root,F(root));  u; A- _: |, p, H% P7 Z+ ]; u
}
" K1 ]: c% Y# c) x& ?, _4 I

结果:

根为: 6.457232, 函数值为:0.000004
. E$ E1 Y$ j: J) F根为: -6.894969,函数值为:-0.000008
/ |) e, Y' H9 i2 ~7 P

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

需要注意些啥?
  • 求斜率可能为0,如为0时,则可能找到了函数的极值,比如:


    ) \+ g" h$ J- }, [4 Z% _  E
微信图片_20200915100718.png
  • 如果选择的初始猜测根的接近方程f(x)=0中函数f(x)的拐点 ,Newton-Raphson方法可能开始偏离根。然后,它可能会又收敛回到根。例如:) D, C) |# B- x+ z) D9 [" K- `
微信图片_20200915100721.png 5435866
  • 如果选择的初值不合适,可能会跳掉一些根,比如:

    微信图片_20200915100725.png
    ( ~7 A: w2 }% E2 M" u+ s

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

有哪些应用?
  • 比如知道某系统的传递函数,求解传函的参数,可以将上述方法推而广之,求解多维变量方程组,求导就变成求偏导了
  • 又比如设计一电路测量某物质的阻抗
  • ....8 _( x7 C, W% a3 {! r
总结一下

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

  ~4 ]( |: i  B. d
" K7 ?6 v) h+ H  `. e, J! `4 j
收藏 1 评论1 发布时间:2020-9-15 10:14

举报

1个回答
Kevin_G 回答时间:2020-9-15 11:40:00
大大的点赞,好厉害!

所属标签

关于
我们是谁
投资者关系
意法半导体可持续发展举措
创新与技术
意法半导体官网
联系我们
联系ST分支机构
寻找销售人员和分销渠道
社区
媒体中心
活动与培训
隐私策略
隐私策略
Cookies管理
行使您的权利
官方最新发布
STM32N6 AI生态系统
STM32MCU,MPU高性能GUI
ST ACEPACK电源模块
意法半导体生物传感器
STM32Cube扩展软件包
关注我们
st-img 微信公众号
st-img 手机版