[R] solving integral equations with undefined parameters using multiroot

Fri May 7 03:56:40 CEST 2021

#using vF1 function
#from my previous posts
u <- seq (0, 0.25,, 200)
cl <- contourLines (u, u, outer (u, u, vF1),, 0)[[1]]
plot (cl$x, cl$y, type="l")


On Thu, May 6, 2021 at 10:18 PM Ursula Trigos-Raczkowski
<utrigos using umich.edu> wrote:
>
> Thanks for your reply. Unfortunately the code doesn't work even when I change the parameters to ensure I have "different" equations.
> Using mathematica I do see that my two equations form planes, intersecting in a line of infinite solutions but it is not very accurate, I was hoping R would be more accurate and tell me what this line is, or at least a set of solutions.
>
> On Thu, May 6, 2021 at 5:28 AM Abbs Spurdle <spurdle.a using gmail.com> wrote:
>>
>> Just realized five minutes after posting that I misinterpreted your
>> question, slightly.
>> However, after comparing the solution sets for *both* equations, I
>> can't see any obvious difference between the two.
>> If there is any difference, presumably that difference is extremely small.
>>
>>
>> On Thu, May 6, 2021 at 8:39 PM Abbs Spurdle <spurdle.a using gmail.com> wrote:
>> >
>> > Hi Ursula,
>> >
>> > If I'm not mistaken, there's an infinite number of solutions, which
>> > form a straight (or near straight) line.
>> > Refer to the following code, and attached plot.
>> >
>> > ----begin code---
>> > library (barsurf)
>> > vF1 <- function (u, v)
>> > {   n <- length (u)
>> >     k <- numeric (n)
>> >     for (i in seq_len (n) )
>> >         k [i] <- intfun1 (c (u [i], v [i]) )
>> >     k
>> > }
>> > plotf_cfield (vF1, c (0, 0.2), fb = (-2:2) / 10,
>> >     main="(integral_1 - 1)",
>> >     xlab="S[1]", ylab="S[2]",
>> >     n=40, raster=TRUE, theme="heat", contour.labels=TRUE)
>> > ----end code----
>> >
>> > I'm not familiar with the RootSolve package.
>> > Nor am I quite sure what you're trying to compute, given the apparent
>> > infinite set of solutions.
>> >
>> > So, for now at least, I'll leave comments on the root finding to someone who is.
>> >
>> >
>> > Abby
>> >
>> >
>> > On Thu, May 6, 2021 at 8:46 AM Ursula Trigos-Raczkowski
>> > <utrigos using umich.edu> wrote:
>> > >
>> > > 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]]
>> > >
>> > > ______________________________________________
>> > > R-help using 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.
>
>
>
> --
> Ursula Trigos-Raczkowski (she/her/hers)
> Ph.D. student, University of Michigan
> Applied and Interdisciplinary Mathematics
> 5828 East Hall
> 530 Church St.
> Ann Arbor, MI 48109-1085
> utrigos using umich.edu
>