読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

メモ:可逆反応の微分方程式

微分方程式 数値計算 R

可逆反応微分方程式

反応物 A が反応物 B を生成し, 生成物から反応物に戻る反応もあるとする.

微分方程式,

\displaystyle  \frac{dA}{dt}=-k_1 A + k_2 B

\displaystyle \frac{dB}{dt} =k_1A - k_2B

は解,

\displaystyle  A(t)  =  A  _{0}\frac{1}{k_{1}+k_{2}}\left( k_{2}+k_{ 1}e^{-\left( k_{ 1}+k_{2} \right)t} \right)+  B _{0} \frac{k_{2}}{k_{1}+k_{2}}\left( 1-e^{-\left( k_{ 1}+k_{2} \right)t} \right)

\displaystyle B(t) = A _{0}\frac{k_{1}}{k_{1}+k_{2}}\left( 1-e^{-\left( k_{1}+k_{2} \right)t} \right)+ B _{0}\frac{1}{k_{1}+k_{2}}\left( k_{1}+k_{2}e^{-\left( k_{1}+k_{2} \right)t} \right)

を持つ.

R で計算

実線が解析解, 点が数値解.

f:id:abrahamcow:20150307081712p:plain

rxnrate=function(t,c,parms){
  
  # rate constant passed through a list called parms
  k1=parms$k1
  k2=parms$k2
  
  # c is the concentration of species
  
  # derivatives dc/dt are computed below
  r=rep(0,length(c))
  r[1]=-k1*c["A"] + k2*c["B"]
  r[2]=k1*c["A"] - k2*c["B"] 
  
  # the computed derivatives are returned as a list
  # order of derivatives needs to be the same as the order of species in c
  return(list(r))
  
}
cinit=c(A=1,B=0)
parms=list(k1=1,k2=0.1)

times <- seq(0, 10, by=0.1)

out=ode(y=cinit,times=times,func=rxnrate,parms=parms)

matplot(out[,1],out[,-1],pch=16,lty=1,lwd=2,xlab="time",ylab="concentration",
        col=c("royalblue","orange2"))
Afn <- function(t,k1,k2,A0,B0){
  ( A0/(k1+k2) )*(k2+k1*exp(-(k1+k2)*t)) + B0*(k2/(k1+k2))*(1-exp(-(k1+k2)*t)) 
}

Bfn <- function(t,k1,k2,A0,B0){
  ( A0 * (k1/(k1+k2)) )*(1- exp(-(k1+k2)*t)) + (B0/(k1+k2))*(k1+exp(-(k1+k2)*t)) 
}

curve(Afn(t=x,A0=cinit["A"],B0=cinit["B"],k1=parms$k1,k2=parms$k2),add=TRUE,col="royalblue",lwd=3)
curve(Bfn(t=x,A0=cinit["A"],B0=cinit["B"],k1=parms$k1,k2=parms$k2),add=TRUE,col="orange2",lwd=3)

legend(locator(1),legend=c("A","B"),lty=1,lwd=3,pch=16,
       col=c("royalblue","orange2"),bty = "n")