方程的求根大家应该在高中就已深入骨髓,今天给大家介绍下在R语言中如何实现方程的求根以及方程中参数的确定。我们需要借助R包rootSolve开始我们的教程。包的安装就不再赘述了。直接进入主题,在此包中求根函数涉及以下三个:
1. uniroot 求一个方程式的一个根。
其中主要的参数是extendInt,在指定范围无法获得相对应的根时,可以利用此参数进行扩展x上下限从而得到对应的根。默认是no,如果我们不确定我们的根是否在我们填写的范围,那么可以设置本参数为yes。
我们来看下实例:
代码语言:javascript复制f1 <- function(x) (121 - x^2)/(x^2 1)
f2 <- function(x) exp(-x)*(x - 12)
try(uniroot(f1, c(0,10)))
try(uniroot(f2, c(0, 2)))
代码语言:javascript复制u1 <- uniroot(f1, c(0,10),extendInt="yes", trace=1)
u2 <- uniroot(f2, c(0,2), extendInt="yes", trace=2)
2. polyroot 求一个带有复杂参数多项式的根。比如f(x)=z1 z2*x z3*x^2 … zn*x^(n-1)。
只需要将所有的参数按照x的排序逐渐增大的顺序将参数以次输入函数即可得到对应的f(x)=0的根。
3. multiroot 求多个方程集合的根。
其中主要的就是model构建以及parms参数的输入。
来看下下面的实例:
代码语言:javascript复制#构建模型
model2<- function(x, parms)
c(F1 = x[1] x[2] x[3]^2 - parms[1],
F2 = x[1]^2 - x[2] x[3] - parms[2],
F3 = 2 * x[1] - x[2]^2 x[3] -parms[3])
# 求解
parms<- c(12, 2, 1)
multiroot(model2,c(1, 1, 1), parms = parms)
代码语言:javascript复制#处理parms求解
multiroot(model2, c(0, 0, 0), parms = parms*2)
以上结果中:
Root指的方程根的位置的x的值;
f.root指的方程根位置的y值。
Iter指的迭代的次数。
Estim.pricis指的评估的精准度,主要指的f.root的绝对值的均值
在后期的更新之后此包具有了获取系数值的功能,涉及两个函数
4. runsteady 动态计算参数直至其导数不再发生变化。
具体的应用直接看实例:
代码语言:javascript复制#构建模型
model<-function(t,y, pars) {
with(as.list(c(y, pars)),{
Min = r*OM
oxicmin = Min*(O2/(O2 ks))
anoxicmin = Min*(1-O2/(O2 ks))* SO4/(SO4 ks2)
dOM =Flux - oxicmin - anoxicmin
dO2 =-oxicmin -2*rox*HS*(O2/(O2 ks)) D*(BO2-O2)
dSO4 = -0.5*anoxicmin rox*HS*(O2/(O2 ks)) D*(BSO4-SO4)
dHS =0.5*anoxicmin -rox*HS*(O2/(O2 ks)) D*(BHS-HS)
list(c(dOM, dO2, dSO4, dHS), SumS = SO4 HS)
})
}
#参数的设置
pars<- c(D = 1, Flux = 100, r = 0.1, rox = 1,
ks = 1, ks2 = 1, BO2 = 100, BSO4 = 10000, BHS= 0)
#需要计算的参数
y <- c(OM = 1, O2 = 1, SO4 = 1, HS = 1)
#计算
print(system.time(
ST2 <- runsteady(y = y, func = model,parms = pars, times = c(0, 1000))
))
stode, stodes, steady,steady.1D, steady.2D, steady.3D, 和 steady.band函数都是基于牛顿迭代法获取对应的参数值。
其中stode和stodes的用法一样,其有时候和runsteady有同样的效果,区别就是为了适应生物学的方程stode和stodes函数有一个关键的参数pos可以保证参数的正性:
代码语言:javascript复制print(system.time(
ST <- stode(y = y, func = model, parms =pars, pos = TRUE)
))
有结果可以看到具有一定的一致性。
steady.1D, steady.2D, steady.3D主要是针对1维,2维,3维偏微分方程的求解。
这一块我也是一头雾水,那位大神懂的可以留言,让大家都膜拜膜拜