steady.1D(rootSolve)
steady.1D()所属R语言包：rootSolve
Steady-state solver for multicomponent 1-D ordinary differential equations
稳态求解器，多组分1-D常微分方程
译者：生物统计家园网 机器人LoveR
描述----------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----------
steady.1D(y, time = NULL, func, parms = NULL,
nspec = NULL, dimens = NULL,
names = NULL, method = "stode",
cyclicBnd = NULL, bandwidth = 1, ...)
参数----------Arguments----------
参数：y
the initial guess of (state) values for the ODE system, a vector.
（州）的初始估计值的ODE系统的一个向量。
参数：time
time for which steady-state is wanted; the default is time=0 (for method = "stode" or method = "stodes", and time = c(0,Inf) for method = "runsteady".
时间的稳定状态被通缉，默认的是time=0（method = "stode"或method = "stodes"和time = c(0,Inf)method = "runsteady"。
参数：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.
可以是R-函数的衍生物中的值，计算在时间的颂歌系统（模型defininition）time，或给人的编译函数在一个动态加载的共享库的名称的字符串。如果func是一个R-函数，它必须被定义为：yprime = func(t, y, parms,...)。 t是当前时间点的整合，y是目前估计的ODE系统中的变量。如果初始值y有一个名称属性，名称将是提供内func。 parms是一个向量或参数列表; （可选）传递给函数的任何其他参数。 func的返回值应该是一个列表的第一个元素是一个向量包含y的衍生工具就time，而其下一个元素是全球性的值，其稳态值也是必需的。衍生工具应以相同的顺序指定的作为状态变量y。
参数：parms
parameters passed to func.
参数传递到func。
参数：nspec
the number of *species* (components) in the model. If NULL, then dimens should be specified.
的数目*物种*（组件）在模型中。如果NULL，那么dimens应该被指定。
参数：dimens
the number of *boxes* in the model. If NULL, then nspec should be specified.
在模型中的数目*盒*。如果为NULL，那么nspec应指定。
参数：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".
的溶液的方法，其中一个“stode”，“stodes”或“runsteady的”。
参数：cyclicBnd
if a cyclic boundary exists, a value of 1 else NULL; see details.
如果存在一个循环边界，1其他NULL;详情请参阅值。
参数：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".
相邻的框的传输发生的数目。通常等于1（盒我只相互作用箱的i-1和i +1）。大于1的值不会与method = "stodes"。
参数：...
additional arguments passed to the solver function as defined by method.
额外的参数传递给求解器的功能所定义的method。
Details
详细信息----------Details----------
This is the method of choice for multi-species 1-dimensional models, that are only subjected to transport between adjacent layers
这是于多品种的1维模型给出的方法，即只进行运送的相邻层之间
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))
A [1]，A [2]，A [3]，.... B [1]，B [2]，B [3]，...。 （物种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.
默认的方法重新排列状态变量A [1]，B [1]，... A [2]，B [2]，... A [3]，B [3]，...此重新导致一个带状的雅可比矩阵（上部和下部）的一半带宽=物种数。然后函数stode解决了条带状的问题。
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.
第二种方法是使用功能stodes。基于尺寸的问题，该方法首先计算出的雅可比矩阵的稀疏性图案，运输只发生在相邻层之间的假设下。然后stodes被调用来解决这个问题。
As stodes is used to estimate steady-state, it may be necessary to specify the length of the real work array, lrw.
stodes是用来估计的稳定状态，可能有必要指定长度的真正的工作数组，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.
虽然一个合理的猜测lrw制成，它是可能的，这将是太低了。在这种情况下，steady.1D将返回一个错误消息，告知实际需要的工作数组的大小。在第二次尝试，然后lrw等于这个数字。
For single-species 1-D models, use steady.band.
对于单种1-D模型，使用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
如果被安排作为状态变量（例如A [1]，B [1]，A [2]，B [2]，A [3]，B [3]，...那么该模型应解决的 X>
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.
在某些情况下，环状的边界条件存在。这是时的第一个框的交互，最后一个盒子，反之亦然。在这种情况下，会有额外的非零条纹中的雅可比矩阵，需要加以考虑。循环边界的发生，可以打开通过指定参数cyclicBnd=1。如果是这种情况，那么稳定状态将被估计使用stodes。默认情况下是没有循环的界限。
值----------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*.
names如果不：一个向量，其在估算过程中的稳定状态方程的系统状态变量值的最后一次迭代。 names如果，与一列的矩阵，每一个稳态组件。
参数：...
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.
的输出属性steady，该函数返回TRUE，如果已达到稳定状态和属性precis在每次迭代中达到的精度。
注意----------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))
这是可取的，但不是强制性的指定两个nspec和dimens。在这种情况下，求解器可以检查输入是否是有道理的（也就是说，如果nspec *梦诗=长度（Y））
（作者）----------Author(s)----------
Karline Soetaert <karline.soetaert@nioz.nl>
参见----------See Also----------
plot.steady1D for plotting the output of steady.1D
plot.steady1D策划的输出steady.1D
steady, for a general interface to most of the steady-state solvers
steady，一个通用的接口，大部分的稳态解算器
steady.band, to find the steady-state of ODE models with a banded Jacobian
steady.band，找到稳定状态的ODE模型的一个带状的雅可比
steady.2D, steady.3D, steady-state solvers for 2-D and 3-D partial differential equations.
steady.2D，steady.3D，2-D和3-D偏微分方程的稳态解算器。
stode, iterative steady-state solver for ODEs with full or banded Jacobian.
stode，常微分方程全部或带状雅可比迭代稳态解算器。
stodes, iterative steady-state solver for ODEs with arbitrary sparse Jacobian.
stodes，常微分方程任意稀疏Jacobian迭代稳态解算器。
runsteady, steady-state solver by dynamically running to steady-state
runsteady，通过动态运行到稳定状态的稳态求解器
实例----------Examples----------
## =======================================================================[＃================================================= ======================]
## EXAMPLE 1: BOD + O2 [例1：BOD + O2]
## =======================================================================[＃================================================= ======================]
## Biochemical Oxygen Demand (BOD) and oxygen (O2) dynamics[＃生化需氧量（BOD）和氧气（O2）动态]
## in a river[＃在河]
#==================#[==================＃]
# Model equations #[模型方程＃]
#==================#[==================＃]
O2BOD <- function(t, state, pars) {
BOD <- state[1:N]
O2 <- state[(N+1) 2*N)]
# BOD dynamics[BOD动态]
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[1，莫诺次消费中的氧]
#rate of change = flux gradient - consumption + reaeration (O2)[变化率=通量梯度 - 消费+复氧（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[END O2BOD]
#==================#[==================＃]
# Model application#[典型应用＃]
#==================#[==================＃]
# parameters[参数]
dx <- 100 # grid size, meters[网格大小，米]
v <- 1e2 # velocity, m/day[速度，m /天]
x <- seq(dx/2, 10000, by = dx) # m, distance from river[米，距离河]
N <- length(x)
r <- 0.1 # /day, first-order decay of BOD[/天，一阶衰减的BOD]
p <- 0.1 # /day, air-sea exchange rate[/天，海 - 气交换率]
O2sat <- 300 # mmol/m3 saturated oxygen conc[mmol/m3饱和的氧气浓度]
O2_0 <- 50 # mmol/m3 riverine oxygen conc[mmol/m3河流氧浓]
BOD_0 <- 1500 # mmol/m3 riverine BOD concentration[mmol/m3河流BOD浓度]
# 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":[一次绘制所有的状态变量，对“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 [例2：硅酸盐成岩作用]
## =======================================================================[＃================================================= ======================]
## Example from the book:[＃从书上的例子：]
## Soetaert and Herman (2009).[＃Soetaert和赫尔曼（2009年）。]
## a practical guide to ecological modelling -[＃生态模型的实用指南 - ]
## using R as a simulation platform.[使用R仿真平台。]
## Springer[＃施普林格]
#====================#[====================＃]
# Model equations #[模型方程＃]
#====================#[====================＃]
SiDIAmodel <- function (time = 0, # time, not used here[时间，而不是用在这里]
Conc, # concentrations: BSi, DSi[浓度：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[施加的上限浓度（bwDSi），较低的零梯度]
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 [BSI解散]
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 [DSi的沉积物 - 水交换率]
DSIDeepFlux = DSiFlux[N+1], # DSi deep-water (burial) flux[DSi的深水（埋葬）光通量]
BSiDeepFlux = BSiFlux[N+1])) # BSi deep-water (burial) flux[BSI深水（埋葬）光通量]
}
#====================#[====================＃]
# 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[cm2/day - 生物扰动系数]
dBcoeff <- 2
mixdepth <- 5 # cm[厘米]
Db <- pmin(dB0, dB0*exp(-(Intdepth-mixdepth)*dBcoeff))
# biogeochemical parameters[生物地球化学参数]
SedDisp <- 0.4 # diffusion coefficient, cm2/d[扩散系数，平方厘米/ D]
rDissSi <- 0.005 # dissolution rate, /day[溶解率，/天]
EquilSi <- 800 # equilibrium concentration[平衡浓度]
pow <- 1
BSidepo <- 0.2*100 # nmol/cm2/day[nmol/cm2/day]
bwDSi <- 150 # mmol/m3[mmol/m3]
# initial guess of state variables-just random numbers between 0,1[初始猜测的状态变量的随机数在0.1]
Conc <- runif(2*N)
# three runs with different deposition rates [三个不同的沉积速率运行]
BSidepo <- 0.2*100 # nmol/cm2/day[nmol/cm2/day]
sol <- steady.1D (Conc, func = SiDIAmodel, parms = NULL, nspec = 2,
names = c("DSi", "BSi"))
BSidepo <- 2*100 # nmol/cm2/day[nmol/cm2/day]
sol2 <- steady.1D (Conc, func = SiDIAmodel, parms = NULL, nspec = 2,
names = c("DSi", "BSi"))
BSidepo <- 3*100 # nmol/cm2/day[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 [图3运行]
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 = " orosity", 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)
转载请注明:出自 生物统计家园网(http://www.biostatistic.net)。
注：
注1：为了方便大家学习，本文档为生物统计家园网机器人LoveR翻译而成，仅供个人R语言学习参考使用，生物统计家园保留版权。
注2：由于是机器人自动翻译，难免有不准确之处，使用时仔细对照中、英文内容进行反复理解，可以帮助R语言的学习。
注3：如遇到不准确之处，请在本贴的后面进行回帖，我们会逐渐进行修订。
|
发表于 2012-9-27 23:00:19
