Рассчитать строки на основе цикла

Question

Я пытаюсь реализовать итерационный метод для расчета количества осадков при использовании уравнения инфильтрации Хортона. Код, описанный ниже, реализует первую строку (time = 0).

#Entrance
f0=6
f1=1
k=2
dt=0.25
f=0
time= seq(from=0, to=2, by=dt)
inc_rainfall=c(0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.4, 0.6, 0.6)
rainfall_int=inc_rainfall/dt

#fc
gfc = function(fc) {
  f - ((f0 - fc) / k) + ((f1 / k) * log((fc - f1) / (f0 - f1)))}
fc=uniroot(gfc, interval = c(f0, f1))
fc=fc$root

#F'
fl=if(fc>rainfall_int[1]) sum(f,inc_rainfall[1]) else 0

#fc'
gfcl = function(fcl) {
  fl - ((f0 - fcl) / k) + ((f1 / k) * log((fcl - f1) / (f0 - f1)))}
fcl=uniroot(gfcl, interval = c(f0, f1))
fcl=fcl$root

#Fp ou Fs
fporfs=ifelse(fc<rainfall_int[1],f,
       ifelse(fcl<rainfall_int[1],
              (f0-rainfall_int[1])/k-
               f1/k*log((rainfall_int[1]-f1)/(f0-f1)), 0))

#dt'
dtl=ifelse((fc>rainfall_int[1]), 
    ifelse(fcl<rainfall_int[1],
           (fporfs-f)/rainfall_int[1],0),0)

#ts
ts=ifelse(fc<rainfall_int[1],time[1],
   ifelse(fcl<rainfall_int[1],dtl+time[1], 0))

#to
hto = function(to) {
  fporfs-f1*(ts-to)-(f0-f1)/k*(1-exp(-k*(ts-to)))}
to=uniroot(hto, interval = c(0, 1))
to=to$root

#Ft+Dt
ftdt=ifelse(to==0, fl, f1*(time[1]-to)+(f0-f1)/k*(1-exp(-k*(time[1]-to)))) #This value will be the "f" on next row

#Infiltration
infiltr=ftdt+f

#Runoff in line 1
runoff1=inc_rainfall[1]-(ftdt)

#Runoff in line 2 to n
runoffn=inc_rainfall[1]-(ftdt[2]-ftdt[1])

out=as.data.frame(cbind(time[1], inc_rainfall, rainfall_int, f, fc, fl, fcl, 
                        fporfs, dtl, ts, to, ftdt, infiltr, runoff1))

colnames(out)= c("Time", "Incremental Rainfall", "Rainfall Intensity", "F", "fc", "Fl", "Fcl", "Fp or Fs", "dt", "ts", "to", "Ftdt", "Infiltration", "Runoff"    )

out

Как мне поступить, чтобы вычислить следующую строку, используя в качестве начального значения (f) значение ftdt предыдущей строки?

Обратите внимание, что в последнем столбце (runoff) также есть изменение функции первого ряда на другие. Помимо этого, вероятно, также потребуется добавить дополнительную точку time (2.25), чтобы вычислить переменную ftdt последней строки.

Ожидаемый результат можно увидеть здесь: http://hydrology.usu.edu/RRP/userdata/4/87/RainfallRunoffProcesses.pdf на стр. 109.

Answer 1

Это должно приблизить вас.Я поместил ваши уравнения в функцию, которую я использую в цикле for.Мне пришлось изменить некоторые ваши формулы, чтобы заставить его работать.Я ничего не знаю о гидрологии, но изменения имели смысл после взгляда на книгу, на которую вы ссылаетесь, и небольшую Википедию.Изменения также привели к получению данных, которые больше напоминают данные в таблице.Ищите #!!!!!!! <COMMENT> !!!!!!!#, чтобы увидеть, что я изменил.

