牛顿迭代法

Author : zbzhen,        Modified : Wed Dec 25 19:43:47 2024

1. 基本理论

根据微分近似公式

f(x)f(x0)+(xx0)f(x0)f(x)\approx f(x_0)+(x-x_0)f'(x_0)

如果f(x)=0f(x)=0,于是有

xx0f(x0)f(x0)x \approx x_0-\dfrac{f(x_0)}{f'(x_0)}

从而推得迭代公式

x(n+1)=x(n)f(x(n))f(x(n))x^{(n+1)}=x^{(n)}-\dfrac{f(x^{(n)})}{f'(x^{(n)})}

需要注意两点:

  • 迭代公式的分母不能为零
  • ff是多元函数, 则f(x(n))f'(x^{(n)})是一个矩阵, 通常记为f\nabla f

求梯度的简单理解

主要难点是多元问题, 除了可以简单的通过定义, 大部分情况下可以如下考虑:

直接通过例子说明:

f(x,y)=[2x+3xy+4yx2+y2sin(xy)]f(x,y) = \left[ \begin{array}{c} 2x+ 3xy + 4y \\ x^2 + y^2 \\ \sin(xy) \end{array} \right]

f(x,y)=[2+3y3x+42x2yysin(xy)xsin(xy)]\nabla f(x,y) = \left[ \begin{array}{c} 2 + 3y & 3x + 4 \\ 2x & 2y\\ y\sin(xy) &x\sin(xy)\\ \end{array} \right]

通常并不会显示给出f\nabla f, 而是直接得到它与向量的乘法

f[uv]=[2+3y3x+42x2yysin(xy)xsin(xy)][uv]=[2u+3yu+3xv+4v2ux+2vyuysin(xy)+vxsin(xy)]\nabla f \begin{bmatrix} u\\ v \end{bmatrix} = \left[ \begin{array}{c} 2 + 3y & 3x + 4 \\ 2x & 2y\\ y\sin(xy) &x\sin(xy)\\ \end{array} \right] \begin{bmatrix} u\\ v \end{bmatrix} = \begin{bmatrix} 2u+3yu+3xv + 4v \\ 2ux + 2vy \\ uy\sin(xy) +vx\sin(xy)\\ \end{bmatrix}

2. 问题

2.1. 问题一

牛顿迭代求一元函数f(x)=x22f(x)=x^2-2的一个正零点

2.2. python代码

def f(x):
    return abs(x*x - 2)
    
def df(x, h):
    return (f(x+h)-f(x))/h

def ite(x, h):
    return x - f(x)/df(x, h)

x = 1
h = 0.001
for i in range(10):
    x = ite(x, h) 
print(x)

2.3. 问题二

牛顿迭代求多元函数f(x)=[5x1+34x022sin(x1x2)x1x21.5]\bm f(\bm x)= \begin{bmatrix} 5x_1+3\\ 4x_0^2-2\sin(x_1x_2)\\ x_1x_2-1.5 \end{bmatrix}的一个零点

import numpy as np
def f(x):
    x0, x1, x2 = x
    ans= np.array([5*x1 + 3, 4*x0**2 - 2*np.sin(x1*x2), x1*x2 - 1.5])
    return ans

def df(x, h=0.00001):
    x = x[:, None]
    hmpx = x + np.eye(len(x))*h
    ans = np.array(f(hmpx) - f(x))/h
    return ans

def ite(x, h):
    return x - np.linalg.solve(df(x, h), f(x))

x = np.array([2, 2, -2])
h = 0.000001
for i in range(10):
    x = ite(x, h) 
print(x)
print(f(x))