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.

> # Steady-State solution for reaction rate equations
> # Shacham homotopy method (discrete changing of one or more parameters)
> # M. Shacham: Numerical Solution of Constrained Non-linear algebriac equations
> # International Journal for Numerical Methods in Engineering, 1986, pp.1455--1481.
> 
> # solution should always be > 0
> 
> library(nleqslv)
> 
> RNGkind(kind="Wichmann-Hill")
> set.seed(123)
> 
> # Problem 2, page 1463/1464
> 
> cutlip <- function(x) {
+     # paper has wrong order of parameters
+     # use the Fortran program to get the correct values
+ 
+     # parameter set 2
+     k1 <-  17.721
+     k2 <-  3.483
+     k3 <-  505.051
+     kr1<-  0.118
+     kr2<-  0.033
+ 
+     r <- numeric(6)
+ 
+     r[1] = 1 - x[1] - k1*x[1]*x[6] + kr1 * x[4]
+     r[2] = 1 - x[2] - k2*x[2]*x[6] + kr2 * x[5]
+     r[3] = -x[3] + 2*k3*x[4]*x[5]
+     r[4] = k1*x[1]*x[6] - kr1*x[4] - k3*x[4]*x[5]
+     r[5] = 1.5*(k2*x[2]*x[6] - kr2*x[5]) - k3*x[4]*x[5]
+     r[6] = 1 - x[4] - x[5] - x[6]
+ 
+     r
+ }
> 
> 
> Nrep <- 50
> xstart <- matrix(0,nrow=Nrep, ncol=6)
> xstart[,1] <- runif(Nrep,min=0,max=2)
> xstart[,2] <- runif(Nrep,min=0,max=1)
> xstart[,3] <- runif(Nrep,min=0,max=2)
> xstart[,4] <- runif(Nrep,min=0,max=1)
> xstart[,5] <- runif(Nrep,min=0,max=1)
> xstart[,6] <- runif(Nrep,min=0,max=1)
> 
> ans <- searchZeros(xstart,cutlip, method="Broyden",global="dbldog")
> nrow(ans$x)==4
[1] TRUE
> all(ans$xfnorm <= 1e-10)
[1] TRUE
> 
> zans <- searchZeros(ans$xstart,cutlip, method="Broyden",global="dbldog")
> length(zans$idxcvg)==4
[1] TRUE
> all(ans$xfnorm == zans$xfnorm)
[1] TRUE
> 
> proc.time()
   user  system elapsed 
   0.15    0.07    0.21