create_calc_DO: Create a function to compute the numerical integration of a...

In USGS-R/streamMetabolizer: Models for Estimating Aquatic Photosynthesis and Respiration

Create a function to compute the numerical integration of a dDOdt function

Description

Create a function to compute the numerical integration of a dDOdt function

Usage

create_calc_DO(
  calc_dDOdt,
  ode_method = environment(calc_dDOdt)$ode_method,
  err.obs = 0
)

Arguments

calc_dDOdt	a function as from create_calc_dDOdt
ode_method	character. The method to use in solving the ordinary differential equation for DO. Options: euler, formerly Euler: the final change in DO from t=1 to t=2 is solely a function of GPP, ER, DO, etc. at t=1 trapezoid, formerly pairmeans: the final change in DO from t=1 to t=2 is a function of the mean values of GPP, ER, etc. across t=1 and t=2. for type='mle', options also include rk2 and any character method accepted by ode in the deSolve package (lsoda, lsode, lsodes, lsodar, vode, daspk, rk4, ode23, ode45, radau, bdf, bdf_d, adams, impAdams, and impAdams_d; note that many of these have not been well tested in the context of streamMetabolizer models)
err.obs	optional numerical vector of length nrow(data) in units of gO2 m^3. Appropriate for simulation, when this vector of observation errors will be added to the calculated DO values to simulate observation error. But usually (for MLE or prediction from a fitted MLE/Bayesian/nighttime regression model) err.obs should be missing or 0

Value

a function that will return a negative log likelihood of the data given a set of metab.pars

Examples

## Not run: 
# prepare data for examples
data <- data_metab('3','30')[97:144,][seq(1,48,by=2),]
# preds.init <- list(GPP.daily=2.82,ER.daily=-2.12,K600.daily=31.27)
preds.init <- as.list(dplyr::select(
  predict_metab(metab(specs(mm_name('mle', ode_method='trapezoid')), data=data)),
  GPP.daily=GPP, ER.daily=ER, K600.daily=K600))
DOtime <- data$solar.time
dDOtime <- data$solar.time[-nrow(data)] + (data$solar.time[2] - data$solar.time[1])/2

# integration of dDOdt by euler, trapezoid, rk2, rk4, and lsoda methods
plot(x=DOtime, y=data$DO.obs, pch=3, cex=0.6)
# euler
dDOdt <- create_calc_dDOdt(data, ode_method='euler', GPP_fun='linlight',
  ER_fun='constant', deficit_src='DO_mod')
DO <- create_calc_DO(dDOdt, ode_method='euler')
DO.mod.euler <- DO(metab.pars=preds.init)
lines(x=DOtime, y=DO.mod.euler, type='l', col='chartreuse3')
# trapezoid=pairmeans
dDOdt <- create_calc_dDOdt(data, ode_method='trapezoid', GPP_fun='linlight',
  ER_fun='constant', deficit_src='DO_mod')
DO <- create_calc_DO(dDOdt, ode_method='trapezoid')
DO.mod.trap <- DO(metab.pars=preds.init)
lines(x=DOtime, y=DO.mod.trap, type='l', col='gold')
# lsoda
dDOdt <- create_calc_dDOdt(data, ode_method='lsoda', GPP_fun='linlight',
  ER_fun='constant', deficit_src='DO_mod')
DO <- create_calc_DO(dDOdt, ode_method='lsoda')
DO.mod <- DO(metab.pars=preds.init)
lines(x=DOtime, y=DO.mod, type='l', col='navy')
# rk2
dDOdt <- create_calc_dDOdt(data, ode_method='rk2', GPP_fun='linlight',
  ER_fun='constant', deficit_src='DO_mod')
DO <- create_calc_DO(dDOdt, ode_method='rk2')
DO.mod <- DO(metab.pars=preds.init)
lines(x=DOtime, y=DO.mod, type='l', col='blue')
# rk4
dDOdt <- create_calc_dDOdt(data, ode_method='rk4', GPP_fun='linlight',
  ER_fun='constant', deficit_src='DO_mod')
DO <- create_calc_DO(dDOdt, ode_method='rk4')
DO.mod <- DO(metab.pars=preds.init)
lines(x=DOtime, y=DO.mod, type='l', col='magenta')

# with observation and/or process error
plot(x=DOtime, y=data$DO.obs, col='black', pch=3, cex=0.6)
dDOdt <- create_calc_dDOdt(data, ode_method='trapezoid', GPP_fun='linlight',
  ER_fun='constant', deficit_src='DO_mod')
DO <- create_calc_DO(dDOdt, ode_method='trapezoid', err.obs=rnorm(nrow(data), 0, 0))
DO.mod.obs <- DO(preds.init)
lines(x=DOtime, y=DO.mod.obs, col='black', lty=1)
dDOdt <- create_calc_dDOdt(data, ode_method='trapezoid', GPP_fun='linlight',
  ER_fun='constant', deficit_src='DO_mod')
DO <- create_calc_DO(dDOdt, ode_method='trapezoid', err.obs=rnorm(nrow(data), 0, 0.3))
DO.mod.obserr <- DO(preds.init)
lines(x=DOtime, y=DO.mod.obserr, col='blue', lty=2)
dDOdt <- create_calc_dDOdt(data, ode_method='trapezoid', GPP_fun='linlight',
  ER_fun='constant', deficit_src='DO_mod', err.proc=rnorm(nrow(data), 0, 2))
DO <- create_calc_DO(dDOdt, ode_method='trapezoid', err.obs=rnorm(nrow(data), 0, 0.1))
DO.mod.operr <- DO(preds.init)
lines(x=DOtime, y=DO.mod.operr, col='red', lty=2)

## End(Not run)

USGS-R/streamMetabolizer documentation built on Aug. 15, 2023, 7:50 a.m.