数值求根:R 中的二分法
Numerical Root Finding: Bisection Method in R
我需要使用二分法求数值根,并打印每次迭代涉及的变量的值,直到达到某个值。
bisection <- function(x1, x2){
l <- vector(mode="integer")
l[1] <- x1
r <- vector(mode="integer")
r[1] <- x2
m <- vector(mode="integer")
gl <- vector(mode="integer")
gr <- vector(mode="integer")
gm <- vector(mode="integer")
root <- 5e-8
i <- 1
repeat{
m[i] <- (l[i]+r[i])/2
gl[i] <- gx(l[i])
gr[i] <- gx(r[i])
gm[i] <- gx(m[i])
if (isTRUE(abs(gm[i]) > root) && isTRUE(gl[i]*gm[i] < 0)){
l[i+1] <- l[i]
r[i+1] <- m[i]
}
if (isTRUE(gm[i] > root) && isTRUE(gr[i]*gm[i] < 0)){
l[i+1] <- m[i]
r[i+1] <- r[i]
}
else if (isTRUE(abs(gm[i]) <= root)){
j <- c(0:(length(gm)-1))
df <- data.frame(j, l,r,m,gl,gr,gm)
names(df) <- c("i", "xl","xr","xm", "gxl","gxr", "gxm")
print(df)
break
}
}
}
当我用 bisection(1,1.5)
尝试 运行 这个函数时,它的输出只有一行迭代,即使手动解决它也会导致至少 12 次迭代。它也挂起(?)。
我不知道我哪里错了。请帮忙。
编辑说 gx 函数是这样的:gx <- function(x){x^3-x-1}
您的代码有两个问题:
- 在你的第二个
if
语句中你忘记了 abs
,它必须是 if (isTRUE(abs(gm[i]) > root)
- 您忘记增加计数器的值,即
i
仍为 1。
此外,我对您的代码进行了小幅调整。不要使用 if-if-else-if,首先检查是否找到了根,如果没有,则使用一个 if-else 更新 l 和 r。使其更易于阅读和理解。
bisection <- function(x1, x2){
l <- vector(mode="integer")
l[1] <- x1
r <- vector(mode="integer")
r[1] <- x2
m <- vector(mode="integer")
gl <- vector(mode="integer")
gr <- vector(mode="integer")
gm <- vector(mode="integer")
root <- 5e-8
i <- 1
repeat{
m[i] <- (l[i]+r[i])/2
gl[i] <- gx(l[i])
gr[i] <- gx(r[i])
gm[i] <- gx(m[i])
if (isTRUE(abs(gm[i]) <= root)){
j <- c(0:(length(gm)-1))
df <- data.frame(j, l,r,m,gl,gr,gm)
names(df) <- c("i", "xl","xr","xm", "gxl","gxr", "gxm")
print(df)
break
}
if (isTRUE(abs(gm[i]) > root) && isTRUE(gl[i]*gm[i] < 0)){
l[i+1] <- l[i]
r[i+1] <- m[i]
} else {
l[i+1] <- m[i]
r[i+1] <- r[i]
}
i <- i +1
}
}
gx <- function(x) x^2 - 2
bisection(1, 1.5)
#> i xl xr xm gxl gxr gxm
#> 1 0 1.000000 1.500000 1.250000 -1.000000e+00 2.500000e-01 -4.375000e-01
#> 2 1 1.250000 1.500000 1.375000 -4.375000e-01 2.500000e-01 -1.093750e-01
#> 3 2 1.375000 1.500000 1.437500 -1.093750e-01 2.500000e-01 6.640625e-02
#> 4 3 1.375000 1.437500 1.406250 -1.093750e-01 6.640625e-02 -2.246094e-02
#> 5 4 1.406250 1.437500 1.421875 -2.246094e-02 6.640625e-02 2.172852e-02
#> 6 5 1.406250 1.421875 1.414062 -2.246094e-02 2.172852e-02 -4.272461e-04
#> 7 6 1.414062 1.421875 1.417969 -4.272461e-04 2.172852e-02 1.063538e-02
#> 8 7 1.414062 1.417969 1.416016 -4.272461e-04 1.063538e-02 5.100250e-03
#> 9 8 1.414062 1.416016 1.415039 -4.272461e-04 5.100250e-03 2.335548e-03
#> 10 9 1.414062 1.415039 1.414551 -4.272461e-04 2.335548e-03 9.539127e-04
#> 11 10 1.414062 1.414551 1.414307 -4.272461e-04 9.539127e-04 2.632737e-04
#> 12 11 1.414062 1.414307 1.414185 -4.272461e-04 2.632737e-04 -8.200109e-05
#> 13 12 1.414185 1.414307 1.414246 -8.200109e-05 2.632737e-04 9.063259e-05
#> 14 13 1.414185 1.414246 1.414215 -8.200109e-05 9.063259e-05 4.314817e-06
#> 15 14 1.414185 1.414215 1.414200 -8.200109e-05 4.314817e-06 -3.884337e-05
#> 16 15 1.414200 1.414215 1.414207 -3.884337e-05 4.314817e-06 -1.726433e-05
#> 17 16 1.414207 1.414215 1.414211 -1.726433e-05 4.314817e-06 -6.474773e-06
#> 18 17 1.414211 1.414215 1.414213 -6.474773e-06 4.314817e-06 -1.079981e-06
#> 19 18 1.414213 1.414215 1.414214 -1.079981e-06 4.314817e-06 1.617417e-06
#> 20 19 1.414213 1.414214 1.414214 -1.079981e-06 1.617417e-06 2.687177e-07
#> 21 20 1.414213 1.414214 1.414213 -1.079981e-06 2.687177e-07 -4.056319e-07
#> 22 21 1.414213 1.414214 1.414214 -4.056319e-07 2.687177e-07 -6.845708e-08
#> 23 22 1.414214 1.414214 1.414214 -6.845708e-08 2.687177e-07 1.001303e-07
#> 24 23 1.414214 1.414214 1.414214 -6.845708e-08 1.001303e-07 1.583661e-08
我需要使用二分法求数值根,并打印每次迭代涉及的变量的值,直到达到某个值。
bisection <- function(x1, x2){
l <- vector(mode="integer")
l[1] <- x1
r <- vector(mode="integer")
r[1] <- x2
m <- vector(mode="integer")
gl <- vector(mode="integer")
gr <- vector(mode="integer")
gm <- vector(mode="integer")
root <- 5e-8
i <- 1
repeat{
m[i] <- (l[i]+r[i])/2
gl[i] <- gx(l[i])
gr[i] <- gx(r[i])
gm[i] <- gx(m[i])
if (isTRUE(abs(gm[i]) > root) && isTRUE(gl[i]*gm[i] < 0)){
l[i+1] <- l[i]
r[i+1] <- m[i]
}
if (isTRUE(gm[i] > root) && isTRUE(gr[i]*gm[i] < 0)){
l[i+1] <- m[i]
r[i+1] <- r[i]
}
else if (isTRUE(abs(gm[i]) <= root)){
j <- c(0:(length(gm)-1))
df <- data.frame(j, l,r,m,gl,gr,gm)
names(df) <- c("i", "xl","xr","xm", "gxl","gxr", "gxm")
print(df)
break
}
}
}
当我用 bisection(1,1.5)
尝试 运行 这个函数时,它的输出只有一行迭代,即使手动解决它也会导致至少 12 次迭代。它也挂起(?)。
我不知道我哪里错了。请帮忙。
编辑说 gx 函数是这样的:gx <- function(x){x^3-x-1}
您的代码有两个问题:
- 在你的第二个
if
语句中你忘记了abs
,它必须是if (isTRUE(abs(gm[i]) > root)
- 您忘记增加计数器的值,即
i
仍为 1。
此外,我对您的代码进行了小幅调整。不要使用 if-if-else-if,首先检查是否找到了根,如果没有,则使用一个 if-else 更新 l 和 r。使其更易于阅读和理解。
bisection <- function(x1, x2){
l <- vector(mode="integer")
l[1] <- x1
r <- vector(mode="integer")
r[1] <- x2
m <- vector(mode="integer")
gl <- vector(mode="integer")
gr <- vector(mode="integer")
gm <- vector(mode="integer")
root <- 5e-8
i <- 1
repeat{
m[i] <- (l[i]+r[i])/2
gl[i] <- gx(l[i])
gr[i] <- gx(r[i])
gm[i] <- gx(m[i])
if (isTRUE(abs(gm[i]) <= root)){
j <- c(0:(length(gm)-1))
df <- data.frame(j, l,r,m,gl,gr,gm)
names(df) <- c("i", "xl","xr","xm", "gxl","gxr", "gxm")
print(df)
break
}
if (isTRUE(abs(gm[i]) > root) && isTRUE(gl[i]*gm[i] < 0)){
l[i+1] <- l[i]
r[i+1] <- m[i]
} else {
l[i+1] <- m[i]
r[i+1] <- r[i]
}
i <- i +1
}
}
gx <- function(x) x^2 - 2
bisection(1, 1.5)
#> i xl xr xm gxl gxr gxm
#> 1 0 1.000000 1.500000 1.250000 -1.000000e+00 2.500000e-01 -4.375000e-01
#> 2 1 1.250000 1.500000 1.375000 -4.375000e-01 2.500000e-01 -1.093750e-01
#> 3 2 1.375000 1.500000 1.437500 -1.093750e-01 2.500000e-01 6.640625e-02
#> 4 3 1.375000 1.437500 1.406250 -1.093750e-01 6.640625e-02 -2.246094e-02
#> 5 4 1.406250 1.437500 1.421875 -2.246094e-02 6.640625e-02 2.172852e-02
#> 6 5 1.406250 1.421875 1.414062 -2.246094e-02 2.172852e-02 -4.272461e-04
#> 7 6 1.414062 1.421875 1.417969 -4.272461e-04 2.172852e-02 1.063538e-02
#> 8 7 1.414062 1.417969 1.416016 -4.272461e-04 1.063538e-02 5.100250e-03
#> 9 8 1.414062 1.416016 1.415039 -4.272461e-04 5.100250e-03 2.335548e-03
#> 10 9 1.414062 1.415039 1.414551 -4.272461e-04 2.335548e-03 9.539127e-04
#> 11 10 1.414062 1.414551 1.414307 -4.272461e-04 9.539127e-04 2.632737e-04
#> 12 11 1.414062 1.414307 1.414185 -4.272461e-04 2.632737e-04 -8.200109e-05
#> 13 12 1.414185 1.414307 1.414246 -8.200109e-05 2.632737e-04 9.063259e-05
#> 14 13 1.414185 1.414246 1.414215 -8.200109e-05 9.063259e-05 4.314817e-06
#> 15 14 1.414185 1.414215 1.414200 -8.200109e-05 4.314817e-06 -3.884337e-05
#> 16 15 1.414200 1.414215 1.414207 -3.884337e-05 4.314817e-06 -1.726433e-05
#> 17 16 1.414207 1.414215 1.414211 -1.726433e-05 4.314817e-06 -6.474773e-06
#> 18 17 1.414211 1.414215 1.414213 -6.474773e-06 4.314817e-06 -1.079981e-06
#> 19 18 1.414213 1.414215 1.414214 -1.079981e-06 4.314817e-06 1.617417e-06
#> 20 19 1.414213 1.414214 1.414214 -1.079981e-06 1.617417e-06 2.687177e-07
#> 21 20 1.414213 1.414214 1.414213 -1.079981e-06 2.687177e-07 -4.056319e-07
#> 22 21 1.414213 1.414214 1.414214 -4.056319e-07 2.687177e-07 -6.845708e-08
#> 23 22 1.414214 1.414214 1.414214 -6.845708e-08 2.687177e-07 1.001303e-07
#> 24 23 1.414214 1.414214 1.414214 -6.845708e-08 1.001303e-07 1.583661e-08