flowchart LR
R[Growth rate r] --> G[Growth flow rP]
P[Population P] --> G
G --> Pflowchart LR R[Growth rate r] --> G[Growth flow rP] P[Population P] --> G G --> P
Learning goals
After working through this module, you should be able to:
- distinguish an amount, a rate of change, and a change over a time step;
- translate a differential equation into a difference equation;
- identify stocks, flows, parameters, initial conditions, and units;
- implement unconstrained growth or decay in R;
- compare an Euler simulation with an analytical solution; and
- explain the effects of time-step size and model assumptions.
Begin with a prediction
A population begins with 100 individuals and grows at a continuous rate of 10% per hour.
- After 8 hours, will it be less than 180, between 180 and 240, or greater than 240?
- Will it gain the same number of individuals every hour?
- What shape should its graph have?
Keep your predictions. We will test them after building the model.
From words to a model
Let \(P(t)\) be the population at time \(t\) and let \(r\) be its per-capita growth rate. Our central assumption is:
The population’s rate of change is proportional to the population currently present.
Mathematically,
\[ \frac{dP}{dt}=rP. \]
| Quantity | Meaning | Example units |
|---|---|---|
| \(P\) | current population | individuals |
| \(P_0\) | initial population | individuals |
| \(r\) | proportional growth rate | per hour |
| \(t\) | time | hours |
| \(\Delta t\) | length of one simulation step | hours |
The population is a stock: an amount accumulated over time. Growth is a flow into that stock. Because the flow \(rP\) depends on the stock itself, this is a positive feedback loop.
The Euler update
Over a short interval, approximate the derivative with
\[ \frac{P_{new}-P_{old}}{\Delta t}\approx rP_{old}. \]
Solving for the new population gives Euler’s update rule:
\[ P_{new}=P_{old}+rP_{old}\Delta t. \]
This has a form we will reuse throughout the course:
\[ \boxed{\text{new amount}=\text{old amount}+\text{rate of change}\times\text{time passed}} \]
Check one step by hand
With \(P_{old}=100\), \(r=0.10\) per hour, and \(\Delta t=1\) hour,
\[ P_{new}=100+(0.10)(100)(1)=110. \]
The next increase will be larger because the growth flow is calculated from a larger population.
Build the simulation in R
1. Define time, the parameter, and the initial condition
start_time <- 0 # hours
end_time <- 8 # hours
delta_t <- 1 # hours per step
growth_rate <- 0.10 # per hour
initial_population <- 100 # individualsKeeping units in the comments makes it easier to catch mistakes.
2. Create the simulated times
time <- seq(from = start_time, to = end_time, by = delta_t)
time[1] 0 1 2 3 4 5 6 7 8Time begins at 0, but R begins counting vector positions at 1. Therefore time[1] contains 0.
3. Create and initialize the stock
population <- numeric(length(time))
population[1] <- initial_populationThe vector will store the population’s entire trajectory. Initially, we know only its first value.
4. Advance through time
for (i in 2:length(time)) {
population_old <- population[i - 1]
growth <- growth_rate * population_old
change <- growth * delta_t
population[i] <- population_old + change
}Each iteration reads the previous stock, calculates its flow, converts the flow into a change, and stores the next stock value.
5. Organize and inspect the results
results <- data.frame(time, population)
results time population
1 0 100.0000
2 1 110.0000
3 2 121.0000
4 3 133.1000
5 4 146.4100
6 5 161.0510
7 6 177.1561
8 7 194.8717
9 8 214.3589Before trusting a graph, verify that the first two rows agree with the hand calculation: 100 at time 0 and 110 at time 1.
6. Graph the trajectory
plot(
results$time,
results$population,
type = "l",
lwd = 3,
col = "#8b3a3a",
xlab = "Time (hours)",
ylab = "Population (individuals)",
main = "Unconstrained population growth"
)
Return to your opening predictions. The curve bends upward because every increase creates a larger population—and therefore an even larger future increase.
Compare with the analytical solution
Calculus gives the exact solution to this differential equation:
\[ P(t)=P_0e^{rt}. \]
Calculate exact values at the same times and measure the numerical error:
results$exact <- initial_population * exp(growth_rate * results$time)
results$error <- results$population - results$exact
results time population exact error
1 0 100.0000 100.0000 0.0000000
2 1 110.0000 110.5171 -0.5170918
3 2 121.0000 122.1403 -1.1402758
4 3 133.1000 134.9859 -1.8858808
5 4 146.4100 149.1825 -2.7724698
6 5 161.0510 164.8721 -3.8211271
7 6 177.1561 182.2119 -5.0557800
8 7 194.8717 201.3753 -6.5035607
9 8 214.3589 222.5541 -8.1952118plot(
results$time,
results$exact,
type = "l",
lwd = 3,
col = "#214f73",
xlab = "Time (hours)",
ylab = "Population (individuals)",
main = "Euler simulation and analytical solution"
)
lines(results$time, results$population, lwd = 3, lty = 2, col = "#b65335")
legend(
"topleft",
legend = c("Analytical solution", "Euler simulation"),
col = c("#214f73", "#b65335"),
lty = c(1, 2),
lwd = 3,
bty = "n"
)
The code can correctly implement Euler’s method without producing the exact continuous solution. The difference is numerical error caused by replacing continuous change with discrete steps.
Refine the time step
To repeat the model easily, place the simulation inside a function:
simulate_growth <- function(initial_population, growth_rate,
end_time, delta_t) {
time <- seq(from = 0, to = end_time, by = delta_t)
population <- numeric(length(time))
population[1] <- initial_population
for (i in 2:length(time)) {
growth <- growth_rate * population[i - 1]
population[i] <- population[i - 1] + growth * delta_t
}
data.frame(time, population)
}Compare several time steps:
time_steps <- c(2, 1, 0.5, 0.25)
exact_at_end <- initial_population * exp(growth_rate * end_time)
final_population <- numeric(length(time_steps))
absolute_error <- numeric(length(time_steps))
for (i in seq_along(time_steps)) {
run <- simulate_growth(
initial_population,
growth_rate,
end_time,
time_steps[i]
)
final_population[i] <- tail(run$population, 1)
absolute_error[i] <- abs(final_population[i] - exact_at_end)
}
data.frame(
delta_t = time_steps,
number_of_steps = end_time / time_steps,
final_population,
absolute_error
) delta_t number_of_steps final_population absolute_error
1 2.00 4 207.3600 15.194093
2 1.00 8 214.3589 8.195212
3 0.50 16 218.2875 4.266634
4 0.25 32 220.3757 2.178399A smaller time step generally reduces Euler error because the rate is recalculated more often. It also requires more computation. Accuracy is therefore connected to both the numerical method and the time-step choice.
Unconstrained decay
The same equation describes decay when \(r<0\). No new algorithm is needed; changing the sign changes the system’s behavior.
decay <- simulate_growth(
initial_population = 100,
growth_rate = -0.10,
end_time = 20,
delta_t = 0.5
)
decay$exact <- 100 * exp(-0.10 * decay$time)plot(
decay$time,
decay$exact,
type = "l",
lwd = 3,
col = "#214f73",
xlab = "Time (hours)",
ylab = "Amount",
main = "Unconstrained decay"
)
lines(decay$time, decay$population, lwd = 3, lty = 2, col = "#b65335")
legend(
"topright",
legend = c("Analytical solution", "Euler simulation"),
col = c("#214f73", "#b65335"),
lty = c(1, 2),
lwd = 3,
bty = "n"
)
The exact decay curve approaches zero without reaching it. Euler’s method can produce an impossible negative amount if the time step is too large relative to the decay rate.
Assumptions and limits
This model assumes:
- the proportional rate remains constant;
- resources never become scarce;
- the population is well mixed;
- individuals are interchangeable;
- continuous rates adequately represent births and deaths; and
- outside influences are negligible.
Those assumptions may be useful over a limited interval. A numerically accurate answer can still be a poor prediction if the model’s assumptions do not fit reality.
Experiments
For each experiment, make a prediction, run the model, and explain what happened.
- Double the initial population while keeping the growth rate fixed.
- Double the growth rate while keeping the initial population fixed.
- Find a time step that keeps the absolute error at 8 hours below 1 individual.
- Find a decay time step that produces a negative amount. Explain why.
- Name one real process that might follow unconstrained growth or decay briefly, but not forever.
What to remember
- A differential equation describes an instantaneous rate of change.
- A difference equation gives an update a computer can execute.
- Euler’s pattern is
new = old + rate * time passed. - A stock’s trajectory can be stored in a vector.
- Smaller steps usually reduce Euler error but require more calculations.
- Verifying an implementation is different from validating a model.
Next: Project 1
In Saving for a Car, the same update pattern controls three quantities: savings, new-car price, and trade-in value.