Steadystate solver for multicomponent 1D ordinary differential equations
Description
Estimates the steadystate condition for a system of ordinary differential equations that result from 1Dimensional partial differential equation models that have been converted to ODEs by numerical differencing.
It is assumed that exchange occurs only between adjacent layers.
Usage
1 2 3 4 
Arguments
y 
the initial guess of (state) values for the ODE system, a vector. 
time, times 
time for which steadystate is wanted;
the default is 
func 
either an Rfunction that computes the values of the
derivatives in the ode system (the model defininition) at time If The return value of The derivatives
should be specified in the same order as the state variables 
parms 
parameters passed to 
nspec 
the number of *species* (components) in the model.
If 
dimens 
the number of *boxes* in the model. If NULL, then

names 
the names of the components; used to label the output, which will be written as a matrix. 
method 
the solution method, one of "stode", "stodes", or "runsteady".
When 
jactype 
the jacobian type  default is a regular 1D structure where layers only interact with adjacent layers.
If the structure does not comply with this, the jacobian can be set equal to 
cyclicBnd 
if a cyclic boundary exists, a value of 
bandwidth 
the number of adjacent boxes over which transport occurs.
Normally equal to 1 (box i only interacts with box i1, and i+1).
Values larger than 1 will not work with 
... 
additional arguments passed to the solver function as defined
by 
Details
This is the method of choice for multispecies 1dimensional models, that are only subjected to transport between adjacent layers
More specifically, this method is to be used if the state variables are arranged per species:
A[1],A[2],A[3],....B[1],B[2],B[3],.... (for species A, B))
Two methods are implemented.
The default method rearranges the state variables as A[1],B[1],...A[2],B[2],...A[3],B[3],.... This reformulation leads to a banded Jacobian with (upper and lower) half bandwidth = number of species. Then function
stode
solves the banded problem.The second method uses function
stodes
. Based on the dimension of the problem, the method first calculates the sparsity pattern of the Jacobian, under the assumption that transport is only occurring between adjacent layers. Thenstodes
is called to solve the problem.As
stodes
is used to estimate steadystate, it may be necessary to specify the length of the real work array,lrw
.Although a reasonable guess of
lrw
is made, it is possible that this will be too low. In this case,steady.1D
will return with an error message telling the size of the work array actually needed. In the second try then, setlrw
equal to this number.
For singlespecies 1D models, use steady.band
.
If state variables are arranged as (e.g. A[1],B[1],A[2],B[2],A[3],B[3],...
then the model should be solved with steady.band
In some cases, a cyclic boundary condition exists. This is when the first
box interacts with the last box and vice versa. In this case, there
will be extra nonzero fringes in the Jacobian which need to be taken
into account. The occurrence of cyclic boundaries can be
toggled on by specifying argument cyclicBnd=1
. If this is the
case, then the steadystate will be estimated using stodes
.
The default is no cyclic boundaries.
Value
A list containing
y 
if 
... 
the number of "global" values returned. 
The output will have the attribute steady
, which returns TRUE, if
steadystate has been reached and the attribute
precis
with the precision attained during each iteration.
Note
It is advisable though not mandatory to specify BOTH nspec
and
dimens
.
In this case, the solver can check whether the input makes sense
(i.e. if nspec*dimens = length(y))
Author(s)
Karline Soetaert <karline.soetaert@nioz.nl>
See Also
plot.steady1D
for plotting the output of steady.1D
steady
, for a general interface to most of the steadystate
solvers
steady.band
, to find the steadystate of ODE models with a
banded Jacobian
steady.2D
,
steady.3D
, steadystate solvers for 2D and 3D
partial differential equations.
stode
, iterative steadystate solver for ODEs with full
or banded Jacobian.
stodes
, iterative steadystate solver for ODEs with arbitrary
sparse Jacobian.
runsteady
, steadystate solver by dynamically running to
steadystate
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200  ## =======================================================================
## EXAMPLE 1: BOD + O2
## =======================================================================
## Biochemical Oxygen Demand (BOD) and oxygen (O2) dynamics
## in a river
#==================#
# Model equations #
#==================#
O2BOD < function(t, state, pars) {
BOD < state[1:N]
O2 < state[(N+1):(2*N)]
# BOD dynamics
FluxBOD < v * c(BOD_0, BOD) # fluxes due to water transport
FluxO2 < v * c(O2_0, O2)
BODrate < r*BOD*O2/(O2+10) # 1st order consumption, Monod in oxygen
#rate of change = flux gradient  consumption + reaeration (O2)
dBOD < diff(FluxBOD)/dx  BODrate
dO2 < diff(FluxO2)/dx  BODrate + p*(O2satO2)
return(list(c(dBOD = dBOD, dO2 = dO2), BODrate = BODrate))
} # END O2BOD
#==================#
# Model application#
#==================#
# parameters
dx < 100 # grid size, meters
v < 1e2 # velocity, m/day
x < seq(dx/2, 10000, by = dx) # m, distance from river
N < length(x)
r < 0.1 # /day, firstorder decay of BOD
p < 0.1 # /day, airsea exchange rate
O2sat < 300 # mmol/m3 saturated oxygen conc
O2_0 < 50 # mmol/m3 riverine oxygen conc
BOD_0 < 1500 # mmol/m3 riverine BOD concentration
# initial guess:
state < c(rep(200, N), rep(200, N))
# running the model
print(system.time(
out < steady.1D (y = state, func = O2BOD, parms = NULL,
nspec = 2, pos = TRUE, names = c("BOD", "O2"))))
#==================#
# Plotting output #
#==================#
mf < par(mfrow = c(2, 2))
plot(x, out$y[ ,"O2"], xlab = "Distance from river",
ylab = "mmol/m3", main = "Oxygen", type = "l")
plot(x, out$y[ ,"BOD"], xlab = "Distance from river",
ylab = "mmol/m3", main = "BOD", type = "l")
plot(x, out$BODrate, xlab = "Distance from river",
ylab = "mmol/m3/d", main = "BOD decay rate", type = "l")
par(mfrow=mf)
# same plot in one command
plot(out, which = c("O2","BOD","BODrate"),xlab = "Distance from river",
ylab = c("mmol/m3","mmol/m3","mmol/m3/d"),
main = c("Oxygen","BOD","BOD decay rate"), type = "l")
# same, but now running dynamically to steadystate
print(system.time(
out2 < steady.1D (y = state, func = O2BOD, parms = NULL, nspec = 2,
time = c(0, 1000), method = "runsteady",
names = c("BOD", "O2"))))
# plot all state variables at once, against "x":
plot(out2, grid=x, xlab = "Distance from river",
ylab = "mmol/m3", type = "l", lwd = 2)
plot(out, out2, grid=x, xlab = "Distance from river", which = "BODrate",
ylab = "mmol/m3", type = "l", lwd = 2)
## =======================================================================
## EXAMPLE 2: Silicate diagenesis
## =======================================================================
## Example from the book:
## Soetaert and Herman (2009).
## a practical guide to ecological modelling 
## using R as a simulation platform.
## Springer
#====================#
# Model equations #
#====================#
SiDIAmodel < function (time = 0, # time, not used here
Conc, # concentrations: BSi, DSi
parms = NULL) # parameter values; not used
{
BSi< Conc[1:N]
DSi< Conc[(N+1):(2*N)]
# transport
# diffusive fluxes at upper interface of each layer
# upper concentration imposed (bwDSi), lower: zero gradient
DSiFlux < SedDisp * IntPor *diff(c(bwDSi ,DSi,DSi[N]))/thick
BSiFlux < Db *(1IntPor)*diff(c(BSi[1],BSi,BSi[N]))/thick
BSiFlux[1] < BSidepo # upper boundary flux is imposed
# BSi dissolution
Dissolution < rDissSi * BSi*(1. DSi/EquilSi )^pow
Dissolution < pmax(0,Dissolution)
# Rate of change= Flux gradient, corrected for porosity + dissolution
dDSi < diff(DSiFlux)/thick/Porosity + # transport
Dissolution * (1Porosity)/Porosity # biogeochemistry
dBSi < diff(BSiFlux)/thick/(1Porosity)  Dissolution
return(list(c(dBSi, dDSi), # Rates of changes
Dissolution = Dissolution, # Profile of dissolution rates
DSiSurfFlux = DSiFlux[1], # DSi sedimentwater exchange rate
DSIDeepFlux = DSiFlux[N+1], # DSi deepwater (burial) flux
BSiDeepFlux = BSiFlux[N+1])) # BSi deepwater (burial) flux
}
#====================#
# Model run #
#====================#
# sediment parameters
thick < 0.1 # thickness of sediment layers (cm)
Intdepth < seq(0, 10, by = thick) # depth at upper interface of layers
Nint < length(Intdepth) # number of interfaces
Depth < 0.5*(Intdepth[Nint] +Intdepth[1]) # depth at middle of layers
N < length(Depth) # number of layers
por0 < 0.9 # surface porosity ()
pordeep < 0.7 # deep porosity ()
porcoef < 2 # porosity decay coefficient (/cm)
# porosity profile, middle of layers
Porosity < pordeep + (por0pordeep)*exp(Depth*porcoef)
# porosity profile, upper interface
IntPor < pordeep + (por0pordeep)*exp(Intdepth*porcoef)
dB0 < 1/365 # cm2/day  bioturbation coefficient
dBcoeff < 2
mixdepth < 5 # cm
Db < pmin(dB0, dB0*exp((Intdepthmixdepth)*dBcoeff))
# biogeochemical parameters
SedDisp < 0.4 # diffusion coefficient, cm2/d
rDissSi < 0.005 # dissolution rate, /day
EquilSi < 800 # equilibrium concentration
pow < 1
BSidepo < 0.2*100 # nmol/cm2/day
bwDSi < 150 # mmol/m3
# initial guess of state variablesjust random numbers between 0,1
Conc < runif(2*N)
# three runs with different deposition rates
BSidepo < 0.2*100 # nmol/cm2/day
sol < steady.1D (Conc, func = SiDIAmodel, parms = NULL, nspec = 2,
names = c("DSi", "BSi"))
BSidepo < 2*100 # nmol/cm2/day
sol2 < steady.1D (Conc, func = SiDIAmodel, parms = NULL, nspec = 2,
names = c("DSi", "BSi"))
BSidepo < 3*100 # nmol/cm2/day
sol3 < steady.1D (Conc, func = SiDIAmodel, parms = NULL, nspec = 2,
names = c("DSi", "BSi"))
#====================#
# plotting output #
#====================#
par(mfrow=c(2,2))
# Plot 3 runs
plot(sol, sol2, sol3, xyswap = TRUE, mfrow = c(2, 2),
xlab = c("mmolSi/m3 liquid", "mmolSi/m3 solid"),
ylab = "Depth", lwd = 2, lty = 1)
legend("bottom", c("0.2", "2", "3"), title = "mmol/m2/d",
lwd = 2, col = 1:3)
plot(Porosity, Depth, ylim = c(10, 0), xlab = "" ,
main = "Porosity", type = "l", lwd = 2)
plot(Db, Intdepth, ylim = c(10, 0), xlab = "cm2/d",
main = "Bioturbation", type = "l", lwd = 2)
mtext(outer = TRUE, side = 3, line = 2, cex = 1.5, "SiDIAmodel")
# similar, but shorter
plot(sol, sol2, sol3, vertical =TRUE,
lwd = 2, lty = 1,
main = c("DSi [mmol/m3 liq]","BSi [mmol/m3 sol]"),
ylab= "depth [cm]")
legend("bottom", c("0.2", "2", "3"), title = "mmol/m2/d",
lwd = 2, col = 1:3)
