nleqslv, R, нелинейная система уравнений - PullRequest
Новая
       4

nleqslv, R, нелинейная система уравнений

0 голосов
/ 29 апреля 2018

У меня есть следующая система нелинейных уравнений: 

r1*r2*(0.25*0.6061-0.5)+r1+(0.25*r2) = 0.6061
r1*r2*(0.25*0.6429-0.5)+r1+(0.25*r2) = 0.6429

Мне нужно найти r1 и r2. Их решение должно быть равно: 

r1 = 0.5754
r2 = 0.6201

, полученный в результате Фортрана.
Я не знаю, как найти решения, используя пакет nleqslv.

Буду признателен, если кто-нибудь сможет мне помочь.

Спасибо.

1 Ответ

0 голосов
/ 29 апреля 2018

Способ использования nleqslv заключается в определении функции, возвращающей 2-мерный вектор, как c(r1, r2).
Эта функция принимает 2-мерный аргумент x. Я заменил r1 и r2 на x[1] и x[2] соответственно.
Я также заменил вашу 1-ую константу уравнения (0.25*0.6061-0.5) на ее значение -0.348475. А во 2-м уравнении (0.25*0.6429-0.5) на -0.339275.
Что касается начальных значений, я сначала попробовал c(0, 0), но это не сработало. Если вы запустите приведенный ниже код с начальными значениями c(0.5, 0.5), вы получите то же решение с точностью до плавающей запятой. 

fn <- function(x){
    c(-0.348475*x[1]*x[2] + x[1] + 0.25*x[2] - 0.6061,
      -0.339275*x[1]*x[2] + x[1] + 0.25*x[2] - 0.6429)
}

nleqslv(c(1, 1), fn)
#$`x`
#[1] 1.000000 3.999999
#
#$fvec
#[1] 4.773959e-14 4.607426e-14
#
#$termcd
#[1] 1
#
#$message
#[1] "Function criterion near zero"
#
#$scalex
#[1] 1 1
#
#$nfcnt
#[1] 2
#
#$njcnt
#[1] 1
#
#$iter
#[1] 2

Обратите внимание, что $fvec близок к нулю, а $termcd равен 1. Это означает, что алгоритм сходится. Чтобы получить только решение, вы можете запустить 

sol <- nleqslv(c(1, 1), fn)
sol$x
#[1] 1.000000 3.999999

Если вы можете вычислить якобиан, и в этом случае это даже очень просто, алгоритм будет более точным. 

jacfn <- function(x) {
    n <- length(x)
    Df <- matrix(numeric(n*n), n, n)
    Df[1,1] <- -0.348475*x[2] + 1
    Df[1,2] <- -0.348475*x[1] + 0.25
    Df[2,1] <- -0.339275*x[2] + 1
    Df[2,2] <- -0.339275*x[1] + 0.25

    Df
}

sol2 <- nleqslv(c(1, 1), fn, jac = jacfn)
sol2$x
#[1] 1 4

Решение такое же.

...