Solving a system of nonlinear equations in R
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
solving a simple (?) system of nonlinear equations
Your function does not reflect properly what you want.
You can see this by evaluating fn(c(0.3,0.1))
as follows.
fn(c(0.3,0.1))
[1] 0.3100255 0.1192029
So the output is very close to the input. You wanted (almost) zero as output.
So you want to solve the system for p
and q
.
What you need to do is to make your function return the difference between the input p
and the expression for pstar
and the difference between the input q
and the expression for qstar
.
So rewrite your function as follows
fn <- function(x, lambda = 1){
p <- x[1]
q <- x[2]
pstar <- exp(lambda * (1*x[2])) / (exp(lambda * (1*x[2])) + exp(lambda * (1 - x[2])))
qstar <- exp(lambda * (1 - x[1])) / (exp(lambda * ((1 - x[1]))) + exp(lambda * (9*x[1])))
return(c(pstar-p,qstar-q))
}
and then call nleqslv
as follows (PLEASE always show all the code you are using. You left out the library(nleqslv)
).
library(nleqslv)
xstart <- c(0.1, 0.3)
nleqslv(xstart, fn)
This will display the full output of the function. Always a good idea to check for succes. Always check $termcd
for succes.
$x
[1] 0.3127804 0.1064237
$fvec
[1] 5.070055e-11 6.547240e-09
$termcd
[1] 1
$message
[1] "Function criterion near zero"
$scalex
[1] 1 1
$nfcnt
[1] 7
$njcnt
[1] 1
$iter
[1] 7
The result for $x
is more what you expect.
Finally please use <-
for assignment. If you don't there will come the day that you will be bitten by R
and its magic.
This is nothing wrong in using nleqslv
for this problem. You only made a small mistake.
Solve system of non-linear equations
Your system of equations has multiple solutions.
I use a different package to solve your system: nleqslv
as follows:
library(nleqslv)
model <- function(x) {
F1 <- sqrt(x[1]^2 + x[3]^2) - 1
F2 <- sqrt(x[2]^2 + x[4]^2) - 1
F3 <- x[1]*x[2] + x[3]*x[4]
F4 <- -0.58*x[2] - 0.19*x[3]
c(F1 = F1, F2 = F2, F3 = F3, F4 = F4)
}
#find solution
xstart <- c(1.5, 0, 0.5, 0)
nleqslv(xstart,model)
This gets the same solution as the answer of Prem.
Your system however has multiple solutions.
Package nleqslv
provides a function to search for solutions given a matrix of different starting values. You can use this
set.seed(13)
xstart <- matrix(runif(400,0,2),ncol=4)
searchZeros(xstart,model)
(Note: different seeds may not find all four solutions)
You will see that there are four different solutions:
$x
[,1] [,2] [,3] [,4]
[1,] -1 -1.869055e-10 5.705536e-10 -1
[2,] -1 4.992198e-13 -1.523934e-12 1
[3,] 1 -1.691309e-10 5.162942e-10 -1
[4,] 1 1.791944e-09 -5.470144e-09 1
.......
This clearly suggests that the exact solutions are as given in the following matrix
xsol <- matrix(c(1,0,0,1,
1,0,0,-1,
-1,0,0,1,
-1,0,0,-1),byrow=TRUE,ncol=4)
And then do
model(xsol[1,])
model(xsol[2,])
model(xsol[3,])
model(xsol[4,])
Confirmed!
I have not tried to find these solutions analytically but you can see that if x[2]
and x[3]
are zero then F3
and F4
are zero. The solutions for x[1]
and x[4]
can then be immediately found.
Related Topics
R - Scaling Numeric Values Only in a Dataframe with Mixed Types
Loop for Reverse Geocoding in R
R: Faceted Bar Chart with Percentages Labels Independent for Each Plot
Difference of Prediction Results in Random Forest Model
Rgdal Installation Difficulty on Ubuntu 16.04 Lts
Error While Using Install_Github | Devtools | Timeout Issue
Shiny: Unwanted Space Added by Plotoutput() And/Or Renderplot()
Directly Adding Titles and Labels to Visnetwork
Using Pivot_Longer with Multiple Paired Columns in the Wide Dataset
Fill in Data Frame with Values from Rows Above
How to Display Strip Labels Below the Plot When Faceting
Flatten Nested List into 1-Deep List
Categorical Scatter Plot with Mean Segments Using Ggplot2 in R
Function Composition in R (And High Level Functions)
Axis Labels for Each Bar and Each Group in Bar Charts with Dodged Groups