# steady.1D: Steady-state solver for multicomponent 1-D ordinary... In rootSolve: Nonlinear Root Finding, Equilibrium and Steady-State Analysis of Ordinary Differential Equations

## Description

Estimates the steady-state condition for a system of ordinary differential equations that result from 1-Dimensional 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``` ```steady.1D(y, time = NULL, func, parms = NULL, nspec = NULL, dimens = NULL, names = NULL, method = "stode", jactype = NULL, cyclicBnd = NULL, bandwidth = 1, times = time, ...) ```

## Arguments

 `y ` the initial guess of (state) values for the ODE system, a vector. `time, times ` time for which steady-state is wanted; the default is `times=0` (for `method = "stode"` or `method = "stodes"`, and `times = c(0,Inf)` for `method = "runsteady"` ). (note- since version 1.7, 'times' has been added as an alias to 'time'). `func ` either an R-function that computes the values of the derivatives in the ode system (the model defininition) at time `time`, or a character string giving the name of a compiled function in a dynamically loaded shared library. If `func` is an R-function, it must be defined as: `yprime = func(t, y, parms,...)`. `t` is the current time point in the integration, `y` is the current estimate of the variables in the ODE system. If the initial values `y` has a names attribute, the names will be available inside `func`. `parms` is a vector or list of parameters; ... (optional) are any other arguments passed to the function. The return value of `func` should be a list, whose first element is a vector containing the derivatives of `y` with respect to `time`, and whose next elements are global values whose steady-state value is also required. The derivatives should be specified in the same order as the state variables `y`. `parms ` parameters passed to `func`. `nspec ` the number of *species* (components) in the model. If `NULL`, then `dimens` should be specified. `dimens` the number of *boxes* in the model. If NULL, then `nspec` should be specified. `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 `method` = 'runsteady', then solver `lsode` is used by default, unless argument `jactype` is set to `"sparse"`, in which case `lsodes` is used and the structure of the jacobian is determined by the solver. `jactype ` the jacobian type - default is a regular 1-D structure where layers only interact with adjacent layers. If the structure does not comply with this, the jacobian can be set equal to `'sparse'` if `method = stodes`. If `jactype = 'sparse'` and `method = runsteady` then `lsodes`, the sparse integrator, will be used. `cyclicBnd ` if a cyclic boundary exists, a value of `1` else `NULL`; see details. `bandwidth ` the number of adjacent boxes over which transport occurs. Normally equal to 1 (box i only interacts with box i-1, and i+1). Values larger than 1 will not work with `method = "stodes"`. `... ` additional arguments passed to the solver function as defined by `method`.

## Details

This is the method of choice for multi-species 1-dimensional 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. Then `stodes` is called to solve the problem.

As `stodes` is used to estimate steady-state, 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, set `lrw` equal to this number.

For single-species 1-D 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 non-zero 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 steady-state will be estimated using `stodes`. The default is no cyclic boundaries.

## Value

A list containing

 `y ` if `names` is not given: a vector with the state variable values from the last iteration during estimation of steady-state condition of the system of equations. if `names` is given, a matrix with one column for every steady-state *component*. `... ` the number of "global" values returned.

The output will have the attribute `steady`, which returns TRUE, if steady-state 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 <[email protected]>

`plot.steady1D` for plotting the output of steady.1D

`steady`, for a general interface to most of the steady-state solvers

`steady.band`, to find the steady-state of ODE models with a banded Jacobian

`steady.2D`, `steady.3D`, steady-state solvers for 2-D and 3-D partial differential equations.

`stode`, iterative steady-state solver for ODEs with full or banded Jacobian.

`stodes`, iterative steady-state solver for ODEs with arbitrary sparse Jacobian.

`runsteady`, steady-state solver by dynamically running to steady-state

## 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) # 1-st order consumption, Monod in oxygen #rate of change = flux gradient - consumption + reaeration (O2) dBOD <- -diff(FluxBOD)/dx - BODrate dO2 <- -diff(FluxO2)/dx - BODrate + p*(O2sat-O2) 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, first-order decay of BOD p <- 0.1 # /day, air-sea 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 steady-state 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 *(1-IntPor)*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 * (1-Porosity)/Porosity # biogeochemistry dBSi <- -diff(BSiFlux)/thick/(1-Porosity) - Dissolution return(list(c(dBSi, dDSi), # Rates of changes Dissolution = Dissolution, # Profile of dissolution rates DSiSurfFlux = DSiFlux[1], # DSi sediment-water exchange rate DSIDeepFlux = DSiFlux[N+1], # DSi deep-water (burial) flux BSiDeepFlux = BSiFlux[N+1])) # BSi deep-water (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 + (por0-pordeep)*exp(-Depth*porcoef) # porosity profile, upper interface IntPor <- pordeep + (por0-pordeep)*exp(-Intdepth*porcoef) dB0 <- 1/365 # cm2/day - bioturbation coefficient dBcoeff <- 2 mixdepth <- 5 # cm Db <- pmin(dB0, dB0*exp(-(Intdepth-mixdepth)*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 variables-just 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) ```

rootSolve documentation built on May 29, 2017, 3:28 p.m.