Решите систему нелинейных уравнений
Я пытаюсь решить следующую систему из четырех уравнений. Я пробовал использовать пакет "rootSolve", но, похоже, я не могу найти решение таким образом.
Я использую следующий код:
model <- function(x) {
F1 <- sqrt(x[1]^2 + x[3]^2) -1
F2 <- sqrt(x[2]^2 + x[4]^2) -1
F3 <- x[1]*x[2] + x[3]*x[4]
F4 <- -0.58*x[2] - 0.19*x[3]
c(F1 = F1, F2 = F2, F3 = F3, F4 = F4)
}
(ss <- multiroot(f = model, start = c(0,0,0,0)))
Но это дает мне следующую ошибку:
Warning messages:
1: In stode(y, times, func, parms = parms, ...) :
error during factorisation of matrix (dgefa); singular matrix
2: In stode(y, times, func, parms = parms, ...) : steady-state not reached
Я изменил начальные значения, как предлагается в другом аналогичном ответе, и для некоторых я могу найти решение. Однако эта система - согласно источнику, который я использую - должна иметь однозначно идентифицированное решение. Есть идеи, как решить эту систему?
Спасибо!
Ответы 2
Вышеупомянутое предупреждение гласит, что при использовании начального значения, которое вы указали для multiroot, он не смог найти оптимального решения.
Давай попробуем это -
library(rootSolve)
model <- function(x) {
F1 <- sqrt(x[1]^2 + x[3]^2) - 1
F2 <- sqrt(x[2]^2 + x[4]^2) - 1
F3 <- x[1]*x[2] + x[3]*x[4]
F4 <- -0.58*x[2] - 0.19*x[3]
c(F1 = F1, F2 = F2, F3 = F3, F4 = F4)
}
#solution
(ss <- multiroot(f = model, start = c(1.5, 0, 0.5, 0)))
это дает
> ss
$root
[1] 1.000000e+00 4.752703e-12 -1.450825e-11 1.000000e+00
$f.root
F1 F2 F3 F4
3.404610e-12 3.494982e-13 -9.755549e-12 1.929753e-20
$iter
[1] 7
$estim.precis
[1] 3.377414e-12
После нескольких испытаний я заметил, что всякий раз, когда я меняю его начальное значение, я получаю почти один и тот же результат (например, 1, 0, 0, 1) каждый раз.
Ваша система уравнений имеет несколько решений.
Я использую другой пакет для решения вашей системы: nleqslv следующим образом:
library(nleqslv)
model <- function(x) {
F1 <- sqrt(x[1]^2 + x[3]^2) - 1
F2 <- sqrt(x[2]^2 + x[4]^2) - 1
F3 <- x[1]*x[2] + x[3]*x[4]
F4 <- -0.58*x[2] - 0.19*x[3]
c(F1 = F1, F2 = F2, F3 = F3, F4 = F4)
}
#find solution
xstart <- c(1.5, 0, 0.5, 0)
nleqslv(xstart,model)
Это дает то же решение, что и ответ Прем.
Однако в вашей системе есть несколько решений.
Пакет nleqslv предоставляет функцию поиска решений по матрице различных начальных значений. Вы можете использовать это
set.seed(13)
xstart <- matrix(runif (400,0,2),ncol=4)
searchZeros(xstart,model)
(Примечание: разные семена могут не найти все четыре решения)
Вы увидите, что есть четыре разных решения:
$x
[,1] [,2] [,3] [,4]
[1,] -1 -1.869055e-10 5.705536e-10 -1
[2,] -1 4.992198e-13 -1.523934e-12 1
[3,] 1 -1.691309e-10 5.162942e-10 -1
[4,] 1 1.791944e-09 -5.470144e-09 1
.......
Это ясно указывает на то, что точные решения представлены в следующей матрице
xsol <- matrix(c(1,0,0,1,
1,0,0,-1,
-1,0,0,1,
-1,0,0,-1),byrow=TRUE,ncol=4)
А потом сделай
model(xsol[1,])
model(xsol[2,])
model(xsol[3,])
model(xsol[4,])
Подтвержденный!
Я не пытался найти эти решения аналитически, но вы можете видеть, что если x[2] и x[3] равны нулю, то F3 и F4 равны нулю. После этого можно сразу найти решения для x[1] и x[4].
Ух ты! Итак, наконец, у нас есть действительно хороший пакет
nleqslv(+1 за ответ).