The numerical integration error must decrease when time steps taken are decreased. Time step convergence means that the model solution converges to the truth, i.e. the analytical solution when it is known. When this solution is not known, it can still be checked if the model follows the expectation that when time-steps are decreased by gradually smaller time steps, the solution should also show a gradual convergence to one solution where large time steps must result in large errors and small time steps in small errors.

A model with a fixed time-step of e.g. 1, takes steps of one time unit. When relevant parameters are expressed in absolute (or relative) value (or amount) per year, then also the rate equations will be expressed in an amount per year. Euler integration adds the rate times the time-step (of 1 year in this case) to the state: a time step of 1 therefore represents a step of one year and a constant rate of change during this whole year. This also means that if these relevant parameters are expressed in absolute (or relative) values (or amounts) per second, one time step will be one second as well and the rates of change will only be kept constant for a second. This solution is obviously more accurate!

Note that a variable time step will automatically adjust its time-step, so then changing parameter values from years to second may change the calculation time (as it may take more or less iterations to find an acceptable time step), but not the actual time steps taken by the model as these are defined by an acceptable integration error in the computed new value of the state.

# Set AL1 here, to automatically change its value when used in the call to the ParmsForSpace function
setA11 <- 0.002 # 0.002 is the default

# Use deSolve to integrate derivatives numerically with Runge-Kutta and a variable time-step
states_ode45 <- ode(y = InitForSpace(),
                    time = seq(from = 0, to = 25000, by = 10),
                    func = RatesForSpace, 
                    parms = ParmsForSpace(A11 = setA11),  
                    method = "ode45")

# Define plotting settings
# Use the help function to understand what this function call using "par" does!
# Plot the states as function of time, with narrow plot margins.
par(mfrow = c(1,2), # nr of rows and columns, filled by row. defaults to c(1,1)
    mar = c(3, 3, 0.25, 0.25), # margins mar = c(bottom, left, top, right), default is c(5, 4, 4, 2) + 0.1
    mgp = c(1.5, 0.5, 0)) # margin line mgp = c(axes text, label, line) in mex units, default is c(3, 1, 0)

# Plot the  states to compare vegetation and herbivores
plot(states_ode45[,"time"], states_ode45[,"P1"], type = "l", col = "black",
     xlim = c(0,25000), ylim = c(0,700), xlab = "Time, d", ylab = "Plants, g biomass / m2")
lines(states_ode45[,"time"], states_ode45[,"P2"], type = "l", col = "red")
legend("topleft", lty = 1, col = c("black","red"), legend = c("P1","P2"), bty = "n")

plot(states_ode45[,"time"], states_ode45[,"H1"], type = "l", col = "black",
     xlim = c(0,25000), ylim = c(0,15), xlab = "Time, d", ylab = "Herbivores, g herbivore / m2")
lines(states_ode45[,"time"], states_ode45[,"H2"],type = "l",col = "red")
legend("topleft", lty = 1, col = c("black","red"), legend = c("H1","H2"), bty = "n")

Indeed, the model did reach an equilibrium. We can also inspect the mass balance:

# compute mass balance in grams
MBForSpace_ode45 <- MassBalanceForSpace(states_ode45, 
                                        ini = InitForSpace(), 
                                        parms = ParmsForSpace(A11 = setA11)) 

# plot mass balance balance after resetting plotting canvas
par(mfrow = c(1,1))
plot(MBForSpace_ode45[,"time"], MBForSpace_ode45[,"BAL"], type = "l", col = "black",
     xlab = "Time, d", ylab = "Balance, g biomass/m2", ylim = c(-5e-10, 5e-10))

The values of the states and the equilibrium is affected by the value for the parameter A11. The model did not yet reach an equilibrium within 25000 days when A11 increased to 0.004 or 0.008 as shown below. The number of herbivores 1 and plant species 1 decreased, while plant species 2 ended up at about the same value. The number of H2 strongly increased as shown below.

Parameter A11 set to value 0.004:

# Set AL1 here, to automatically change its value when used in the call to the ParmsForSpace function
setA11 <- 0.004 # 0.002 is the default

