求解 R 中的非线性方程组
Solving a system of nonlinear equations in R
假设我有以下方程组:
a * b = 5
sqrt(a * b^2) = 10
如何求解 R 中 a 和 b 的这些方程?
我想这个问题可以表述为优化问题,具有以下函数... ?
fn <- function(a, b) {
rate <- a * b
shape <- sqrt(a * b^2)
return(c(rate, shape) )
}
使用这个库。
library("nleqslv")
您需要定义要求解的多元函数。
fn <- function(x) {
rate <- x[1] * x[2] - 5
shape <- sqrt(x[1] * x[2]^2) - 10
return(c(rate, shape))
}
那你就可以开始了。
nleqslv(c(1,5), fn)
总是看详细的结果。数值计算可能很棘手。在这种情况下,我得到了这个:
Warning message:
In sqrt(x[1] * x[2]^2) : NaNs produced
这只是意味着该程序搜索了一个包含 x[1] < 0
的区域,然后可能会回到飞机的右侧。
在评论中发帖人特别询问了使用 solve
和 optim
所以我们展示了如何解决这个问题 (1) 手动,(2) 使用 solve
,(3 ) 使用 optim
和 (4) 固定点迭代。
1) 手写 首先注意如果我们根据第一个等式写 a = 5/b
并将其代入第二个等式我们得到 sqrt(5/b * b^2) = sqrt(5 * b) = 10
所以b = 20 和 a = 0.25.
2) 求解关于solve
的使用,这些方程可以通过取两边的对数转化为线性形式:
log(a) + log(b) = log(5)
0.5 * (loga + 2 * log(b)) = log(10)
可以表示为:
m <- matrix(c(1, .5, 1, 1), 2)
exp(solve(m, log(c(5, 10))))
## [1] 0.25 20.00
3) optim 使用 optim
我们可以这样写,其中 fn
来自问题。 fn2
是通过减去方程的 RHS 并使用 crossprod
来形成平方和。
fn2 <- function(x) crossprod( fn(x[1], x[2]) - c(5, 10))
optim(c(1, 1), fn2)
给予:
$par
[1] 0.2500805 19.9958117
$value
[1] 5.51508e-07
$counts
function gradient
97 NA
$convergence
[1] 0
$message
NULL
4) fixed point 对于这个,将方程重写为不动点形式,即 c(a, b) = f(c(a, b) ) 然后迭代。一般来说,有几种方法可以做到这一点,但并非所有方法都会收敛,但在这种情况下,这似乎可行。我们对 a
和 b
都使用 1 的起始值,并将第一个方程的两边除以 b
得到第一个方程的定点形式,我们将第二个方程的两边相除通过 sqrt(a)
得到定点形式的第二个方程:
a <- b <- 1 # starting values
for(i in 1:100) {
a = 5 / b
b = 10 / sqrt(a)
}
data.frame(a, b)
## a b
## 1 0.25 20
假设我有以下方程组:
a * b = 5
sqrt(a * b^2) = 10
如何求解 R 中 a 和 b 的这些方程?
我想这个问题可以表述为优化问题,具有以下函数... ?
fn <- function(a, b) {
rate <- a * b
shape <- sqrt(a * b^2)
return(c(rate, shape) )
}
使用这个库。
library("nleqslv")
您需要定义要求解的多元函数。
fn <- function(x) {
rate <- x[1] * x[2] - 5
shape <- sqrt(x[1] * x[2]^2) - 10
return(c(rate, shape))
}
那你就可以开始了。
nleqslv(c(1,5), fn)
总是看详细的结果。数值计算可能很棘手。在这种情况下,我得到了这个:
Warning message:
In sqrt(x[1] * x[2]^2) : NaNs produced
这只是意味着该程序搜索了一个包含 x[1] < 0
的区域,然后可能会回到飞机的右侧。
在评论中发帖人特别询问了使用 solve
和 optim
所以我们展示了如何解决这个问题 (1) 手动,(2) 使用 solve
,(3 ) 使用 optim
和 (4) 固定点迭代。
1) 手写 首先注意如果我们根据第一个等式写 a = 5/b
并将其代入第二个等式我们得到 sqrt(5/b * b^2) = sqrt(5 * b) = 10
所以b = 20 和 a = 0.25.
2) 求解关于solve
的使用,这些方程可以通过取两边的对数转化为线性形式:
log(a) + log(b) = log(5)
0.5 * (loga + 2 * log(b)) = log(10)
可以表示为:
m <- matrix(c(1, .5, 1, 1), 2)
exp(solve(m, log(c(5, 10))))
## [1] 0.25 20.00
3) optim 使用 optim
我们可以这样写,其中 fn
来自问题。 fn2
是通过减去方程的 RHS 并使用 crossprod
来形成平方和。
fn2 <- function(x) crossprod( fn(x[1], x[2]) - c(5, 10))
optim(c(1, 1), fn2)
给予:
$par
[1] 0.2500805 19.9958117
$value
[1] 5.51508e-07
$counts
function gradient
97 NA
$convergence
[1] 0
$message
NULL
4) fixed point 对于这个,将方程重写为不动点形式,即 c(a, b) = f(c(a, b) ) 然后迭代。一般来说,有几种方法可以做到这一点,但并非所有方法都会收敛,但在这种情况下,这似乎可行。我们对 a
和 b
都使用 1 的起始值,并将第一个方程的两边除以 b
得到第一个方程的定点形式,我们将第二个方程的两边相除通过 sqrt(a)
得到定点形式的第二个方程:
a <- b <- 1 # starting values
for(i in 1:100) {
a = 5 / b
b = 10 / sqrt(a)
}
data.frame(a, b)
## a b
## 1 0.25 20