求解非线性方程组 R

Question

我尝试使用 R 中的 nleqslv 函数求解一组非线性方程组。不幸的是，我在猜测正确的初始值以使函数成功运行时遇到了麻烦。我有一个值在 0 和 1 之间的向量，称为 c(t)。它们应满足以下等式

c(t)=A*(exp(-mt)+exp(-m(1024-t)))+B^2

使用 t 的三个后续值，我的目标是使用以下代码确定系数 A、B、m

library(nleqslv)

C10 <- c(1.000000e+00,9.754920e-01,9.547681e-01,9.359057e-01,9.182586e-01,9.014674e-01)
system_size <- 1024

for(i in 2:5)
{
  C <- c(C10[i-1],C10[i],C10[i+1],i-2)
  #function
  target <- function(Coeffs){
  y <- numeric(3)
  y[1] <- Coeffs[1]*(exp(-Coeffs[2]*(C[4]-1))+exp(-Coeffs[2]*(system_size-(C[4]-1))))+Coeffs[3]^2-C[1]
  y[2] <- Coeffs[1]*(exp(-Coeffs[2]*(C[4]))+exp(-Coeffs[2]*(system_size-(C[4]))))+Coeffs[3]^2-C[2]
  y[3] <- Coeffs[1]*(exp(-Coeffs[2]*(C[4]+1))+exp(-Coeffs[2]*(system_size-(C[4]+1))))+Coeffs[3]^2-C[3]
  y
 }
 init <- c(0.001,0.01,0)
 sol <- nleqslv(init, target,control=list(btol=.01), method="Broyden")
 }

使用的初始值反映了我在绘制关联值 c(t) 时得到的结果。尽管如此，生成的输出 sol 给出了

chr "Jacobian is ill-conditioned (1/condition=9.0e-18) (see allowSingular option)"

知道出了什么问题以及如何解决这个问题吗？

---

OP 编辑​​：修改代码以具有最小的工作示例：为 C10 添加了前几个值，调整了循环并为 system_size 添加了值

Answer 1

随着您为C10 添加一些数据，如果考虑到循环应该考虑length(C10)，则该示例运行得非常好。 像这样（有一些变化；为什么见下文）：

library(nleqslv)

C10 <- c(1.000000e+00,9.754920e-01,9.547681e-01,9.359057e-01,9.182586e-01,9.014674e-01)
system_size <- 1024

target <- function(Coeffs){
    y <- numeric(3)
    y[1] <- Coeffs[1]*(exp(-Coeffs[2]*(C[4]-1))+exp(-Coeffs[2]*(system_size-(C[4]-1))))+Coeffs[3]^2-C[1]
    y[2] <- Coeffs[1]*(exp(-Coeffs[2]*(C[4]))+exp(-Coeffs[2]*(system_size-(C[4]))))+Coeffs[3]^2-C[2]
    y[3] <- Coeffs[1]*(exp(-Coeffs[2]*(C[4]+1))+exp(-Coeffs[2]*(system_size-(C[4]+1))))+Coeffs[3]^2-C[3]
    y
}

init <- 50*c(0.001,0.01,0)

for(i in 2:min(length(C10)-1,(system_size/2)))
{

    C <- c(C10[i-1],C10[i],C10[i+1],i-2)
  #function
    target <- function(Coeffs){
        y <- numeric(3)
        y[1] <- Coeffs[1]*(exp(-Coeffs[2]*(C[4]-1))+exp(-Coeffs[2]*(system_size-(C[4]-1))))+Coeffs[3]^2-C[1]
        y[2] <- Coeffs[1]*(exp(-Coeffs[2]*(C[4]))+exp(-Coeffs[2]*(system_size-(C[4]))))+Coeffs[3]^2-C[2]
        y[3] <- Coeffs[1]*(exp(-Coeffs[2]*(C[4]+1))+exp(-Coeffs[2]*(system_size-(C[4]+1))))+Coeffs[3]^2-C[3]
        y
    }

    cat("i=",i,"init=",init, "target(init)=",target(init),"\n")
    sol <- nleqslv(init, target,control=list(btol=.01), method="Broyden")
    print(sol)
}

使用您的初始起始值，模型无法解决并给出您提到的错误消息。我通过增加值更改了init 的值。然后找到解决方案直到i=5。具有给定 C10 的较大值将不会运行，因为在循环内 C10[i+1] 被引用（并且它不存在）。

我在调用nleqslv 之前插入了一个cat 语句，在函数调用之后插入了一个print(sol)，这样至少可以看到发生了什么以及是否真的找到了解决方案。

您不需要指定method="Broyden"，因为它是默认值。

如果发生错误，您应该在for 循环结束出口中测试sol$termcd。

像这样的脚本总是在循环内打印东西！

Mislav 是正确的：起始值可能完全错误。

即使起始值没问题，使用的算法也可能会失败。这就是包提供函数testnslv和searchZeros的原因。

我用testnslv 做了一些实验（此处未显示），结论是method="Newton" 是失败的。狗腿式全球战略似乎总是奏效。线搜索策略并不总是有效。