# Use deSolve to integrate derivatives numerically with Runge-Kutta and a variable time-step
states_ode45_004 <- ode(y = InitForSpace(),
                        time = seq(from = 0, to = 25000, by = 10),
                        func = RatesForSpace,
                        parms = ParmsForSpace(A11 = setA11),
                        method = "ode45")

# Define plotting settings
par(mfrow = c(1,2), # nr of rows and columns, filled by row. defaults to c(1,1)
    mar = c(3, 3, 0.25, 0.25), # margins mar = c(bottom, left, top, right), default is c(5, 4, 4, 2) + 0.1
    mgp = c(1.5, 0.5, 0)) # margin line mgp = c(axes text, label, line) in mex units, default is c(3, 1, 0)

# Plot the  states to compare vegetation and herbivores
plot(states_ode45_004[,"time"], states_ode45_004[,"P1"], type = "l", col = "black",
     xlim = c(0,25000), ylim = c(0,700), xlab = "Time, d", ylab = "Plants, g biomass / m2")
lines(states_ode45_004[,"time"], states_ode45_004[,"P2"], type = "l", col = "red")
legend("topleft", lty = 1, col = c("black","red"), legend = c("P1","P2"), bty = "n")

plot(states_ode45_004[,"time"], states_ode45_004[,"H1"], type = "l", col = "black",
     xlim = c(0,25000), ylim = c(0,15), xlab = "Time, d", ylab = "Herbivores, g herbivore / m2")
lines(states_ode45_004[,"time"], states_ode45_004[,"H2"],type = "l",col = "red")
legend("topleft", lty = 1, col = c("black","red"), legend = c("H1","H2"), bty = "n")

Parameter A11 set to value 0.008:

# Set AL1 here, to automatically change its value when used in the call to the ParmsForSpace function
setA11 <- 0.008 # 0.002 is the default

# Use deSolve to integrate derivatives numerically with Runge-Kutta and a variable time-step
states_ode45_008 <- ode(y = InitForSpace(),
                        time = seq(from = 0, to = 25000, by = 10),
                        func = RatesForSpace,
                        parms = ParmsForSpace(A11 = setA11),
                        method = "ode45")

# Define plotting settings
par(mfrow = c(1,2), # nr of rows and columns, filled by row. defaults to c(1,1)
    mar = c(3, 3, 0.25, 0.25), # margins mar = c(bottom, left, top, right), default is c(5, 4, 4, 2) + 0.1
    mgp = c(1.5, 0.5, 0)) # margin line mgp = c(axes text, label, line) in mex units, default is c(3, 1, 0)

# Plot the  states to compare vegetation and herbivores
plot(states_ode45_008[,"time"], states_ode45_008[,"P1"], type = "l", col = "black",
     xlim = c(0,25000), ylim = c(0,700), xlab = "Time, d", ylab = "Plants, g biomass / m2")
lines(states_ode45_008[,"time"], states_ode45_008[,"P2"], type = "l", col = "red")
legend("topleft", lty = 1, col = c("black","red"), legend = c("P1","P2"), bty = "n")

plot(states_ode45_008[,"time"], states_ode45_008[,"H1"], type = "l", col = "black",
     xlim = c(0,25000), ylim = c(0,15), xlab = "Time, d", ylab = "Herbivores, g herbivore / m2")
lines(states_ode45_008[,"time"], states_ode45_008[,"H2"],type = "l",col = "red")
legend("topleft", lty = 1, col = c("black","red"), legend = c("H1","H2"), bty = "n")

# Reset plotting settings (to their default values)
par(mfrow = c(1,1),
    mar = c(5, 4, 4, 2) + 0.1,
    mgp = c(3, 1, 0)) 

Increasing the parameter $A11$ increased the preference for $P1$ by $H1$. Therefore more of $P2$ was available for $H2$. $P1$ is only eaten by $H1$ and this increase of the A11 parameter increased the pressure on $P1$ by $H1$ grazing. The reduced pressure on $P2$ from $H1$ grazing was compensated by more $P2$ grazing. Note that herbivore 2 is more effective in searching for this $P2$ when compared to $H1$ as the parameter value for $A12$ is half the value of parameter $A22$.

