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.