【问题标题】:Solving a system of nonlinear equations in R求解 R 中的非线性方程组
【发布时间】:2018-02-16 18:19:09
【问题描述】:

假设我有以下方程组:

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) )

}

【问题讨论】:

  • @G5W 你有办法把它变成一个线性方程组吗?
  • 看起来很简单,首先将 b = 5/a 代入第二个等式。
  • 当我搜索你的问题标题时,第一个命中是nleqslv 包的文档。所以我会从那个开始。
  • @G5W 好吧,我在问如何以通用方式在 R 中求解这种方程(例如使用 optimsolve)。感谢@ngm 的nleqslv 建议。我浏览了文档,但不确定如何使用它来解决我的问题...
  • R 不做抽象微积分,也不“解决”任何问题——你可能想在 Mathematica 等中做。

标签: r equation-solving


【解决方案1】:

在评论中,发帖者特别询问了如何使用solveoptim,因此我们展示了如何解决这个问题:(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) 定点 为此,将方程改写为定点形式,即 c(a, b) = f(c(a, b)) 的形式,然后进行迭代。一般来说,有几种方法可以做到这一点,并不是所有的方法都会收敛,但在这种情况下,这似乎是可行的。我们对ab 使用起始值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

【讨论】:

  • 这个答案值得更多的支持。真有见地。非常感谢。
【解决方案2】:

使用这个库。

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] &lt; 0 的区域,然后大概一头扎到了飞机的右侧。

【讨论】:

  • 首先需要通过如下代码安装包:install.packages("nleqslv")。然后你就可以使用这个库了。
猜你喜欢
  • 1970-01-01
  • 2019-03-25
  • 2021-12-13
  • 1970-01-01
  • 2022-07-13
  • 1970-01-01
相关资源
最近更新 更多