牛顿迭代法与平方根
原理节选自: 牛顿迭代法实现平方根函数
牛顿迭代法:
$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-point与newton函数,
它将具体的过程抽象出来, 给了它一个名字. 意味着我们可以使用这些函数去求任意函数的不动点,
或是使用牛顿迭代法求零点.
附: 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);
|