牛顿迭代法与平方根

原理节选自: 牛顿迭代法实现平方根函数

牛顿迭代法:

20210828100504

$x_{n+1}= x_n-f(x_n)/f’(x_n)$

假设输入数是S, 那么要求的平方根为x, 满足$S=x^2$, 那么我们就可以定义函数$f(x)=x^2−S$, 最终问题就转换为求方程$f(x)=0=x^2−S$的根。

$x_{n+1}=x_n-(x_n^2-S)/ 2x_n$

化简之后会变成下面这样:

$x_{n+1}=(x_n^2+S)/2x_n$

函数求导问题

但是我们也可以不化简, 直接让程序帮我求出$f’(x_n)$, 我在接触SICP之前, 甚至都没有想过类似的操作. 现在问题转化为, 当我们拥有了$f(x)$如何求得$f’(x)$, 我们知道在数学上$f’(x)$是这么定义的.

$f’(x) = \frac{f(x+\Delta x) - f(x))}{\Delta x}$

那么根据$f(x)$, 我们可以定义出$f’(x)$, 该函数接受一个值, 并返回该值对应的原函数的导数值. 这是在scheme语言中的定义:

1
2
3
4
5
6
7
8
9
10
11
(define dx 0.00001)
(define DERIV
(lambda (f)
(lambda (X)
(/
(- (f (+ X dx)) (f X))
dx
)
)
)
)

不动点(fixed-point)

我们这里得到的公式是:

$x_{n+1}=x_n-(x_n^2-S)/ 2x_n$

如果把$x_{n+1}$看做y, $x_n$看做x, 那么上面的表达式就是在求下面公式的定点

$f(x) = x - (x^2-S)/2x$

我们想要的结果是多次迭代之后, $x_{n+1}$约等于$x_n$, 因此, 我们正好要求的就是这个迭代函数的定点$f(x)=x$, 可以抽出一个求迭代函数定点的工具, 定义如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
(define (fixed-point f first-guess)
(define tolerance 0.00001)
(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)
)

最终代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
(define (newton f guess)
(define DF (DERIV F))
(fixed-point
(lambda (x) (- x (/ (f x) (DF x))))
guess
)
)
(define (sqrt A)
(newton
(lambda (Y) (- A (square Y)))
1
)
)
(sqrt 1234)

总结

我认为值得借鉴的方案是给出了$f(x)$之后, 能使用(DERIV f)给出$f’(x)$, DERIV接受一个函数, 并返回这个函数的导函数. 另外就是对于执行过程的多次抽象, 比如fixed-pointnewton函数, 它将具体的过程抽象出来, 给了它一个名字. 意味着我们可以使用这些函数去求任意函数的不动点, 或是使用牛顿迭代法求零点.

附: JS实现

感觉要写很多return语句, 可能没有原来那么优雅.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
const dx = 0.00001;
function deriv(f) {
return (x) => {
return (f(x+dx)-f(x))/dx;
}
}

function fixed_point(f, first_guess) {
const tolerance = 0.00001;
const close_enough = function(v1, v2) {
return Math.abs(v1-v2) < tolerance;
};
const try_guess = function(guess) {
const next_guess = f(guess);
if (close_enough(guess, next_guess)) {
return next_guess;
}
return try_guess(next_guess);
};
return try_guess(first_guess);
}

function newton(f, guess) {
const DF = deriv(f);

return fixed_point(
(x) => {return x - f(x)/DF(x)},
guess,
)
}

function sqrt(A) {
return newton(
(x) => (A - x*x),
1
)
}
sqrt(1234);

// 35.12833614050059