[数值分析]牛顿迭代法Newton!!!

牛顿迭代法Newton!!!

牛顿迭代法的基本思想

牛顿迭代法是一种用来求解方程的方法，它的基本思想是：如果一个函数在某一点的切线是直线，那么迭代下一次产生的值就是切线与x轴的交点的x坐标，按照这个思想，我们可以不断迭代，直到收敛到方程的根。

牛顿迭代法的数学表达

设函数$f(x)$在$x_0$处可导，且$f(x)$在$x_0$处的一阶导数$f'(x_0)$存在，且$f'(x_0)\neq 0$，则称$x_0$是方程$f(x)=0$的一个牛顿迭代点，记作$x_0\in\mathbb{R}$，且称函数$f(x)$在$x_0$处的切线方程为

$$f(x)=f(x_0)+f'(x_0)(x-x_0)$$

若$x_1$是方程$f(x)=0$的一个牛顿迭代点，且$x_1$满足

$$x_1=x_0-\frac{f(x_0)}{f'(x_0)}$$

则称$x_1$是方程$f(x)=0$的一个牛顿迭代点，记作$x_1\in\mathbb{R}$，且称函数$f(x)$在$x_1$处的切线方程为

$$f(x)=f(x_1)+f'(x_1)(x-x_1)$$

若$x_2$是方程$f(x)=0$的一个牛顿迭代点，则$x_2$也满足

$$x_2=x_1-\frac{f(x_1)}{f'(x_1)}$$

所以，牛顿迭代法的迭代公式为

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

简单例子

求解方程$x^3 - x^2 - x - 1 = 0$的正根，$\varepsilon=10^{-4}$

用牛顿迭代法

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

def f(x):
    return x**3 - x**2 - x - 1

def f1(x):
    return 3*x**2 - 2*x - 1

def newton(x0, eps):
    x = x0
    while abs(f(x)) > eps:
        x = x - f(x)/f1(x)
    return x

newton(1, 1e-4)

newton(2, 1e-4)

from pickle import TRUE
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline


def f(x):
    return x**3 - x**2 - x - 1


def f1(x):
    return 3*x**2 - 2*x - 1


def newton(x0, eps,isprint=True):
    x = x0
    while abs(f(x)/f1(x)) > eps:
        x = x - f(x)/f1(x)
        if isprint:
            print(x)
    x = x - f(x)/f1(x)
    return x


# newton(1, 1e-4) 

x = np.linspace(-2, 3, 100)
plt.plot(x, f(x))
plt.plot(x, f1(x))
plt.grid()
# 显示坐标轴
plt.axhline(y=0, color='k')
plt.axvline(x=0, color='k')

<matplotlib.lines.Line2D at 0x2b08f696cd0>

​

​

f(1)

-2

f(2)

1

newton(2, 1e-4)

1.8571428571428572
1.839544513457557
1.8392868100680193





1.8392867552141636

tmpx = 1.839544513457557
f(tmpx)/f1(tmpx)

0.00025770338953776654

弦截法

def secant(x0, x1, eps):
    x = x1
    while abs(f(x)) > eps:
        x = x - f(x)*(x-x0)/(f(x)-f(x0))
        x0 = x1
        x1 = x
    return x

secant(1, 2, 1e-4)

def secant(x0, x1, eps, isprint=True):
    x = x1
    while abs(f(x)*(x-x0)/(f(x)-f(x0))) > eps:
        x = x - f(x)*(x-x0)/(f(x)-f(x0))
        x0 = x1
        x1 = x
        if isprint:
            print(x)
    x = x - f(x)*(x-x0)/(f(x)-f(x0))
    return x

secant(1, 2, 1e-4)

1.6666666666666667
1.816326530612245
1.8429939112321743
1.8392156332109642





1.839286537939871

小结

相比于弦截法，牛顿迭代法的收敛速度更快，但是牛顿迭代法需要求解一阶导数，而弦截法不需要。

但牛顿迭代法的一阶导数为零时，牛顿迭代法会出现问题，所以牛顿迭代法的收敛速度更快的前提是一阶导数不为零。

我们在后面的课程中会学习到牛顿迭代法的改进方法，这些改进方法可以解决牛顿迭代法一阶导数为零时的问题。

牛顿迭代法解方程组

简单例子

$$\begin{cases} f(x,y) = x + 2y - 3 = 0\\ g(x,y) = 2x^2 + y^2 - 5 = 0\\ (x_0, y_0) = (-1, 2) \end{cases}$$

对上述方程组，我们求其迭代过程中的一阶展开式

$$\begin{bmatrix} f(x_{n+1}, y_{n+1})\\ g(x_{n+1}, y_{n+1}) \end{bmatrix} = \begin{bmatrix} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y}\\ \frac{\partial g}{\partial x} & \frac{\partial g}{\partial y} \end{bmatrix} \begin{bmatrix} x_{n+1} - x_n\\ y_{n+1} - y_n \end{bmatrix} + \begin{bmatrix} f(x_n, y_n)\\ g(x_n, y_n) \end{bmatrix}$$

分别令$f(x_{n+1}, y_{n+1})=0$和$g(x_{n+1}, y_{n+1})=0$，且记$\Delta x=x_{n+1}-x_n$和$\Delta y=y_{n+1}-y_n$，则有

$$\begin{cases} f(x_n, y_n) + \frac{\partial f}{\partial x}\Delta x + \frac{\partial f}{\partial y}\Delta y = 0\\ g(x_n, y_n) + \frac{\partial g}{\partial x}\Delta x + \frac{\partial g}{\partial y}\Delta y = 0 \end{cases}$$

则我们的目标成为了解下述方程

$$\begin{align} \frac{\partial f}{\partial x} \Delta x + \frac{\partial f}{\partial y} \Delta y &= -f(x_n,y_n)\\ \frac{\partial g}{\partial x} \Delta x + \frac{\partial g}{\partial y} \Delta y &= -g(x_n,y_n) \end{align}$$

这是一个二元一次方程组，由高中知识我们可以很容易地解出$\Delta x$和$\Delta y$，从而得到$x_{n+1}$和$y_{n+1}$。

$$\begin{align} \Delta x &= - \frac{\frac{\partial g}{\partial y}f(x_n,y_n) - \frac{\partial f}{\partial y}g(x_n,y_n)}{\frac{\partial f}{\partial x}\frac{\partial g}{\partial y} - \frac{\partial f}{\partial y}\frac{\partial g}{\partial x}}\\ \Delta y &= - \frac{\frac{\partial f}{\partial x}g(x_n,y_n) - \frac{\partial g}{\partial x}f(x_n,y_n)}{\frac{\partial f}{\partial x}\frac{\partial g}{\partial y} - \frac{\partial f}{\partial y}\frac{\partial g}{\partial x}} \end{align}$$

又有$x_{n+1}=x_n+\Delta x$和$y_{n+1}=y_n+\Delta y$，所以我们的迭代公式为

$$\begin{cases} x_{n+1} = x_n - \frac{f(x_n, y_n)g_y(x_n, y_n) - g(x_n, y_n)f_y(x_n, y_n)}{f_x(x_n, y_n)g_y(x_n, y_n) - f_y(x_n, y_n)g_x(x_n, y_n)}\\ y_{n+1} = y_n - \frac{f_x(x_n, y_n)g(x_n, y_n) - g_x(x_n, y_n)f(x_n, y_n)}{f_x(x_n, y_n)g_y(x_n, y_n) - f_y(x_n, y_n)g_x(x_n, y_n)} \end{cases}$$

则对上述具体例子，我们便可以首先求解$f(x,y)$和$g(x,y)$的偏导数，容易得到

$$\begin{align} \frac{\partial f}{\partial x} &= 1\\ \frac{\partial f}{\partial y} &= 2\\ \frac{\partial g}{\partial x} &= 4x\\ \frac{\partial g}{\partial y} &= 2y \end{align}$$

Just show code and the result

def f(x, y):
    return x + 2*y - 3

def g(x, y):
    return 2*x**2 + y**2 - 5

def f1(x, y):
    return 1

def f2(x, y):
    return 2

def g1(x, y):
    return 4*x

def g2(x, y):
    return 2*y

def newton(x0, y0, eps，is_print=True):
    x = x0
    y = y0
    while abs(f(x, y)) > eps or abs(g(x, y)) > eps:
        x = x - (f(x, y)*g2(x, y) - g(x, y)*f2(x, y))/(f1(x, y)*g2(x, y) - f2(x, y)*g1(x, y))
        y = y - (g(x, y)*f1(x, y) - f(x, y)*g1(x, y))/(f1(x, y)*g2(x, y) - f2(x, y)*g1(x, y))
        if is_print:
            print(x, y)
    return x, y

newton(-1, 2, 1e-4)

def f(x, y):
    return x + 2*y - 3


def g(x, y):
    return 2*x**2 + y**2 - 5


def f1(x, y):
    return 1


def f2(x, y):
    return 2


def g1(x, y):
    return 4*x


def g2(x, y):
    return 2*y


def newton(x0, y0, eps, is_print=True):
    x = x0
    y = y0
    while True:
        if is_print:
            print(x, y)
            print(f1(x, y), f2(x, y), f(x, y))
            print(g1(x, y), g2(x, y), g(x, y))
        dx = - (f(x, y)*g2(x, y) - g(x, y)*f2(x, y))/(f1(x, y)*g2(x, y) - f2(x, y)*g1(x, y))
        dy = - (f1(x, y)*g(x, y) - f(x, y)*g1(x, y))/(f1(x, y)*g2(x, y) - f2(x, y)*g1(x, y))
        x = x + dx
        y = y + dy
        if is_print:
            print(dx , dy)
        if abs(dx) < eps and abs(dy) < eps:
            break
    return x, y

newton(-1, 2, 1e-4)

-1 2
1 2 0
-4 4 1
0.16666666666666666 -0.08333333333333333
-0.8333333333333334 1.9166666666666667
1 2 0.0
-3.3333333333333335 3.8333333333333335 0.06250000000000089
0.011904761904762074 -0.005952380952381037
-0.8214285714285713 1.9107142857142858
1 2 4.440892098500626e-16
-3.285714285714285 3.8214285714285716 0.00031887755102033566
6.136475208622433e-05 -3.068237604333421e-05





(-0.8213672066764851, 1.9106836033382424)

改进的牛顿迭代法 -- 牛顿下山法

依然，我们想求解方程$f(x) = 0$，其中$f(x)$是一个可导函数。我们的目标是找到一个$x_0$，使得$f(x_0) = 0$即可；

由于传统的牛顿迭代法可能迭代到某个点之后就发散了，所以我们牛顿迭代法的基础上，引入了一个步长因子$\lambda$，使得迭代公式变为

$$\begin{align} x_{n+1} = x_n - \lambda \frac{f(x_n)}{f'(x_n)}\\ \end{align}$$

这样，我们就可以在迭代过程中，通过调整步长因子$\lambda$，来控制迭代的速度，从而达到更好的收敛效果。调节步长因子$\lambda$的方法有很多，这里我们采用的是二分法。

def f(x):
    return x**2 - 2

def f1(x):
    return 2*x

def newton_downhill(x0, eps, is_print=True):
    x = x0
    while True:
        l = 1.0
        while abs(f(x - l*f(x)/f1(x))) >= abs(f(x)):
            l /= 2
        x = x - l*f(x)/f1(x)
        if abs(- l*f(x)/f1(x)) < eps:
            break
        if is_print:
            print(x)
    return x

newton2(0.5, 1e-4)

def f(x):
    return x**2 - 2


def f1(x):
    return 2*x


def newton_downhill(x0, eps, is_print=True):
    x = x0
    while True:
        l = 1.0
        while abs(f(x - l*f(x)/f1(x))) >= abs(f(x)):
            l /= 2
        x = x - l*f(x)/f1(x)
        if is_print:
            print(x)
            print(abs(- l*f(x)/f1(x)))
        if abs(- l*f(x)/f1(x)) < eps:
            break
    return x

newton_downhill(0.5, 1e-4)

1.375
0.019886363636363636
1.4147727272727273
0.0005590543994158503
1.4142136728733115
1.1050021215020399e-07





1.4142136728733115

f(0.5)

-1.75

f(2.25)

3.0625

f(1.375)

-0.109375

# 假如用传统方法
def newton(x0, eps,isprint=True):
    x = x0
    while True:
        x = x - f(x)/f1(x)
        if isprint:
            print(x)
            print(abs(f(x)/f1(x)))
        if abs(f(x)/f1(x)) < eps:
            break
    return x

newton(0.5, 1e-4)

2.25
0.6805555555555556
1.5694444444444444
0.14755408062930186
1.4218903638151426
0.007656077875069233
1.4142342859400734
2.072341514131297e-05





1.4142342859400734

小结

这个例子中，用牛顿下山法去做，我们可以看到，迭代的过程中，步长因子$\lambda$是在不断减小的，这样，我们就可以在迭代过程中，通过调整步长因子$\lambda$，来控制迭代的速度，从而达到更好的收敛效果。

但是用普通牛顿迭代法达到同样的精度也只是多一次迭代而已，只能说这个例子可能不是很好，但牛顿下山法确实一定程度上可以防止迭代发散的情况。