r - 微分方程式の近似：データセットをRの微分方程式に近似する方法

Question

データセットの場合：

conc <- data.frame(time = c(0.16, 0.5, 1.0, 1.5, 2.0, 2.5, 3), concentration = c(170, 122, 74, 45, 28, 17, 10))

そして、このデータを以下の微分方程式に当てはめたいと思います。

dC/dt= -kC

ここで、Cは濃度、データセットの時間tです。これにより、kの結果も得られます。誰かが私にRでこれを行う方法の手がかりを教えてもらえますか？ありがとう。

Answer 1

まず、変数分離を使用して微分方程式を解きます。これにより、log（C）=-k * t+C0が得られます。

データをプロットします。

plot(log(concentration) ~ time,data=conc)

線形モデルを適合させる：

fit <- lm(log(concentration) ~ time,data=conc)
summary(fit)

# Coefficients:
#               Estimate Std. Error t value Pr(>|t|)    
# (Intercept)   5.299355   0.009787   541.4 4.08e-13 ***
#   time       -0.992208   0.005426  -182.9 9.28e-11 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
# 
# Residual standard error: 0.01388 on 5 degrees of freedom
# Multiple R-squared: 0.9999,  Adjusted R-squared: 0.9998 
# F-statistic: 3.344e+04 on 1 and 5 DF,  p-value: 9.281e-11 

予測値をプロットします。

lines(predict(fit)~conc$time)

抽出k：

k <- -coef(fit)[2]
#0.9922081

Answer 2

これは解決策になる可能性があります：

    require('deSolve')
conc <- data.frame(time <- c(0.16, 0.5, 1.0, 1.5, 2.0, 2.5, 3), concentration <- c(170, 122, 74, 45, 28, 17, 10))

##"Model" with differential equation
model <- function(t, C, k){
  list(-k * C)
}

##Cost function with sum of squared residuals:
cost <- function(k){
  c.start <- 170
  out <- lsoda(c(C=c.start), conc$time, model, k)
  c <- sum( (conc$concentration - out[,"C"])^2)       
  c
}

##Initial value for k
k <- 3
##  Use some optimization procedure
opt <- optim(k, cost, method="Brent", lower=0, upper=10)

k.fitted <- opt$par

lsodaを使用することは、 1つの微分方程式だけで計算するには少しやり過ぎのように思われるので、多分それは少しナイーブです...しかし、それは確かにあなたのkを最適化します。統合のためにCの開始値を確認することをお勧めします。ここでは、170に設定しました。t= 0の値はありませんか？