Nothing
# 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
all(ans$xfnorm <= 1e-10)
zans <- searchZeros(ans$xstart,cutlip, method="Broyden",global="dbldog")
length(zans$idxcvg)==4
all(ans$xfnorm == zans$xfnorm)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.