如何推导“求算数平方根”的递推公式
如何推导“求算数平方根”的递推公式
前言
回头一看已经有将近四年没有写博文了，自从当了社畜就积极性不高了，最近恰逢裁员广进，想着学学go语言，在Tour of Go中看到了这样（z -= (z * z - x) / (2 * z)）的求算数平方根的递推公式，所以想出个推导过程，过程很简单，不需要知道很深的数学知识，知道导数、切线方程、数列足矣。
Tour of Go 平方根递推式地址https://tour.go-zh.org/flowcontrol/8
原文提到这个公式使用牛顿法，感兴趣的同学可以去看一看。
那么我们开始我们的推导过程。
首先，要求一个数的平方根，假定这个数为
a
(
a
>
=
0
)
\ a(a>=0)
a(a>=0)，我们要求得一个
x
\ x
x，使得
x
2
=
a
\ x^2=a
x2=a成立。
那么我们设方程
f
(
x
)
=
x
2
−
a
\ f(x)=x^2-a
f(x)=x2−a，这个一元二次方程图像为：
图中绿点就是我们最终要求的值，那么我们找一个点
P
(
X
n
,
f
(
X
n
)
)
\ P (X_{n}, f(X_{n}))
P(Xn,f(Xn))，做一条切线，切线与x轴交于点
Q
(
X
n
+
1
,
0
)
\ Q(X_{n+1},0)
Q(Xn+1,0)，如下图：
设切线方程为
y
=
k
x
+
b
\ y = kx + b
y=kx+b，易得
k
=
f
′
(
X
n
)
=
2
X
n
\ k=f'(X_n)=2X_n
k=f′(Xn)=2Xn，带入P、Q两点坐标，得到
f
(
X
n
)
=
2
X
n
2
+
b
\ f(X_n)=2X_n^2+b
f(Xn)=2Xn2+b ①
0
=
2
X
n
X
n
+
1
+
b
\ 0 = 2X_nX_{n+1}+b
0=2XnXn+1+b ②
①-②得
f
(
X
n
)
=
2
X
n
2
−
2
X
n
X
n
+
1
\ f(X_n) =2X_n^2- 2X_nX_{n+1}
f(Xn)=2Xn2−2XnXn+1，
又
f
(
X
n
)
=
X
n
2
−
a
\ f(X_n)=X_n^2-a
f(Xn)=Xn2−a,
整理得
X
n
+
1
=
1
2
(
X
n
+
a
X
n
)
\ X_{n+1}=\frac{1}{2}(X_n+\frac{a}{X_n})
Xn+1=21(Xn+Xna)
由图像(紫色线)我们可以看到，当
P
\ P
P趋近于
X
\ X
X时，对于
P
\ P
P点做切线交于
x
\ x
x轴的
Q
\ Q
Q点将近似于
X
\ X
X
那么说
X
n
+
1
=
1
2
(
X
n
+
a
X
n
)
\ X_{n+1}=\frac{1}{2}(X_n+\frac{a}{X_n})
Xn+1=21(Xn+Xna)在
X
n
+
1
2
−
a
≈
0
\ X_{n+1}^2-a\approx0
Xn+12−a≈0时，
X
n
+
1
≈
a
\ X_{n+1}\approx\sqrt{a}
Xn+1≈a
此时，递推公式已经推出来了，那么有小伙伴说，俩公式不一样啊，其实Tour of Go中和这个是一个递推式，请看
X
n
+
1
=
1
2
(
X
n
+
a
X
n
)
=
1
2
X
n
+
a
2
X
n
=
X
n
−
1
2
X
n
+
a
2
X
n
=
X
n
−
X
n
2
−
a
2
X
n
\ X_{n+1}=\frac{1}{2}(X_n+\frac{a}{X_n})=\frac{1}{2}X_n+\frac{a}{2X_n}=Xn-\frac{1}{2}X_n+\frac{a}{2X_n}=X_n-\frac{X_n^2-a}{2X_n}
Xn+1=21(Xn+Xna)=21Xn+2Xna=Xn−21Xn+2Xna=Xn−2XnXn2−a
这样和Tour of Go中提到的形式就一模一样了！
用Go把递推公式实现一下吧
X
n
+
1
=
X
n
−
X
n
2
−
a
2
X
n
\ X_{n+1}=X_n-\frac{X_n^2-a}{2X_n}
Xn+1=Xn−2XnXn2−a
// 递推求x的算数平方根z
x := 100.0
z := x / 2
for {
t := z*z - x
if t < 0.0000001 && t >= 0 {
break
}
z -= t / (2 * z)
}
fmt.Println(z)
执行结果为10.000000000107446
X n + 1 = 1 2 ( X n + a X n ) \ X_{n+1}=\frac{1}{2}(X_n+\frac{a}{X_n}) Xn+1=21(Xn+Xna)
// 递推求x的算数平方根z
x := 100.0
z := x / 2
for {
z = 0.5 * (z + x / z)
if z*z - x <= 0.0000001 {
break
}
}
fmt.Println(z)
执行结果为10.000000000107445