Hello,
I am trying to solve a system of integral equations using multiroot. I have
tried asking on stack exchange and reddit without any luck.
Multiroot uses the library(RootSolve).
I have two integral equations involving constants S[1] and S[2] (which are
free.) I would like to find what *positive* values of S[1] and S[2] make
the resulting
(Integrals-1) = 0.
(I know that the way I have the parameters set up the equations are very
similar but I am interested in changing the parameters once I have the code
working.)
My attempt at code:
```{r}
a11 <- 1 #alpha_{11}
a12 <- 1 #alpha_{12}
a21 <- 1 #alpha_{21}
a22 <- 1 #alpha_{22}
b1 <- 2 #beta1
b2 <- 2 #beta2
d1 <- 1 #delta1
d2 <- 1 #delta2
g <- 0.5 #gamma
integrand1 <- function(x,S) {b1*g/d1*exp(-g*x)*(1-exp(-d1*
x))*exp(-a11*b1*S[1]/d1*(1-exp(-d1*x))-a12*b2*S[2]/d2*(1-exp(-d2*x)))}
integrand2 <- function(x,S) {b2*g/d2*exp(-g*x)*(1-exp(-d2*
x))*exp(-a22*b2*S[2]/d2*(1-exp(-d2*x))-a21*b1*S[1]/d1*(1-exp(-d1*x)))}
#defining equation we would like to solve
intfun1<- function(S) {integrate(function(x) integrand1(x,
S),lower=0,upper=Inf)[[1]]-1}
intfun2<- function(S) {integrate(function(x) integrand2(x,
S),lower=0,upper=Inf)[[1]]-1}
#putting both equations into one term
model <- function(S) c(F1 = intfun1,F2 = intfun2)
#Solving for roots
(ss <-multiroot(f=model, start=c(0,0)))
```
This gives me the error Error in stode(y, times, func, parms = parms, ...) :
REAL() can only be applied to a 'numeric', not a 'list'
However this simpler example works fine:
```{r}
#Defining the functions
model <- function(x) c(F1 = x[1]+ 4*x[2] -8,F2 = x[1]-4*x[2])
#Solving for the roots
(ss <- multiroot(f = model, start = c(0,0)))
```
Giving me the required x_1= 4 and x_2 =1.
I was given some code to perform a least squares analysis on the same
system but I neither understand the code, nor believe that it is doing what
I am looking for as different initial values give wildly different S values.
```{r}
a11 <- 1 #alpha_{11}
a12 <- 1 #alpha_{12}
a21 <- 1 #alpha_{21}
a22 <- 1 #alpha_{22}
b1 <- 2 #beta1
b2 <- 2 #beta2
d1 <- 1 #delta1
d2 <- 1 #delta2
g <- 0.5 #gamma
integrand1 <- function(x,S) {b1*g/d1*exp(-g*x)*(1-exp(-d1*
x))*exp(-a11*b1*S[1]/d1*(1-exp(-d1*x))-a12*b2*S[2]/d2*(1-exp(-d2*x)))}
integrand2 <- function(x,S) {b2*g/d2*exp(-g*x)*(1-exp(-d2*
x))*exp(-a22*b2*S[2]/d2*(1-exp(-d2*x))-a21*b1*S[1]/d1*(1-exp(-d1*x)))}
#defining equation we would like to solve
intfun1<- function(S) {integrate(function(x)integrand1(x,
S),lower=0,upper=Inf)[[1]]-1}
intfun2<- function(S) {integrate(function(x)integrand2(x,
S),lower=0,upper=Inf)[[1]]-1}
#putting both equations into one term
model <- function(S) if(any(S<0))NA else intfun1(S)**2+ intfun2(S)**2
#Solving for roots
optim(c(0,0), model)
```
I appreciate any tips/help as I have been struggling with this for some
weeks now.
thank you,
--
Ursula
Ph.D. student, University of Michigan
Applied and Interdisciplinary Mathematics
utri...@umich.edu
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
