Astrodynamics

非线性方程数值解法

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
}

重根迭代改进

Newton迭代法在处理重根运算时候往往效果不是太好，有时候收敛速度会非常的慢。为了处理这一问题，往往我们需要对迭代法进行一些改造。

如果采取以下的迭代方法，计算效果会很不错，但是需要计算二次导数。 \[{x_k} = {x_{k - 1}} - \frac{{F({x_{k - 1}})}}{{F'({x_{k - 1}})}}\] 其中 \[F(x) = \frac{{f(x)}}{{f'(x)}}\] 展开上式，便得到迭代公式 \[{x_k} = {x_{k - 1}} - \frac{{f({x_{k - 1}})f'({x_{k - 1}})}}{{{{[f'({x_{k - 1}})]}^2} - f({x_{k - 1}})f''({x_{k - 1}})}}\] 这就是重根时的改进迭代格式。

用重根迭代改进方法计算方程 \[f(x) = {x^4} - 4{x^2} + 4 = 0\] 在\(x=1.5\)附近的根。下面给出go语言代码。


// multiroot
package main
import (
	"fmt"
	"math"
)
func main() {
	/*------------------------------------------------------
	!  Author  : Song Yezhi      
	!  verison : 2020-5-18 10:29
	!  -----------------------------------------------------
	!  Input  Parameters :
	!
	!  Output Parameters :
	!
	------------------------------------------------------*/
	var x0 float64 = 1.5    
    fmt.Printf("newton result: \n")
	_, _  = newton(x0)
	fmt.Printf("-------------- \n")	
	fmt.Printf("multiroot iter : \n")
	_, _  = multiRoot(x0)
}

