The problem is a stiff system of 66 ordinary differential equations.
The 'Mathematics and the Environment'
project group at CWI contributed this problem to the test set.
The software part of the
problem is in the file emep.f available at [MM08].
A matrix of class deSolve with up to as many rows as elements in
times and as many
columns as elements in yini, plus an additional column (the first)
for the time value.
There will be one row for each element in times unless the
solver returns with an unrecoverable error. If
yini has a names attribute, it will be used to label the columns
of the output value.
[MM08] F. Mazzia and C. Magherini. Test Set for Initial Value Problem Solvers,
release 2.4. Department
of Mathematics, University of Bari and INdAM, Research Unit of Bari, February
2008.
[SASJ93] D. Simpson, Y. Andersson-Skold, and M.E. Jenkin. Updating the chemical scheme for the
EMEP MSC-W model: Current status. Report EMEP MSC-W Note 2/93, The Norwegian
Meteorological Institute, Oslo, 1993.
[Sim93] D. Simpson. Photochemical model calculations over Europe for two extended summer
periods: 1985 and 1989. model results and comparisons with observations. Atmospheric
Environment, 27A:921-943, 1993.
[VS94] J.G. Verwer and D. Simpson. Explicit methods for stiff ODEs from atmospheric chemistry.
Report NM-R9409, CWI, Amsterdam, 1994.
Examples
out <- emep()
plot(out, lwd = 2, col = "darkblue",
which = c("NO", "NO2", "SO2", "CH4", "O3", "N2O5"))
plot(out, col = "darkblue", lwd = 2, which = 1:16)
mtext(side = 3,line = -1.5, "emep", cex = 1.25, outer = TRUE)
# compare with reference solution (component 36 and 38 not included)
refsol <- reference("emep")
inderr <- c(1:35,37,39:66)
max(abs(out[nrow(out),inderr+1] - refsol[inderr])/refsol[inderr])