Как обойти ошибку «особой матрицы» при численном решении многомерной системы уравнений с использованием R rootSolve :: multiroot - PullRequest
Новая
       13

Как обойти ошибку «особой матрицы» при численном решении многомерной системы уравнений с использованием R rootSolve :: multiroot

0 голосов
/ 26 декабря 2018

Моя система уравнений содержит условие ifelse;Я не знаю, является ли это причиной единственной проблемы с матрицей, потому что она иногда работает, а иногда не работает.

Это работает: 

zf <- function(p,parms) {
  with(as.list(c(p,parms)),{
    eq = NULL
    eq = c(eq,N_A-A-S_AC)
    eq = c(eq,N_C-C-S_AC)
    eq = c(eq,ifelse(S_AC==0,0,A*C-K_AC))
    return(eq)})
}

require(rootSolve)

multiroot(f=zf,start=c(A=0.001,C=0.001,S_AC=0.01),parms=c(N_A=0.0001,N_C=0.01,K_AC=1e-10),positive=TRUE,atol=1.e-50)
#$root
#           A            C         S_AC 
#1.010101e-08 9.900010e-03 9.998990e-05 
#
#$f.root
#[1] 0.000000e+00 2.981556e-19 6.802063e-17
#
#$iter
#[1] 4
#
#$estim.precis
#[1] 2.277293e-17

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

Это, с другой стороны: 

zf2 <- function(p,parms) {
  with(as.list(c(p,parms)),{
    eq = NULL
    eq = c(eq,N_A-A-S_AC-S_AK)
    eq = c(eq,N_C-C-S_AC)
    eq = c(eq,N_K-K-S_AK)
    eq = c(eq,ifelse(S_AC==0,0,A*C-K_AC))
    eq = c(eq,ifelse(S_AK==0,0,A*K-K_AK))
    return(eq)})
}

иногда работает: 

multiroot(f=zf2,start=c(A=0.001,C=0.001,K=0.0005,S_AC=0.01,S_AK=0.005),parms=c(N_A=0.015,N_C=0.01,N_K=0.005,K_AC=1e-10,K_AK=1e-7),positive=TRUE,atol=1.e-50)
#$root
#           A            C            K         S_AC         S_AK 
#3.247940e-04 2.536305e-05 2.994309e-04 9.974637e-03 4.700569e-03 
#
#$f.root
#[1]  0.000000e+00  0.000000e+00 -8.673617e-19  8.137766e-09 -2.746643e-09
#
#$iter
#[1] 4
#
#$estim.precis
#[1] 2.176882e-09

и иногда выдает ошибку:

multiroot(f=zf2,start=c(A=0.001,C=0.001,K=0.0005,S_AC=0.01,S_AK=0.005),parms=c(N_A=0.0001,N_C=0.01,N_K=0.005,K_AC=1e-10,K_AK=1e-7),positive=TRUE,atol=1.e-50)
#diagonal element is zero 
#[1] 4
#$root
#           A            C            K         S_AC         S_AK 
#           0            0 230465362144 230465362144            0 
#
#$f.root
#[1] -2.304654e+11 -2.304654e+11 -2.304654e+11 -1.000000e-10  0.000000e+00
#
#$iter
#[1] 3
#
#$estim.precis
#[1] 138279217286
#
#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

Из других постов по optim и аналогичным функциям я, кажется, понял, что начальные условия (start) могут определять особенности матрицы.
Но тогда я бы не сталзнать, как выбрать start таким образом, чтобы избежать этой проблемы.

Я попытался заменить уравнение ifelse на что-то вроде: 

max(0,A*C-K_AC)

, потому что это уравнение действительно толькоприменимо, когда A*C больше, чем K_AC.
Это, однако, привело к еще более частым ошибкам «единственной матрицы».

У вас есть какие-либо предложения о том, как избежать этой проблемы?
Есть ливы видите что-то не так с моей формулойФункция или что-то еще?

Спасибо!

...