На пятой итерации / строке формула ftdt=ifelse(to==0, fl, f1*(time-to)+(f0-f1)/k*(1-exp(-k*(time-to)))) дает неверный результат, что приводит к ошибкам в других результатах.Я заменил f1 на fc, потому что это имело смысл после прочтения соответствующих страниц книги, и результаты были ближе к тому, что вы хотите.Единственное, в чем я не уверен, это переменная времени t , которая в вашей формуле равна time - to, т.е. время минус смещение времени.Я думаю, что это может быть проблемой.

compute_horton <- function (f, time, inc_rainfall, rainfall_int) {
    #->->-> calculates runoff rainfall using the Horton infiltration equation

    # constants
    f0=6
    f1=1
    k=2

    #fc
    gfc = function(fc) {
        f - ((f0 - fc) / k) + ((f1 / k) * log((fc - f1) / (f0 - f1)))
    }
    fc=uniroot(gfc, interval = c(f0, f1))
    fc=fc$root

    #F'
    fl=if(fc>rainfall_int) sum(f,inc_rainfall) else 0

    #fc'
    gfcl = function(fcl) {
        fl - ((f0 - fcl) / k) + ((f1 / k) * log((fcl - f1) / (f0 - f1)))
    }
    fcl=uniroot(gfcl, interval = c(f0, f1))
    fcl=fcl$root

    #Fp ou Fs
    fporfs=ifelse(fc<rainfall_int,f,
                  ifelse(fcl<rainfall_int,
                         (f0-rainfall_int)/k-
                             f1/k*log((rainfall_int-f1)/(f0-f1)), 0))

    #dt'
    dtl=ifelse((fc>rainfall_int), 
               ifelse(fcl<rainfall_int,
                      (fporfs-f)/rainfall_int,0),0)

    #ts
    ts=ifelse(fc<rainfall_int,time,
              ifelse(fcl<rainfall_int,dtl+time, 0))

    #to
    hto = function(to) {
        fporfs-f1*(ts-to)-(f0-f1)/k*(1-exp(-k*(ts-to)))
        }
    to=uniroot(hto, interval = c(0, 1))
    to=to$root

    #Ft+Dt - This value will be the "f" on next row
    #!!!!!!! I THINK YOU NEED FC AND NOT F1 !!!!!!!#
    #!!!!!!! EVEN SO, FORMULA DOESN'T WORK QUITE RIGHT !!!!!!!#
    #ftdt=ifelse(to==0, fl, f1*(time-to)+(f0-f1)/k*(1-exp(-k*(time-to))))
    ftdt=ifelse(to==0, fl, fc*(time-to)+(f0-fc)/k*(1-exp(-k*(time-to))))

    #Infiltration
    #!!!!!!! I THINK InFILTRATION IS DIFFERENCE OF FTDT AND F !!!!!!!#
    infiltr=ftdt-f

    #Runoff in line 1
    #!!!!!!! I THINK RUNOFF IS DIFFERENCE OF RAINFALL AND INFILTRATION !!!!!!!#
    runoff1=round(inc_rainfall-infiltr, 3)

    #Runoff in line 2 to n
    #!!!!!!! THIS ISN'T USED ANYWHERE !!!!!!!#
    runoffn=inc_rainfall-(ftdt[2]-ftdt[1]) 

    #### OUTPUT ####
    c(time, inc_rainfall, rainfall_int, f, fc, fl, fcl, fporfs, dtl, ts, to,
      ftdt, infiltr, runoff1
      )
}


# Function input
f=0
dt=0.25
time= seq(from=0, to=2, by=dt)
inc_rainfall=c(0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.4, 0.6, 0.6)
rainfall_int=inc_rainfall/dt

# Matrix column names
cnames <- c("Time", "Incremental Rainfall", "Rainfall Intensity", "F", "fc",
           "Fl", "Fcl", "Fp or Fs", "dt", "ts", "to", "Ftdt", "Infiltration",
           "Runoff"
           )

# Initialize matrix and add column names
hort_mat <- matrix(0, nrow = length(inc_rainfall), ncol = length(cnames))
colnames(hort_mat) <- cnames

