程振兴 @czxa.top
截止今天,我已经在本博客上写了604.4k个字了!
引言
$\sqrt{2}$的求解使用软件是非常简单的事情,但是其求解过程是如何实现的呢?本文介绍了3种简单的实现方法。
最优化方法
这种方法最为简单。下面是推导过程:
$$
\sqrt{2} = x \\
x^2 = 2 \\
x^2 - 2 = 0 \\
So, we\ just\ need\ to\ minimize{(x^2-2)^2}
$$
代码如下,R的optimize()
函数可以实现最优化,一个函数对象、一个初始求解区间即可:1
2
3
4obj <- function(x){
return((x^2-2)^2)
}
optimize(obj, interval = c(0, 2))
牛顿迭代法
牛顿迭代法的思想如下图:1
2
3
4
5
6
7
8
9
10
11
12
13
14x <- seq(1, 10, 0.1)
y <- x^2 - 2
plot(x, y, type = "l", ann = F)
arrows(8, 0, 8, 62)
arrows(8, 62, 4.125, 0)
arrows(4.125, 0, 4.125, 15.0156)
arrows(4.125, 15.0156, 2.305, 0)
y1 <- rep(0,91)
lines(x, y1, lty = 1)
points(sqrt(2), 0, pch = 9)
text(7, 70, expression(x[2] == x[1] - frac((x[1]^2 - 2), (2%.%x[1]))))
title("牛顿迭代法", family="STSongti-SC-Regular", xlab = 'x', ylab = 'y')
axis(1, family="STSongti-SC-Regular")
axis(2, family="STSongti-SC-Regular")
其中,$x_1$和$x_2$的关系推导方法如下:
$$
\frac{x_1^2-2}{x_1-x_2} = x_2 \\
x_2 = x_1 - \frac{x_1^2-2}{2x_1}
$$
用R实现这一过程的代码为:1
2
3
4
5
6
7
8obj1 <- function(x){
return(x^2-2)
}
x <- 0.1
while(abs(obj1(x)) > 0.0000000001){
x <- x - (x^2-2)/(2*x)
}
print(x)
丁文亮同学还给出了一段看似不同,实际相同的代码:1
2
3
4
5
6
7
8
9
10
11options(digits = 22)
i <- 0
A <- 2
x <- 2
result <- c()
while(i < 5) {
x <- 0.5*(x + A/x)
i <- i + 1
result <- c(result, x)
}
cbind(result, error = (result-sqrt(2)))
二分法
这个方法应该是最早认识的一种方法了,记得小学的时候就用手算过$\sqrt{3}$。就是先选定一个恰当的区间,然后把中间值代进入,和0比较,如果大于0,取左边的子区间(假设左边界函数值小于0,右边界函数值大于0),反之则取右边的子区间。R的实现代码如下:1
2
3
4
5
6
7
8
9
10
11
12
13
14obj1 <- function(x){
return(x^2-2)
}
y <- c(1, 2)
while(abs(obj1((y[1]+y[2])/2)) > 0.00000000001){
ymean <- obj1((y[1]+y[2])/2)
if (ymean > 0) {
y <- c(y[1], (y[1]+y[2])/2)
}
else {
y <- c((y[1]+y[2])/2, y[2])
}
}
print((y[1]+y[2])/2)
小结
实际上这个问题的解决有两个过程,第一是寻找迭代关系,第二是寻找能使得结果收敛的尽可能快的迭代关系。