Comparison of ode45 and Euler with a time step of 1, 10 and 100 days. Indeed, with smaller time step during Euler integration, the model is converging with solutions ever closer to the “ode45” solution. This is expected, as we expect that errors are bigger if the rates are kept constant for longer time steps, while in reality they change continuously, even at infinitely small timesteps. Note that also the “ode45” solution will have some error when compared to a analytic solution, although this error will be very small.

# Set parameter A11 to default value
setA11 <- 0.002 # 0.002 is the default

# Solve with Euler and varyiing time steps
states_euler1 <- ode(y = InitForSpace(),
                    time = seq(from = 0, to = 2500, by = 1),
                    func = RatesForSpace, 
                    parms = ParmsForSpace(A11 = setA11),  
                    method = "euler")
states_euler10 <- ode(y = InitForSpace(),
                    time = seq(from = 0, to = 2500, by = 10),
                    func = RatesForSpace, 
                    parms = ParmsForSpace(A11 = setA11),  
                    method = "euler")
states_euler100 <- ode(y = InitForSpace(),
                    time = seq(from = 0, to = 2500, by = 100),
                    func = RatesForSpace, 
                    parms = ParmsForSpace(A11 = setA11),  
                    method = "euler")

# Define plotting settings
# Use the help function to understand what this function call using "par" does!
# Plot the states as function of time, with narrow plot margins.
par(mfrow = c(1,2), # nr of rows and columns, filled by row. defaults to c(1,1)
    mar = c(3, 3, 0.25, 0.25), # margins mar = c(bottom, left, top, right), default is c(5, 4, 4, 2) + 0.1
    mgp = c(1.5, 0.5, 0)) # margin line mgp = c(axes text, label, line) in mex units, default is c(3, 1, 0)

# Plot the  states to compare vegetation and herbivores
plot(states_ode45[,"time"], states_ode45[,"P1"], type = "l", col = "black",
     xlim = c(0,2500), ylim = c(0,800), xlab = "Time, d", ylab = "Plants, g biomass / m2")
lines(states_ode45[,"time"], states_ode45[,"P2"], type = "l", col = "red")
# Add euler solutions with different line types (but same colours)
# 1 time unit
lines(states_euler1[,"time"], states_euler1[,"P1"], type = "l", col = "black", lty = 2)
lines(states_euler1[,"time"], states_euler1[,"P2"], type = "l", col = "red", lty = 2)
# 10 time units
lines(states_euler10[,"time"], states_euler10[,"P1"], type = "l", col = "black", lty = 3)
lines(states_euler10[,"time"], states_euler10[,"P2"], type = "l", col = "red", lty = 3)
# 100 time unit
lines(states_euler100[,"time"], states_euler100[,"P1"], type = "l", col = "black", lty = 4)
lines(states_euler100[,"time"], states_euler100[,"P2"], type = "l", col = "red", lty = 4)
# Add legend (colour for state, linetype for method)
legend("topleft", title = "Solver", col = "grey", 
       legend = c("ode45","Euler 1", "Euler 10", "Euler 100"), lty = 1:4, bty = "n")
legend("topright", title = "State", col = c("black","red"), 
       legend = c("P1", "P2"), lty = 1, bty = "n")

# Plot the  states to compare vegetation and herbivores
plot(states_ode45[,"time"], states_ode45[,"H1"], type = "l", col = "black",
     xlim = c(0,2500), ylim = c(0,15), xlab = "Time, d", ylab = "Plants, g biomass / m2")
lines(states_ode45[,"time"], states_ode45[,"H2"], type = "l", col = "red")
# Add euler solutions with different line types (but same colours)
# 1 time unit
lines(states_euler1[,"time"], states_euler1[,"H1"], type = "l", col = "black", lty = 2)
lines(states_euler1[,"time"], states_euler1[,"H2"], type = "l", col = "red", lty = 2)
# 10 time units
lines(states_euler10[,"time"], states_euler10[,"H1"], type = "l", col = "black", lty = 3)
lines(states_euler10[,"time"], states_euler10[,"H2"], type = "l", col = "red", lty = 3)
# 100 time unit
lines(states_euler100[,"time"], states_euler100[,"H1"], type = "l", col = "black", lty = 4)
lines(states_euler100[,"time"], states_euler100[,"H2"], type = "l", col = "red", lty = 4)
# Add legend (colour for state, linetype for method)
legend("topleft", title = "Solver", col = "grey", 
       legend = c("ode45","Euler 1", "Euler 10", "Euler 100"), lty = 1:4, bty = "n")
