牛顿迭代法

Posted by

如果x \mapsto f(x)是一个可微函数,那么方程g(x) = 0的一个解就是函数x \mapsto f(x)的一个不动点,其中:

f(x) = x - \frac{g(x)}{Dg(x)}                     (1)

Dg(x)是g对x的导数,对于许多函数(有一些函数或者特殊情况会使牛顿法无法求解),以及充分好的初始猜测x,牛顿法都能很快收敛到g(x)=0的一个解.

一般而言,如果g是一个函数而dx是一个很小的数,那么g的导数在任一数值x的值由下面函数(作为很小的数dx的极限)给出:

Dg(x) = \frac{g(x+dx) -g(x)}{dx}

对于(1)进行推导.

x_n的切线方程为g(x_n) = f'(x_n) ( x - x_n ) + f(x_n), x_{n+1}g(x_n) = 0的解,

f'(x_n) ( x - x_n ) + f(x_n) = 0

x = x_n - \frac{f(x_n)}{f'(x_n)}

x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}

用图像表示牛顿迭代法的过程


对曲线f(x) = \frac{(x - 5) ^ 3}{10} - 1以x = 10为初始点的应用牛顿迭代法求零点,其中红色为函数切线,蓝色为垂线.B,D,F为三个切点,A,C,E为切点在x轴上的投影,可以发现每一次投影的点都在向曲线与x轴的交点接近.

用scheme写的牛顿迭代法求开方代码

(define dx 0.00001 )

(define tolerance 0.00001 )

(define (fixed-point f first-guess )
    (define (close-enough? v1 v2  )
        ( < ( abs ( - v1 v2 ) ) tolerance ))
    (define (try guess  )
        (let ( (next ( f guess ) ) )
            (if ( close-enough? guess next )
                next
                ( try next )
            )
        )
    )
    ( try first-guess )
)

(define (cube x ) ( * x x x ))

(define (deriv g )
    (lambda (x) 
        ( / ( -  (g ( + x dx ) ) (g x ) ) dx ) 
    )
)

(define (newton-transform g )
    (lambda (x)
        ( - x ( / ( g x ) (( deriv g ) x) ) )
    )
)

(define (newton-method g guess )
    ( fixed-point ( newton-transform g ) guess )
)


(define (fixed-point-of-transform g transform guess )
    ( fixed-point ( transform g ) guess )
)

(define (sqrt3 x )
    (fixed-point-of-transform 
        (lambda (y) ( - ( square y ) x )  )
        newton-transform
        1.0
    )
)

( sqrt3 9 )

用c写的牛顿迭代法求开方的代码

#include <bits/stdc++.h>

using namespace std;
#define dx 0.000001
#define minx 0.0000001

double f( double x ) {return x * x - 2;}

double Df( double x ) {return (f( x + dx ) - f(x)) / dx;}

double next_point ( double x ) { return x - f(x) / Df(x); }

double find_ans(double x)
{
    double point1 = x;
    double point2;
    for ( ;; )
    {
        point2 = next_point(point1);
        if ( abs(point1 - point2) < minx ) break;
        point1 = point2;
    }
    return point2;
}

int main()
{
    cout << find_ans(2)<<endl;
    return 0;
}

一些牛顿法不能求解的例子

在迭代过程中遇到了驻点

y = |x| ^ {1/2}的两个切点之间循环震荡的情况

y = x^{1/3}的切线投影与函数零点越来越远,不收敛于零点的情况

诸如这些的还有在多个零点情况下不能求出多解,或者求出非期望的零点等等,不再一一举例.

Leave a Reply

电子邮件地址不会被公开。 必填项已用*标注