[R] solving integral equations with undefined parameters using multiroot

Abbs Spurdle @purd|e@@ @end|ng |rom gm@||@com
Thu May 6 10:39:00 CEST 2021


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.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: solution.png
Type: image/png
Size: 36474 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20210506/561a8bc4/attachment.png>


More information about the R-help mailing list