Abbs Spurdle
2021-May-06 09:27 UTC
[R] solving integral equations with undefined parameters using multiroot
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 at 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 at 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 at umich.edu > > > > [[alternative HTML version deleted]] > > > > ______________________________________________ > > R-help at 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
2021-May-06 10:17 UTC
[R] solving integral equations with undefined parameters using multiroot
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 at 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 at 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 at 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 at umich.edu > > > > > > [[alternative HTML version deleted]] > > > > > > ______________________________________________ > > > R-help at 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 at umich.edu [[alternative HTML version deleted]]