I have two simple equations.
46.85 = r/k
8646.709 = r/(k^2)
I am trying to solve for r and k and I tried setting up my equations as follows
model <- function(r,k) { c(46.85 = r/k,
8646.709 = r/(k^2))}
ss <- multiroot(f = model, start = c(1, 1))
I am seeing some errors. Not sure where I am going wrong. Any advice on how to solve this equation for r and k is much appreciated.
You appear to be using package rootSolve. The way you have specified your function is very wrong. You should have something like this
library(rootSolve)
model <- function(par) {
y <- numeric(2)
r <- par[1]
k <- par[2]
y[1] <- 46.85 - r/k
y[2] <- 8646.709 - r/(k^2)
y
}
and then try this:
xstart <- c(1,1)
multiroot(model, xstart)
multiroot cannot find a solution. Output is
$root
[1] -203307927145 -204397435886
$f.root
[1] 45.85533 8646.70900
$iter
[1] 3
$estim.precis
[1] 4346.282
which is not a solution. Trying other starting values doesn't help.
I will demonstrate two ways of solving your system of equations: manually and using another package. You can solve your equations manually like this: From the first equation we have
r <- 46.85 * k
Substitute this in the second equation and simplfy and we get
k <- 46.85/8646.709
Insert the value found for k in the equation for r
and display the values for r and k
> r
[1] 0.2538448
> k
[1] 0.005418246
and then run the model function like this
model(c(r,k))
giving output
[1] 0 0
The second method involves using another package namely nleqslv.
It has more methods and strategies for finding solutions of system of nonlinear equations.
It also has a testfunction for investigating which method and/or strategy works. Like this
library(nleqslv)
xstart <- c(1,1)
testnslv(xstart,model)
Output is
Call:
testnslv(x = xstart, fn = model)
Results:
Method Global termcd Fcnt Jcnt Iter Message Fnorm
1 Newton cline 1 68 13 13 Fcrit 1.009e-19
2 Newton qline 1 75 17 17 Fcrit 4.896e-19
3 Newton gline 3 90 11 11 Stalled 2.952e+07
4 Newton pwldog 1 91 50 50 Fcrit 2.382e-22
5 Newton dbldog 1 97 54 54 Fcrit 3.243e-22
6 Newton hook 1 72 41 41 Fcrit 6.359e-21
7 Newton none 5 1 2 2 Illcond 3.738e+07
8 Broyden cline 2 88 4 20 Xcrit 8.744e-14
9 Broyden qline 2 153 4 32 Xcrit 1.111e-11
10 Broyden gline 3 156 7 16 Stalled 1.415e+07
11 Broyden pwldog 4 261 13 150 Maxiter 1.462e+03
12 Broyden dbldog 4 288 20 150 Maxiter 1.618e+03
13 Broyden hook 1 192 7 90 Fcrit 1.340e-22
14 Broyden none 5 2 1 3 Illcond 3.738e+07
Read the help page for function nleqslv from start to finish to see what the methods and global columns mean. From the output it seems that the Broyden method is not very successful.
You can also see that a standard Newton method (with none in the global column) cannot find a solution.
So let's try the cline global strategy in a call to nleqslv.
nleqslv(xstart,model,method="Newton",global="cline")
with output
$x
[1] 0.253844844 0.005418246
$fvec
[1] 8.100187e-13 4.492904e-10
$termcd
[1] 1
$message
[1] "Function criterion near zero"
$scalex
[1] 1 1
$nfcnt
[1] 68
$njcnt
[1] 13
$iter
[1] 13
Comparing the solution found with the manually calculated solution shows that nleqslv found the solution.
You can try the rootSolve package with a little restructuring of the equations:
library(rootSolve)
model <- function(x) {
F1 <- x[1] - 46.85*x[2]
F2 <- x[1] - 8646.709 * (x[2]^2)
c(F1 = F1, F2 = F2)
}
ss <- multiroot(f = model, start = c(1, 1))
Considering x[1] as r and x[2] as k, this gives :
print(ss)
#$root
#[1] 0.253844844 0.005418246
#
#$f.root
# F1 F2
#-2.164935e-15 -6.191553e-11
#
#$iter
#[1] 13
#
#$estim.precis
#[1] 3.095885e-11
