song.yz@foxmail.com wechat: math-box

计算数学



非线性方程数值解法



Newton-Raphson迭代法

Newton-Rahpson迭代方法是解非线性方程比较著名而且也比较有效的方法之一。 如果初值比较接近根,收敛速度是很快的。Newton-Rahpson迭代法也是工程上广泛采用的方法。 Newton迭代法有比较直观的几何意义。

函数方程\(f(x) = 0\)的根是曲线\(y = f(x)\) 与\(x\)轴的交点的横坐标。 通过当前的点做曲线的切线(与一次导数有关),切线方程为 \[y = f({x_k}) + f'({x_k})(x - {x_k})\]

在上式中令\(y=0\)得到切线与\(x\)轴交点的横坐标 \[x = {x_k} - \frac{{f({x_k})}}{{f'({x_k})}}\] 于是我们可以构造迭代公式 \[{x_{k + 1}} = {x_k} - \frac{{f({x_k})}}{{f'({x_k})}}\] 如此反复迭代便可以比较快的得到方程的根,所以Newton法又叫切线法。

关于Newton迭代法如果函数 \(f(x)\)有二阶以上连续倒数,\(p\)是方程\(f(x) = 0\)的单根, 则当\(x\) 充分接近\(p\)时,Newton迭代法收敛,而且最少二阶收敛。

[例] 用Newton-Raphson方法计算方程 \[f(x) = {x^3} + 2{x^2} + 10x - 20 = 0\] 在1.5附近的根。下面给出go语言源代码。


// newton
package main

import (
	"fmt"
	"math"
)

func main() {
	/*------------------------------------------------------
	!  Author  : Song Yezhi      
	!  verison : 2020-5-18 10:29
	!  -----------------------------------------------------
	!  Input  Parameters :
	!
	!  Output Parameters :
	!
	------------------------------------------------------*/
	var x0 float64 = 0.0

	x, fx := newton(x0,myfx,dmyfx)
    //把函数作为参数传入
	fmt.Printf("the result: \n")
	fmt.Printf("x= %12.7f, f(x)= %12.7f \n", x, fx)
}

func myfx(x float64) float64 {
	/*------------------------------------------------------
	  !  Author  : Song Yezhi      
	  !  verison : 2020-5-18 10:28
	  !
	  ------------------------------------------------------*/
	fx := x*x*x + 2.0*x*x + 10.0*x - 20.0

	return fx
}

func dmyfx(x float64) float64 {
	/*------------------------------------------------------
	!  Author  : Song Yezhi      
	!  verison : 2020-5-18 10:29
	!
	------------------------------------------------------*/

	df := 3.0*x*x + 4.0*x + 10.0
	return df
}

func newton(x0 float64,funcX func(float64) float64,
          dfuncX func(float64) float64) (x1, fx float64) {
	/*------------------------------------------------------
	!  Author  : Song Yezhi      
	!  verison : 2020-5-18 10:30
	!      牛顿法计算方程根
	!  -----------------------------------------------------
	!  Input  Parameters :
	!       x0-----初值
    !       funcX------目标函数
    !       dfuncX-----导函数
	!  Output Parameters :
	!
	-------------------------------------------------------*/
	var imax int = 20
	var tol float64 = 1e-8

	for i := 0; i < imax; i++ {
		x1 = x0 - funcX(x0)/dfuncX(x0)
		fx = funcX(x0)
		dfx := dfuncX(x0)
		x1 = x0 - fx/dfx
		dx := math.Abs(x1 - x0)
		if dx < tol {
			break
		}
		x0 = x1
		fmt.Printf("i= %4d  x= %12.7f  f(x)=%12.7f \n", i, x1, fx)
	}
	return x1,fx
}


C语言版本如下

#include 
#include 

void func(double x, double &fx);
void dfunc(double x,double &df);
void   newtonRaphson( double x0, double &x, 
     double &root,   void(*func)(double x,double &fx),
     void (*dfun)(double x,double &df));

/*------------------------------------------------------
   Created  :  Song Yezhi   2023-3-18 11:41
 
   solve nonlinear eqation by newtonRaphson's method 
 	 
-------------------------------------------------------*/
int main()
{
    double x0 = 1.5;
    double x,fx ;
    newtonRaphson(x0,x,fx,func,dfunc);
}

/*------------------------------------------------------
   Created  :  Song Yezhi   2023-3-18 11:43
 
-------------------------------------------------------*/
void func(double x, double &fx)
{
    fx = x*x*x + 2.0*x*x +10.0*x -20 ;
}

/*------------------------------------------------------
   Created  :  Song Yezhi   2023-3-18 11:43
 
-------------------------------------------------------*/
void dfunc(double x,double &df)
{
    df = 3.0*x*x +4.0*x +10.0 ;
}

/*------------------------------------------------------
   Created  :  Song Yezhi   2023-3-18 11:43
       solve nonlinear eqation by newton Raphson's methods
--------------------------------------------------------
   Input Parameters   :
        
   Output Parameters  :

-------------------------------------------------------*/

void   newtonRaphson( double x0, double &x, 
     double &root,   void(*func)(double x,double &fx),
     void (*dfun)(double x,double &df))
{
    x =x0;
    int imax = 40;
    double tol = 1e-8;
    int i ;
    double fx, df ;
    for(i=0;i 小于 imax ;i++ )
    {
        func(x,fx);
        dfun(x,df);
        double dx = - fx/df ;
        printf("%3d %18.6lf %18.6lf \n",i, x,fx);
        if (fabs(dx) 小于 tol) break ;
        x = x + dx ;
    }
}