Solving simple linear function in R

Question

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.

Answer 1

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.

Answer 2

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