for (i in 1:nrow(hort_mat)) {
    hort_mat[i,] <- compute_horton(f = hort_mat[ifelse(i-1 > 0, i-1, i), "Ftdt"],
                                   time = time[i],
                                   inc_rainfall = inc_rainfall[i],
                                   rainfall_int = rainfall_int[i]
                                   )
}

Answer 2

Спасибо за ваш вклад, @ gersht .Основываясь на этом, я очистил код и вставил ifelse функции там, где их раньше не было.

Однако мне все еще трудно использовать вычисленное значение ftdt в первом ряду для запуска второго, наloop.

Следуйте установленному коду.Я думал о for для каждой строки, но я не могу придумать способ сделать это эффективно, сохраняя код чистым.

#Constants
f0=6
f1=1
k=2
dt=0.25
f=0
time= seq(from=0, to=2, by=dt)
inc_rainfall=c(0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.4, 0.6, 0.6)
rainfall_int=inc_rainfall/dt

horton=data.frame(matrix(ncol = 14, nrow = 9))
names(horton)=c("Time", "Incremental Rainfall", "Rainfall Intensity", "F", "fc",
                "Fl", "Fcl", "Fp or Fs", "dt", "ts", "to", "Ftdt", "Infiltration",
                "Runoff")

#g(Fc)
gfc=ifelse (f==0, 

  function(fc) {
  f - ((f0 - fc) / k) + ((f1 / k) * log((fc - f1) / (f0 - f1)))}, 

  function(fc) {
  ftdt[i-1] - ((f0 - fc) / k) + ((f1 / k) * log((fc - f1) / (f0 - f1)))})

fc=uniroot(gfc, interval = c(f0, f1))
fc=fc$root

#F'
fl=ifelse(fc>rainfall_int, f+inc_rainfall, 0)

#fc'
gfcl= ifelse (fl==0,

  0,

  function(fcl) {
  fl - ((f0 - fcl) / k) + ((f1 / k) * log((fcl - f1) / (f0 - f1)))})

fcl=uniroot(gfcl, interval = c(f0, f1))
fcl=fcl$root

#Fp or Fs
fporfs=ifelse(fc<rainfall_int[1],f,
              ifelse(fcl<rainfall_int[1],
                     (f0-rainfall_int[1])/k-
                       f1/k*log((rainfall_int[1]-f1)/(f0-f1)), 0))

#dt'
dtl=ifelse((fc>rainfall_int[1]), 
           ifelse(fcl<rainfall_int[1],
                  (fporfs-f)/rainfall_int[1],0),0)

#ts
ts=ifelse(fc<rainfall_int[1],time[1],
          ifelse(fcl<rainfall_int[1],dtl+time[1], 0))

#to
hto=function(to) {
  fporfs-f1*(ts-to)-(f0-f1)/k*(1-exp(-k*(ts-to)))}

to=uniroot(hto, interval = c(0, 1))
to=to$root

#Ft+Dt
ftdt=ifelse(to==0, fl, f1*(time-to)+(f0-f1)/k*(1-exp(-k*(time-to)))) #Value that starts next row at "gFc"

#Infiltration
infiltr=ifelse(to==0,ftdt-f, ftdt[i]-ftdt[i-1])

runoff=ifelse(time==0, inc_rainfall[i]-(ftdt),  inc_rainfall[i]-(ftdt[i-1]-ftdt[i]))


horton$Time=time
horton$`Incremental Rainfall`=inc_rainfall
horton$`Rainfall Intensity`=rainfall_int
horton$F[1]=f
horton$fc[1]=fc
horton$Fl[1]=fl
horton$Fcl[1]=fcl
horton$`Fp or Fs`[1]=fporfs
horton$dt[1]=dtl
horton$ts[1]=ts
horton$to[1]=to
horton$Ftdt[1]=ftdt
horton$Infiltration[1]=infiltr
horton$Runoff[1]=runoff

horton