I want to obtain new values for each step and use them inside an equation for further steps. Normally, I can perform loop but for this problem, I need to use past values, too. For instance, I have flow data like this:

q<-c(10, 15.83333, 21.66667)

I created a loop manually:

    z1<-190 #initial elevation
      s1<-24011 #initial storage
      in1<-q[1] #initial inflow
      out1<-1.86*sqrt((z1-110)*19.62) #outflow
      z2<-  z1 0.3*((in1-out1)/s1) #elevation at second step
      
      in2<-q[2] #second inflow
      out2<-1.86*sqrt((z2-110)*19.62) #outflow at z2 elevation
    ds2<-0.3*((in1 in2)/2-(out1 out2)/2) #change in storage 
     s2<-s1 ds2 #net storage value
     z3<-z2 0.3*((in2-out2)/s2) #elevation at third step
      
     in3<-q[3]
     out3<-1.86*sqrt((z3-110)*19.62)
     ds3<-0.3*((in2 in3)/2-(out2 out3)/2)
     s3<-s2 ds3 
.
.
.
     z4<-z3 0.3*((in3-out3)/s3) 

Briefly, I am calculating z value using previous values of in,out,s. What I need to find is z values considering q values.

Expected result is:

        z        q     outflows  storages
[1,] 190.0000 10.00000 73.68981 24011.00
[2,] 189.9992 15.83333 73.68944 23992.77
[3,] 189.9985 21.66667 73.68911 23976.29

CodePudding user response：

1. Init values and create a result table
2. Add the current state values to the table
3. Simulate new state using old or new states
4. Set the new state to all variables

library(tidyverse)

data <- tibble(step = numeric(), out = numeric(), y = numeric(), z = numeric())

# Initialization
z <- 190
y <- 1
out <- NA

for (step in seq(5)) {
  # save current state
  data <- data %>% add_row(step = step, out = out, z = z, y = y)

  # use old state of z
  new_out <- z / 2
  
  # use old state of y
  new_z <- y   1 
  
  # use new state of out
  new_y <- new_out
  
  # Lastly, update all new variables
  out <- new_out
  y <- new_y
  z <- new_z
}
data
#> # A tibble: 5 x 4
#>    step   out     y     z
#>   <dbl> <dbl> <dbl> <dbl>
#> 1     1    NA     1   190
#> 2     2    95    95     2
#> 3     3     1     1    96
#> 4     4    48    48     2
#> 5     5     1     1    49

^{Created on 2021-11-10 by the reprex package (v2.0.1)}

CodePudding user response：

I'll extend my comment here.

You can overwrite the variables. Following your implementation you could for instance create a temporal variable for your out2:

q = c(10, 15.83333, 21.66667)

#Results storage array
results = array(numeric(),c(length(q),4)) 
colnames(results) = c("z", "q", "outflows", "storages")

z = 190
s = 24011
infl = q[1]
out = 1.86*sqrt((z-110)*19.62)

#Save init values
results[1,1] = z
results[,2] = q
results[1,3] = out
results[1,4] = s

for (n in 2:length(q)) {

  z = z 0.3*((infl-out)/s)
  out_tmp = 1.86*sqrt((z-110)*19.62)
  ds = 0.3*((infl q[n])/2-(out out_tmp)/2)
  s = s ds
  infl = q[n] 
  out = out_tmp

  results[n,1] = z
  results[n,3] = out
  results[n,4] = s

}

View(results)

If you want to avoid to create the temporal variable, you can try something like this:

q = c(10, 15.83333, 21.66667)

results = array(numeric(),c(length(q),4)) 
colnames(results) = c("z", "q", "outflows", "storages")
z = 190
s = 24011
infl = q[1]
out = 1.86*sqrt((z-110)*19.62)

#Save init values
results[1,1] = z
results[,2] = q
results[1,3] = out
results[1,4] = s

for (n in 2:length(q)) {

  z = z 0.3*((infl-out)/s)
  out = 1.86*sqrt((z-110)*19.62)
  ds = 0.3*((infl q[n])/2-(results[n-1,3] out)/2)
  s = s ds
  infl = q[n] 

  results[n,1] = z
  results[n,3] = out
  results[n,4] = s

}

View(results)