R Under development (unstable) (2023-11-25 r85635 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # http://stackoverflow.com/questions/29134996/solving-nonlinear-equation-in-r
> 
> # wants to know if system has closed form solution
> # I want to see how nleqslv behaves
> 
> set.seed(29)
> 
> library(nleqslv)
> 
> print.result <- function(z, do.print.xf=FALSE) {
+     if( do.print.xf ) {
+         print(z$x)
+         print(z$fvec)
+     }
+     print(z$message)
+     print(all(abs(z$fvec)<=1e-8))
+ }
> 
> f <- function(X, a, b, c1, c2, c3) {
+     Y <- numeric(3)
+     x <- X[1]
+     y <- X[2]
+     z <- X[3]
+     Y[1] <- x + y - x*y - c1
+     Y[2] <- x + z - x*z - c2
+     Y[3] <- a*y + b*z - c3
+     return(Y)
+ }
> 
> Jac <- function(X, a, b, c1, c2, c3) {
+     J <- matrix(0,nrow=3,ncol=3)
+     x <- X[1]
+     y <- X[2]
+     z <- X[3]
+ 
+     J[1,1] <- 1-y
+     J[2,1] <- 1-z
+     J[3,1] <- 0
+     J[1,2] <- 1-x
+     J[2,2] <- 0
+     J[3,2] <- a
+     J[1,3] <- 0
+     J[2,3] <- 1-x
+     J[3,3] <- b
+     J
+ }
> 
> a <- 1
> b <- 1
> c1 <- 2
> c2 <- 3
> c3 <- 4
> 
> # exact solution
> x <- (a*c1+b*c2-c3)/(a+b-c3)
> y <- (b*c1-b*c2-c1*c3+c3)/(-a*c1+a-b*c2+b)
> z <- (a*(c1-c2)+(c2-1)*c3)/(a*(c1-1)+b*(c2-1))
> xsol <- c(x,y,z)
> 
> X.start <- c(1,2,3)
> z1 <- nleqslv(X.start,f,Jac,a=a,b=b,c1=c1,c2=c2,c3=c3,
+                 method="Newton",control=list(trace=0,allowSingular=TRUE))
> 
> z2 <- nleqslv(X.start,f,Jac,a=a,b=b,c1=c1,c2=c2,c3=c3,
+         method="Broyden",control=list(trace=0,allowSingular=TRUE))
> 
> all.equal(z1$x,xsol)
[1] TRUE
> all.equal(z2$x,xsol)
[1] TRUE
> print.result(z1)
[1] "Function criterion near zero"
[1] TRUE
> print.result(z2)
[1] "Function criterion near zero"
[1] TRUE
> 
> proc.time()
   user  system elapsed 
   0.12    0.07    0.18