reference <- 6.239
approximation <- 6.24
absolute_error <- abs(approximation - reference)
relative_error <- absolute_error / abs(reference)
absolute_error[1] 0.001100 * relative_error # percent[1] 0.01602821Meeting 1: Find and measure error
Start with a real failure
NASA’s Mars Climate Orbiter was lost after ground software used English units while onboard software used metric units. Before we discuss any floating-point arithmetic, ask: Would a more precise computer have fixed that problem? Why or why not?
Errors can enter at different stages:
| Type | Example | What to check |
|---|---|---|
| Data | An instrument reports the wrong temperature | Calibration and source quality |
| Model | We omit an important way a disease spreads | Assumptions and validation against observations |
| Implementation | One program supplies pounds while another expects newtons | Code, units, and interface tests |
| Numerical | A finite time step approximates continuous change | Precision and sensitivity to the method or step size |
The last category is our main coding focus, but fixing numerical error cannot repair a wrong model or incorrect units.
How large is an error?
Suppose the reference value is \(6.239\), but we report \(6.24\). Predict the difference before running R.
\[ \text{absolute error}=|\text{approximation}-\text{reference}|, \qquad \text{relative error}=\frac{\text{absolute error}}{|\text{reference}|}. \]
The relative error is undefined when the reference value is zero. Multiply it by 100 to express it as a percentage.
Which measure is easier to interpret if the number represents meters? Which helps compare errors for values measured on very different scales?
Why does 0.1 + 0.2 surprise us?
Predict whether the final line is TRUE or FALSE:
0.1 + 0.2[1] 0.30.3[1] 0.30.1 + 0.2 == 0.3[1] FALSER prints a rounded display. Ask it to show more digits:
sprintf("%.17f", 0.1 + 0.2)[1] "0.30000000000000004"sprintf("%.17f", 0.3)[1] "0.29999999999999999"(0.1 + 0.2) - 0.3[1] 5.551115e-17Most decimal fractions cannot be represented exactly with a finite number of binary digits. The computer stores nearby values; the difference here is tiny, but exact equality sees it.
When our question is whether two computed values are close enough, choose a tolerance appropriate to the scale of the problem:
tolerance <- 1e-12
abs((0.1 + 0.2) - 0.3) < tolerance[1] TRUEThis tolerance is suitable for this small arithmetic example; it is not a universal threshold for every model or unit. R also provides all.equal() for approximate comparisons.
Does the order of subtraction matter?
Try these two expressions before reading the output. In exact arithmetic, both are zero. R evaluates each one from left to right:
9.1 - 9 - 0.1[1] -3.608225e-169.1 - 0.1 - 9[1] 0On our R setup, the first leaves a tiny negative remainder (about \(-3.6\times10^{-16}\)), while the second gives exactly zero. The intermediate calculations round differently, even though the expressions are algebraically equivalent. Which error measure can we use when the correct answer is zero: absolute or relative error? Would either result matter if these were measurements in meters?
Can tiny errors accumulate?
Our simulations have used a vector of times made with seq(). What might happen if we repeatedly add \(0.1\) instead? Predict whether both methods below give exactly 10,000.
number_of_steps <- 100000
delta_t <- 0.1
time_by_addition <- 0
for (i in 1:number_of_steps) {
time_by_addition <- time_by_addition + delta_t
}
time_from_step_number <- number_of_steps * delta_t
time_by_addition[1] 10000time_from_step_number[1] 10000time_by_addition - time_from_step_number[1] 1.884837e-08The accumulated difference is small here. The general lesson is to measure the effect, not to assume every floating-point difference matters. When time is determined by the step number, calculating it from that number avoids repeated addition of the time step.
Meeting 2: Test the limits of a computation
Today we will read numbers the way Module 5.2 does, then make three predictions and test them in R: Can a number be too large or too small? Does arithmetic order matter? What changes when an Euler simulation uses a smaller time step?
First, recall yesterday’s distinction: a wrong unit or a missing biological mechanism will not be repaired by a more precise calculation. Today’s examples focus on what the computation itself can change.
Reading a number: notation, precision, and magnitude
The book uses normalized exponential notation with the decimal point before the first nonzero digit. This differs from the scientific notation many of us learned, which puts one nonzero digit before the decimal point. Both notations name the same value:
\[ 0.0004500 = 0.4500\times10^{-3} = 4.500\times10^{-4}. \]
The middle expression is normalized in the book’s convention. In \(0.4500\times10^{-3}\), the book calls 4500 the significand (the digits after dropping the decimal point), precision is 4 significant digits, and magnitude is \(10^{-3}\)—not merely the exponent \(-3\). The leading zeros do not count; the final two zeros after the decimal point do.
Try one together without R: write \(3{,}704{,}000\) in the book’s normalized notation. Which digits count as significant? What are its precision and magnitude? The book’s convention treats trailing zeros in an integer written without a decimal point as not significant, so the answer is \(0.3704\times10^7\), precision 4, magnitude \(10^7\).
These definitions help us discuss how many digits a stored or reported result carries. They do not tell us whether a model is scientifically accurate.
Finite range: overflow and underflow
A computer stores numbers within a finite range. Predict the two results:
1e308 * 10[1] Inf1e-300 / 1e100[1] 0In ordinary R numeric arithmetic, the first becomes Inf (too large) and the second becomes 0 (too small to represent here). Both should trigger questions if they appear unexpectedly in a simulation. In ESA’s account of Ariane 5 Flight 501, a value exceeded the capacity of a converted representation; that was not the same mechanism as 0.1 + 0.2.
Order of arithmetic can matter
Our \(9.1\) subtraction produced a tiny discrepancy. Mathematically, \((a+b)+c=a+(b+c)\) too—but finite precision can make the difference much larger. Predict these R results:
a <- 1e16
b <- -1e16
c <- 1
(a + b) + c[1] 1a + (b + c)[1] 0Why can the 1 be lost in one order? This is another consequence of finite precision. The textbook also discusses adding small quantities before large ones when feasible.
A finite-step simulation is also an approximation
Return to the growth model from Module 2.2: the starting population is 100, and its growth rate is 10% of its current size per hour. How many individuals does a one-hour Euler step add at the start? What would the first half-hour step add instead?
The first one-hour step adds \(1(0.1)(100)=10\), reaching 110. The first half-hour step adds \(0.5(0.1)(100)=5\), reaching 105. After that half hour, the rate must be calculated again from 105, not 100. Smaller steps let the rate change more often.
Code the familiar loop with a one-hour step. The exact value at hour 8 is available for this particular growth model, so we can measure the numerical error.
delta_t <- 1
times <- seq(0, 8, by = delta_t)
population <- numeric(length(times))
population[1] <- 100
for (i in 2:length(times)) {
old_population <- population[i - 1]
growth_rate <- 0.1 * old_population
population[i] <- old_population + delta_t * growth_rate
}
one_hour_result <- population[length(population)]
exact_result <- 100 * exp(0.1 * 8)
one_hour_result[1] 214.3589exact_result[1] 222.5541abs(one_hour_result - exact_result)[1] 8.195212Now change only the time step to half an hour. Predict whether the result will move closer to or farther from the exact value before running this chunk.
delta_t <- 0.5
times <- seq(0, 8, by = delta_t)
population <- numeric(length(times))
population[1] <- 100
for (i in 2:length(times)) {
old_population <- population[i - 1]
growth_rate <- 0.1 * old_population
population[i] <- old_population + delta_t * growth_rate
}
half_hour_result <- population[length(population)]
half_hour_result[1] 218.2875abs(half_hour_result - exact_result)[1] 4.266634The half-hour result is closer here. If time permits, change delta_t to 0.1 and then 0.01 in that second chunk. Does the answer appear to settle down? A smaller step can reduce time-step approximation error, but it cannot fix a wrong growth assumption or incorrect data. It also does not make floating-point roundoff disappear.
What would make you trust a result more?
For one of our growth, predator–prey, or SIR simulations, propose three checks. Consider units, conservation or other invariants, a different time step, and whether the assumptions fit the question. What would each check reveal—and what could it not reveal?
Before leaving: If the half-hour answer is closer to the exact solution, have we shown that the growth model accurately predicts a real population? Explain in one sentence.