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]]