[R] solving integral equations with undefined parameters using multiroot

[Ursula Trigos-Raczkowski]

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
utrigos using umich.edu

	[[alternative HTML version deleted]]