文章目录
- Newton求根法
- 算法
- 求解平方根
Newton求根法
1664年Newton提出了一种迭代求根的方法。这种方法有时也被叫做Raphson方法。虽然Newton更早发现了这种方法,但Raphson首先在自己的文章中介绍了这种方法。
该方法解决的问题如下。
给定一个方程:
f
(
x
)
=
0
f(x) = 0
f(x)=0
求解该方程组。更具体地来说,我们需要找到这个方程的一个根(假设这个根存在)。
f
(
x
)
f(x)
f(x)在
[
a
,
b
]
[a,b]
[a,b]上是可导的。
算法
算法的输入包括函数 f ( x ) f(x) f(x)和一个初始的近似值 x 0 x_0 x0,这个 x 0 x_0 x0就是算法开始的地方。
假设 x i x_i xi已知,那么需要按照接下来的方法求 x i + 1 x_{i+1} xi+1。画出 f ( x ) f(x) f(x)在点 x = x i x = x_i x=xi处的切线,这条线与 x x x轴的交点,这个交点的横坐标就是 x i + 1 x_{i+1} xi+1,然后重复整个过程。
不难得到下面的这个等式:
x
i
+
1
=
x
i
−
f
(
x
i
)
f
′
(
x
i
)
x_{i+1} = x_i - \frac{f(x_i)}{f^\prime(x_i)}
xi+1=xi−f′(xi)f(xi)
首先,计算出
f
(
x
)
f(x)
f(x)的导数,斜率
f
′
(
x
)
f'(x)
f′(x),确定切线的方程为:
y
−
f
(
x
i
)
=
f
′
(
x
i
)
(
x
i
+
1
−
x
i
)
y - f(x_i) = f'(x_i)(x_{i+1} - x_i)
y−f(xi)=f′(xi)(xi+1−xi)
通过解这个方程可以得到
x
i
+
1
x_{i+1}
xi+1的值。
注意到,如果方程 f ( x ) f(x) f(x)是平滑的,并且 x i x_i xi离根足够接近,那么 x i + 1 x_{i+1} xi+1与所求的根更加接近。
收敛速率是二次的。
求解平方根
求解平方根是应用Newton方法的一个例子。
求解
n
\sqrt{n}
n即求方程
f
(
x
)
=
x
2
−
n
f(x) = x^2 - n
f(x)=x2−n的根,带入上面的式子化简,可以得到
x
i
+
1
=
x
i
+
n
x
i
2
x_{i+1} = \frac{x_i + \frac{n}{x_i}}{2}
xi+1=2xi+xin
这个问题的第一个形式时,当
n
n
n是一个有理数时,给定一个eps
,可以求出它的根:
double sqrt_newton(double n) {
const double eps = 1E-15;
double x = 1;
for (;;) {
double nx = (x + n / x) / 2;
if (abs(x - nx) < eps)
break;
x = nx;
}
return x;
}
另一个常见的问题形式是要求求出它的整数根(给定一个 n n n,求出最大的 x x x使得 x 2 ≤ n x^2 \leq n x2≤n)。这样的化就需要稍微改变一下这个算法的结束条件,因为当 x x x开始接近答案时,可能会“跳动”。因此需要添加一个条件,如果 x x x在上一次迭代中减小了,在当前的迭代中又增加了,这个时候算法必须停止。
int isqrt_newton(int n) {
int x = 1;
bool decreased = false;
for (;;) {
int nx = (x + n / x) >> 1;
if (x == nx || nx > x && decreased)
break;
decreased = nx < x;
x = nx;
}
return x;
}
最后,这个问题还有一个解决非常大的数字(bignum arithmetic)的形式。由于数字 n n n可能会非常大,所以算法开始时有理由对答案进行一个初始的估计。很显然,这个估计值离答案跃进,算法收敛的越快。一般我们会选择数字 2 b i t s / 2 2^{bits/2} 2bits/2作为这个估计值,其中 b i t s bits bits是 n n n的二进制位数。下面是这个问题的java实现。
public static BigInteger isqrtNewton(BigInteger n) {
BigInteger a = BigInteger.ONE.shiftLeft(n.bitLength() / 2);
boolean p_dec = false;
for (;;) {
BigInteger b = n.divide(a).add(a).shiftRight(1);
if (a.compareTo(b) == 0 || a.compareTo(b) < 0 && p_dec)
break;
p_dec = a.compareTo(b) > 0;
a = b;
}
return a;
}