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

廿TT

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

タンクモデルのシミュレーション

R 微分方程式

タンクモデルとは

を参照。

f:id:abrahamcow:20150913000846p:plain
(上記ページから引用)

シンプルといえばシンプルだけどこんなモデルどうやってパラメータ推定するんだろう。

都市域:1段タンクモデル

#R のコード
library(deSolve)
library(tidyr)
library(ggplot2)
library(pipeR)

tunk1 <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dS1 = (1-beta1)*S1  + R
    return(list(c(dS1)))
  })
}
times <- seq(0, 5, by = 0.1)
pars  <- c(beta1=5,R=50)
ini  <- c(S1 = 20)
out1   <- ode(ini,times, tunk1, pars)

S1;タンクの貯留高
β1:タンクの浸透流出孔の浸透係数
R:単位時間あたりの降雨量

dS1 = S1 - beta1*S1 + R
(タンクの貯留高の増加量)=(タンクの貯留高)−(タンクの浸透流出孔からの流出)+(単位時間あたりの降雨量)

#プロット
as.data.frame(out1) %>>%
  gather(tunk,value,-time) %>>%
  ggplot(aes(time,value))+
  geom_line()

f:id:abrahamcow:20150913002044p:plain

非都市域:直列3段タンクモデル

弱雨時

tunk2 <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dS1 = (1-beta1)*S1  + R
    dS2 = (1-beta2)*S2  + beta1*S1
    dS3 = (1-beta3)*S3  + beta2*S2 - q3
    return(list(c(dS1,dS2,dS3)))
  })
}
pars  <- c(beta1=5,beta2=4,beta3=3,R=50,q3=30)
ini  <- c(S1 = 20, S2=10, S3=0)
out2   <- ode(ini,times, tunk2, pars)

f:id:abrahamcow:20150913002106p:plain

強雨時

tunk3 <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dS1 = (1-beta1)*S1 - q1 + R
    dS2 = (1-beta2)*S2 - q2 + beta1*S1
    dS3 = (1-beta3)*S3 - q3 + beta2*S2
    return(list(c(dS1,dS2,dS3)))
  })
}
pars  <- c(beta1=5,beta2=4,beta3=3,q1=10,q2=20,q3=30,R=50)
ini  <- c(S1 = 20, S2=10, S3=0)
out3   <- ode(ini,times, tunk3, pars)

as.data.frame(out3) %>>%
  gather(tunk,value,-time) %>>%
  ggplot(aes(time,value))+
  geom_line()+
  facet_wrap(~tunk)

f:id:abrahamcow:20150913002216p:plain