6

Suppose I have the following system of equations:

a * b = 5
sqrt(a * b^2) = 10

How can I solve these equations for a and b in R ?

I guess this problem can be stated as an optimisation problem, with the following function... ?

fn <- function(a, b) {

    rate <- a * b
    shape <- sqrt(a * b^2)

    return(c(rate, shape) )

}
Ladislas Nalborczyk
  • 725
  • 2
  • 5
  • 20
  • @G5W Do you have a way to turn this into a system of linear equations? – Dason Feb 16 '18 at 18:33
  • looks like pretty simple math beginning with substituting b = 5/a into the second equation. – G5W Feb 16 '18 at 18:33
  • When I googled your question title, the first hit was documentation for the `nleqslv` package. So I'd start with that. – ngm Feb 16 '18 at 18:33
  • @G5W Well I am asking how to solve this kind of equation in R in a generic way (like using `optim` or `solve`). Thank you @ngm for the `nleqslv` suggestion. I have glanced at the doc but I'm not sure how to use it in order to solve my problem... – Ladislas Nalborczyk Feb 16 '18 at 18:44
  • R doesn't do abstract calculus and doesn't "solve for" anything - you probably want to do it in Mathematica or the like. – gented Feb 16 '18 at 18:45

2 Answers2

20

In a comment the poster specifically asks about using solve and optim so we show how to solve this (1) by hand, (2) using solve, (3) using optim and (4) a fixed point iteration.

1) by hand First note that if we write a = 5/b based on the first equation and substitute that into the second equation we get sqrt(5/b * b^2) = sqrt(5 * b) = 10 so b = 20 and a = 0.25.

2) solve Regarding the use of solve these equations can be transformed into linear form by taking the log of both sides giving:

log(a) + log(b) = log(5)
0.5 * (loga + 2 * log(b)) = log(10)

which can be expressed as:

m <- matrix(c(1, .5, 1, 1), 2)
exp(solve(m, log(c(5, 10))))
## [1]  0.25 20.00

3) optim Using optim we can write this where fn is from the question. fn2 is formed by subtracting off the RHS of the equations and using crossprod to form the sum of squares.

fn2 <- function(x) crossprod( fn(x[1], x[2]) - c(5, 10))
optim(c(1, 1), fn2)

giving:

$par
[1]  0.2500805 19.9958117

$value
[1] 5.51508e-07

$counts
function gradient 
      97       NA 

$convergence
[1] 0

$message
NULL

4) fixed point For this one rewrite the equations in a fixed point form, i.e. in the form c(a, b) = f(c(a, b)) and then iterate. In general, there will be several ways to do this and not all of them will converge but in this case this seems to work. We use starting values of 1 for both a and b and divide both side of the first equation by b to get the first equation in fixed point form and we divide both sides of the second equation by sqrt(a) to get the second equation in fixed point form:

a <- b <- 1  # starting values
for(i in 1:100) {
  a = 5 / b
  b = 10 / sqrt(a)
}

data.frame(a, b)
##      a  b
## 1 0.25 20
G. Grothendieck
  • 254,981
  • 17
  • 203
  • 341
13

Use this library.

library("nleqslv")

You need to define the multivariate function you want to solve for.

fn <- function(x) {

  rate <- x[1] * x[2] - 5
  shape <- sqrt(x[1] * x[2]^2) - 10

  return(c(rate, shape))

}

Then you're good to go.

nleqslv(c(1,5), fn)

Always look at the detailed results. Numerical calculations can be tricky. In this case I got this:

Warning message:
In sqrt(x[1] * x[2]^2) : NaNs produced

That just means the procedure searched a region that included x[1] < 0 and then presumably noped the heck back to the right hand side of the plane.

ngm
  • 2,539
  • 8
  • 18
  • Firstly, you need to install the package by the following code: install.packages("nleqslv"). Then you could use the library. – M Shafaei N Oct 26 '21 at 05:53