# R/ex_pollution.R In deTestSet: Test Set for Differential Equations

#### Documented in pollution

```## =============================================================================
##
## Pollution problem, chemistry
##
## This code is derived from the Test Set for IVP solvers
##     http://www.dm.uniba.it/~testset/
##     ODE of dimension 20
##
## =============================================================================

pollution <- function(times = seq(0, 60, 1), yini = NULL,
parms = list(), printmescd = TRUE,
method = mebdfi, atol = 1e-6, rtol = 1e-6,  ...) {

### check input
parameter <- c(
k1 = .35, k2 = .266e2, k3 = .123e5, k4 = .86e-3, k5 = .82e-3,
k6 = .15e5, k7 = .13e-3, k8 = .24e5,k9 = .165e5,
k10 = .9e4, k11 = .22e-1, k12 = .12e5, k13 = .188e1,
k14 = .163e5, k15 = .48e7, k16 = .35e-3, k17 = .175e-1,
k18 = .1e9, k19 = .444e12, k20 = .124e4, k21 = .21e1,
k22 = .578e1, k23 = .474e-1, k24 = .178e4, k25 = .312e1)

parameter <- overrulepar(parameter, parms, 25)

if (is.null(yini))  {
yini <- rep(0, 20)
yini[2]  <- 0.2
yini[4]  <- 0.04
yini[7]  <- 0.1
yini[8]  <- 0.3
yini[9]  <- 0.01
yini[17] <- 0.007
}
checkini(20, yini)

if (is.null(names(yini)))
names(yini) <- c("NO2", "NO", "O3P", "O3", "HO2", "OH", "HCHO",
"CO", "ALD", "MEO2", "C2O3", "CO2", "PAN", "CH3O",
"HNO3", "O1D", "SO2", "SO4", "NO3", "N2O5")

prob <- polluprob()
### solve
useres <- FALSE
if (is.character(method)) {
if (method %in% c("mebdfi", "daspk"))
useres <- TRUE
} else  if("res" %in% names(formals(method)))
useres <- TRUE
if (useres){

#  out <- ode(y = yini, times = times, func = "polfunc",
# 	dllname = "deTestSet", initfunc = "polpar", method=method,
# 	parms = parameter, ...)
dyini <- rep(0,20)

checkini(20, yini, dyini)
out <-dae(y = yini, dy = dyini, times = times, res = "polres",
dllname = "deTestSet", jacres = "poljacres",initfunc = "polpar",
parms = parameter, method=method,atol=atol, rtol=rtol,  ...)
}
else
out <- ode(y = yini, times = times, func = "polfunc", jacfunc = "poljac",
dllname = "deTestSet", initfunc = "polpar", method=method,
parms = parameter,atol=atol, rtol=rtol, ...)
if(printmescd)
out <- printpr (out, prob, "pollution", rtol, atol)
return(out)
}

polluprob <- function(){
fullnm <- 'Pollution problem'
problm <- 'pollu'
type   <- 'ODE'
neqn   <- 20
t <- matrix(1,2)
t[1]   <- 0
t[2]   <- 60
numjac <- FALSE
mljac  <- neqn
mujac  <- neqn
return(list(fullnm=fullnm, problm=problm,type=type,neqn=neqn,
t=t,numjac=numjac,mljac=mljac,mujac=mujac))
}
```

## Try the deTestSet package in your browser

Any scripts or data that you put into this service are public.

deTestSet documentation built on May 29, 2017, 2:56 p.m.