legend("topright", title = "State", col = c("black","red"), 
       legend = c("H1", "H2"), lty = 1, bty = "n")

When you set inputs to 0 the rate equations for $INF$, $IF$, $SICK$, $HOSP$ and $IM$ are reflecting a situation with exponential decrease. For $UP$ however this is not the case, here outflow depends on the number of $IF$ persons, while for $DEAD$ there is no outflow.

The TC is a system characteristic and does not depend on the inflow. With input set to 0, equilibrium values are 0 for $INF$, $IF$, $SICK$, $HOSP$ and $IM$. The inputs in the rate equations then also become 0 and fall away in these equations.

For $IF$: $TC=1/rinfdr$, as $TC=|\frac{IFeq - IF}{rinfdr * IF}|$ and $TC = \frac{IF}{rinfdr*IF}=\frac{1}{rinfdr}$.

This approach works as well for the other states.

For $INF$: $TC=\frac{1}{rinfr}$

For $IM$: $TC=\frac{1}{rdimr}$

For $SICK$, without inflow the rate equation simplifies rate equation to:

$\frac{dSICK}{dt} = -(rdrs + rsrh + rims) * SICK$

$TC = |\frac{1}{rdrs + rsrh + rims} |$

Likewise for $HOSP$: $TC=\frac{1}{rdrh+rimh}$.

For exponential equations, ART = TC and with equations from question 9.3b for IF, ART = 1/(1/300) = 300 days.

For exponential equations, ART = TC. With an ART of 3, TC = 3 and rinfdr = 1/3 d^-1

For $HOSP$ in equilibrium the state $SICK$ also needs to be in equilibrium:

$$\begin{aligned} \frac{}{}dHOSPdt &= rsrh * SICK - rdrh * HOSP - rimh * HOSP = 0 \\ - rsrh * SICKeq &= - (rdrh + rimh) * HOSPeq \\ HOSPeq &= \frac{rsrh * SICKeq}{rdrh + rimh} \\ \end{aligned}$$

For $SICK$:

$$\begin{aligned} \frac{dSICK}{dt} &= rinfdr_s * IF - rdrs * SICK - rsrh * SICK - rims * SICK = 0 \\ - rinfdr_s * IF_{eq} &= - rdrs *SICK_{eq} - rsrh*SICK_{eq} - rims *SICK_{eq} \\ rinfdr_s * IF_{eq} &= (rdrs + rsrh + rims) * SICK_{eq} \\ SICK_{eq} &= \frac{rinfdr_s * IF_{eq}}{rdrs + rsrh + rims} \end{aligned}$$

The number of people that will be infected by 1 person can be determined by looking at a component of the rate equation for $IF$. Note that there was only one person infectious that remains infectious for a short period of time. The average residence time for $IF = 3$, so a person is infectious for on average only 3 days (see 9.3d). Note that the parameter $rinfdr$ had a value of 1/3. With one infected person $IF$ equals 1 and we can use this rate equation:

$$\frac{dINF}{dt} = curr_{nen} * rifr * IF$$

The total number of people that will be infected by this person, the R number, can be found by integrating this rate equation from 0 to 3 days:

$$R = \int_0^3 curr_{nen} × rifr \text{ } dt = | curr_{nen} × rifr × t |_0^3 = curr_{nen} × rifr × 3$$

It can be easily understood that with e.g. 60 person to person encounters per day and a relative infection rate of 0.01 each day 0.6 other persons are infected resulting in an R value of 1.8. A new virus strain with an $rifr$ of 0.02 will result in an R value of 60 * 0.02 * 3 = 3.6.