func funcX(x float64) float64 {
	/*------------------------------------------------------
	  !  Author  : Song Yezhi      
	  !  verison : 2021.10.04
	  !
	  ------------------------------------------------------*/
	fx := math.Pow(x,4)-4*x*x+4
	return fx
}
func dfuncX(x float64) float64 {
	/*------------------------------------------------------
	!  Author  : Song Yezhi      
	!  verison : 2020-5-18 10:29
	!
	------------------------------------------------------*/
	df := 4.0*x*x*x - 8.0*x 
	return df
}
func d2funcX(x float64) float64{
    /*------------------------------------------------------
     Author  : Song Yezhi                       
     verison : 2021-10-4 18:04
     go build -gcflags "-N -l"  
     
    ------------------------------------------------------*/
    d2f:=12.0*x*x-8.0
    return d2f    
}
func newton(x0 float64) (x1, fx float64) {
	/*------------------------------------------------------
	!  Author  : Song Yezhi      
	!  verison : 2020-5-18 10:30
	!      牛顿法计算方程根
	!  -----------------------------------------------------
	!  Input  Parameters :
	!       x0-----初值
	!  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
}
func  multiRoot(x0 float64)(x1,fx float64){
/*------------------------------------------------------
  Author  : Song Yezhi                       
  verison : 2021-10-4 18:08
  go build -gcflags "-N -l"  
     
  -----------------------------------------------------  
  Input  Parameters :
     x0  --- 初值
  Output Parameters : 
     x1----- 根
     fx -----函数值
------------------------------------------------------*/
    imax:= 200
    tol:= 1e-8
    
    var tmp1,tmp2  float64
    var dx float64
    for i:=0;i < imax; i++ {
        tmp1 = funcX(x0)*dfuncX(x0)
        tmp2 = dfuncX(x0)*dfuncX(x0)-funcX(x0)*d2funcX(x0)
        x1 = x0 -tmp1/tmp2
        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
}

Muhammad两步三阶迭代格式

对于非线性方程\({F\left( x \right) = 0}\)，对方程函数在\({{x_n}}\)处泰勒展开得： \[\begin{gathered} F\left( x \right) = F\left( {{x_n}} \right) + \cdots + \frac{1}{{\left( {k - 1} \right)!}}{F^{\left( {k - 1} \right)}}\left( {{x_n}} \right){\left( {x - {x_n}} \right)^{\left( {k - 1} \right)}} + \hfill \\ \quad \quad \quad \int_0^1 {\frac{{{{\left( {1 - t} \right)}^{\left( {k - 1} \right)}}}}{{\left( {k - 1} \right)!}}{F^{\left( k \right)}}\left[ {{x_n} + t\left( {x - {x_n}} \right)} \right]} {\left( {x - {x_n}} \right)^k}dt \hfill \\ \end{gathered} \] 当\(k=1\)时，有 \[F\left( x \right) = F\left( {{x_n}} \right) + \int_0^1 {F'\left[ {{x_n} + t\left( {x - {x_n}} \right)} \right]} \left( {x - {x_n}} \right)dt\] 对上式采用不同的数值积分方法可以得到一系列迭代公式。 Muhammad Aslam Noor给出了如下两种的积分公式 \[\int_0^1 {F'\left[ {{x_n} + t\left( {x - {x_n}} \right)} \right]} \left( {x - {x_n}} \right)dt \approx \left[ {\frac{1}{4}F'\left( {{x_n}} \right) + \frac{3}{4}F'\left( {\frac{{{x_n} + 2x}}{3}} \right)} \right]\left( {x - {x_n}} \right)\] \[\int_0^1 {F'\left[ {{x_n} + t\left( {x - {x_n}} \right)} \right]} \left( {x - {x_n}} \right)dt \approx \left[ {\frac{3}{4}F'\left( {\frac{{2{x_n} + x}}{3}} \right) + \frac{1}{4}F'\left( {{x_n}} \right)} \right]\left( {x - {x_n}} \right)\] 由牛顿法作为第一步，由此可以给出如下两种解算非线性方程的两步迭代格式： \[\left\{ \begin{gathered} {y_n} = {x_n} - F'{\left( {{x_n}} \right)^{ - 1}}F\left( {{x_n}} \right) \hfill \\ {x_{n + 1}} = {x_n} - \frac{{F\left( {{x_n}} \right)}}{{\frac{1}{4}F'\left( {{x_n}} \right) + \frac{3}{4}F'\left( {\frac{{{x_n} + 2{y_n}}}{3}} \right)}} \hfill \\ \end{gathered} \right.\quad \quad n = 1,2, \cdots \] \[\left\{ \begin{gathered} {y_n} = {x_n} - F'{\left( {{x_n}} \right)^{ - 1}}F\left( {{x_n}} \right) \hfill \\ {x_{n + 1}} = {x_n} - \frac{{F\left( {{x_n}} \right)}}{{\frac{3}{4}F'\left( {\frac{{2{x_n} + {y_n}}}{3}} \right) + \frac{1}{4}F'\left( {{y_n}} \right)}}\quad \hfill \\ \end{gathered} \right.n = 1,2, \cdots \] 以上格式为三阶收敛，且同样适用于非线性方程组的计算。计算非线性方程 \[f\left( x \right) = 3{x^5} - 2{x^3} + 6x - 8\] 的根。容易给出 \[\frac{{df\left( x \right)}}{{dx}} = 15{x^4} - 6{x^2} + 6\] 下面给出go语言代码。


// muhammad
package main

import (
	"fmt"
	"math"
)

func main() {
	/*------------------------------------------------------
	!  Author  : Song Yezhi      
	!  verison : 2021-10-04 19:02:15
	!  -----------------------------------------------------
	!  Input  Parameters :
	!
	!  Output Parameters :
	!
	------------------------------------------------------*/
	var x0 float64 = 1.5
	fmt.Printf("muhammad 1 method: \n")
	_, _= muhammad1(x0)
	
	fmt.Printf("------------------ \n")
	fmt.Printf("muhammad 2 method: \n")
	_, _= muhammad2(x0)
}
func muhammad1(x0 float64)(x1,fx float64){
/*------------------------------------------------------
  Author  : Song Yezhi                       
  verison : 2021-10-4 18:51
  go build -gcflags "-N -l"  
     
  -----------------------------------------------------  
  Input  Parameters :
 
  Output Parameters :
 
------------------------------------------------------*/
    imax := 200
    tol := 1e-8
    x1 = x0
    
    var y1,tmp,dx,x2 float64
    for i:=0;i < imax ; i++ {
        fx = funcX(x1)
        y1 = x1 - fx/dfuncX(x1)
        tmp = 0.25*dfuncX(x1) +0.75 * dfuncX((x1+2*y1)/3.0)
        x2 = x1 - fx/tmp
        dx = math.Abs(x2-x1)
        if dx < tol {
            break
        }
        x1 = x2 
        fmt.Printf("i= %4d  x= %12.7f  f(x)=%12.7f \n", i, x2, fx)       
    }
    return x2,fx
}

func muhammad2(x0 float64)(x1,fx float64){
/*------------------------------------------------------
  Author  : Song Yezhi                       
  verison : 2021.10.04 
  go build -gcflags "-N -l"  
     
  -----------------------------------------------------  
  Input  Parameters :
 
  Output Parameters :
 
------------------------------------------------------*/
    imax := 200
    tol := 1e-8
    x1 = x0
    
    var y1,tmp,dx,x2 float64
    for i:=0;i < imax ; i++ {
        fx = funcX(x1)
        y1 = x1 - fx/dfuncX(x1)
        tmp = 0.75*dfuncX((2.0*x1+y1)/3.0) +0.25 * dfuncX(y1)
        x2 = x1 - fx/tmp
        dx = math.Abs(x2-x1)
        if dx < tol {
            break
        }
        x1 = x2 
        fmt.Printf("i= %4d  x= %12.7f  f(x)=%12.7f \n", i, x2, fx)       
    }
    return x2,fx
}
func funcX(x float64) float64 {
	/*------------------------------------------------------
	  !  Author  : Song Yezhi      
	  !  verison : 2021.10.04
	  !
	  ------------------------------------------------------*/
	fx := 3.0*math.Pow(x,5)-2.0*x*x*x + 6*x -8.0
	return fx
}
func dfuncX(x float64) float64 {
	/*------------------------------------------------------
	!  Author  : Song Yezhi      
	!  verison : 2021.10.04
	!
	------------------------------------------------------*/
	df := 15.0*math.Pow(x,4)-6.0*x*x + 6.0
	return df
}

计算数学

numerical analysis

非线性方程数值解法

Newton-Raphson迭代法

重根迭代改进

Muhammad两步三阶迭代格式