Euler Method

Euler Method in Health Economic Models

The Euler method solves a differential equation one time step at a time. Starting from known values, it multiplies the current rate of change by the step length and adds the result to the current state, then repeats, so a continuous-time model becomes a row-by-row calculation. In health economics the method is a common engine inside spreadsheet versions of compartmental transmission models, system dynamics models and other continuous models, where each row of the sheet holds one step. Every number can be checked by hand. Its weaknesses are an error that shrinks only in proportion to the step and a stability limit that a long step can breach, and either can shift the costs and QALYs behind an incremental cost-effectiveness ratio.

The update rule and its derivation

An initial value problem gives the rate of change of a state together with its value at the start. In a cohort or population model the state is the number of people in each health state, and its rate of change is the sum of flows in and out. Euler's method, as set out in Süli's Oxford lecture notes, divides the time horizon into mesh points a fixed step apart and advances by the slope at the start of each step.

$$y_{n+1} = y_n + h\,f(t_n, y_n), \qquad t_n = t_0 + nh$$

where $y_n$ is the numerical approximation to the true state $y(t_n)$, $h$ is the step size, $f(t, y)$ is the rate of change given by the differential equation $dy/dt = f(t, y)$, $t_0$ is the start time and $y_0$ is the known initial value.

Süli derives the rule by integrating the differential equation across one step and replacing the integral of the slope with the rectangle rule, which uses the value at the left end of the interval for the whole interval. Euler's method is therefore the left rectangle rule of numerical integration applied to a slope that itself depends on the unknown state. The right end gives the implicit, or backward, Euler method.

$$y_{n+1} = y_n + h\,f(t_{n+1}, y_{n+1})$$

where the symbols are as before and the unknown $y_{n+1}$ appears on both sides, so an equation has to be solved at every step. Süli treats the explicit, or forward, Euler method as the one-stage member of the Runge-Kutta family, whose higher-order members gain accuracy by evaluating the slope at points inside the step.

How fast the error shrinks with the step

Süli measures accuracy through the truncation error, the amount by which the true solution fails to satisfy the Euler update per unit step. A Taylor series expansion shows that it equals half the step times the second derivative of the solution at a point inside the step, so the method is first-order accurate. For a smooth problem with an exact starting value, Süli bounds the global error left after many steps as follows.

$$|e_n| \le \frac{M_2\, h}{2L}\left(e^{L(t_n - t_0)} - 1\right)$$

where $e_n = y(t_n) - y_n$ is the global error at step $n$, $M_2$ is the largest absolute value of the second derivative of the solution over the interval and $L$ is the Lipschitz constant of $f$, a bound on how sharply the slope changes as the state changes.

The bound is linear in $h$, so halving a small step should roughly halve the error. The bound also grows exponentially with elapsed time, so a step that is adequate over a trial period may not be adequate over a lifetime time horizon.

Stability limits on the step

A long step can produce behaviour the true solution never shows. Süli analyses this with the test equation $y' = \lambda y$ for real negative $\lambda$, whose solution decays to zero, and shows that Euler's method then gives the following.

$$y_n = (1 + h\lambda)^{n}\, y_0$$

where $\lambda$ is the rate constant (negative for an outflow), $h$ is the step and $y_0$ is the starting value.

The numerical values decay only when $|1 + h\lambda| < 1$, so $h\lambda$ must lie between minus 2 and 0, the interval of absolute stability. For a health state that people leave at rate $r$ per unit time, $\lambda = -r$ and the step must satisfy $h < 2/r$. When $hr$ lies between 1 and 2 the values alternate in sign while shrinking, above 2 they alternate and grow, and any step with $hr$ above 1 moves more people out of the state than it contains.

With several linked states and constant rates the same condition applies to each eigenvalue of the rate matrix, and Süli shows that $h\lambda$ must lie inside a disc of radius 1 centred at minus 1 for every eigenvalue. Soetaert, Petzoldt and Setzer, describing the deSolve package for R, call a system stiff when it changes on a wide variety of time scales, and note that unless such problems are solved with specially designed methods they need very small time steps to satisfy stability requirements. A disease model with an acute phase measured in days inside a horizon measured in decades has this structure. Süli shows that the stability region of the implicit Euler method includes the whole left half of the complex plane (A-stability), at the price of solving an equation at each step.

Euler steps in transmission and system dynamics models

For dynamic transmission models written as differential equations, the ISPOR-SMDM task force report by Pitman and colleagues states that in a spreadsheet the Euler method is often used to project the system forward in time, although it is not as accurate as other numerical techniques. Crielaard and colleagues make the same point for system dynamics, observing that a difference equation implicitly uses Euler integration with a fixed time step and that the step size and solver can significantly affect the numerical error.

For an SIR model with frequency-dependent transmission, the Euler update computes every flow from the current state and then applies all of them together.

$$S_{n+1} = S_n - h\beta \frac{S_n I_n}{N}, \qquad I_{n+1} = I_n + h\beta \frac{S_n I_n}{N} - h\gamma I_n, \qquad R_{n+1} = R_n + h\gamma I_n$$

where $S_n$, $I_n$ and $R_n$ are the susceptible, infectious and recovered numbers at step $n$, $N$ is the population size, $\beta$ is the transmission rate and $\gamma$ is the recovery rate, both per unit time.

The changes sum to zero, so the total stays at $N$ apart from rounding. Using the already updated $S_{n+1}$ when computing the infectious flow mixes two time points and changes the method. A step for which $h\beta I_n / N$ exceeds 1 drives the susceptible number below zero, which signals a step too long for the force of infection, and clipping the value at zero changes the equations and breaks the population total.

Cohort cycles read as Euler steps

A Markov cohort model with a fixed cycle length has the same row-by-row form, and for a model built from rates the link is exact. One Euler step for a state left at a constant rate gives a simple multiplier.

$$y_{n+1} = (1 - rh)\, y_n$$

where $y_n$ is the number in the state at the start of cycle $n$, $r$ is the exit rate per unit time and $h$ is the cycle length.

Treating $rh$ as the per-cycle transition probability is therefore a forward Euler step, while the exact solution over one cycle gives $p = 1 - e^{-rh}$, the usual rate-to-probability conversion. For a rate of 0.4 per year and annual cycles the two probabilities are 0.400 and 0.330, and the gap narrows as the cycle shortens.

Costs, QALYs and discounting along the path

Model outcomes are integrals of the state over time, weighted by costs or utilities and discounted. Counting each step's outcomes from its starting state is the left rectangle rule and matches forward Euler. With continuous discounting the left-rectangle total cost of a state is as follows.

$$C = \sum_{n=0}^{K-1} h\, c\, y_n\, e^{-\rho t_n}$$

where $C$ is the discounted total cost, $K$ is the number of steps, $c$ is the cost per person per unit time in the state, $y_n$ is the number in the state at the start of step $n$, $\rho$ is the continuous discount rate and $t_n$ is the start time of step $n$.

Elbasha and Chhatwal defined seven methods from numerical integration for accumulating outcomes in discrete-time state-transition models and compared them with a continuous-time gold standard. They found that the standard half-cycle correction gave the same results as the trapezoidal rule, that errors need not cancel when incremental outcomes are calculated and that a less accurate method can lead to a wrong decision, with every method improving under a shorter cycle. The error in the state path and the error in the outcome rule are separate, so a model report needs to state the step used for the states and the rule used to accumulate QALYs and costs, including how discounting is timed.

Worked example with a declining cohort and a short hospital stay

All numbers in this example are illustrative. A cohort of 1,000 people starts in a treated disease state that people leave, through progression or death, at a constant rate of 0.4 per year. The exact number remaining after two years is 1,000 × $e^{-0.8}$ = 449.33.

1. Euler factor per step. Each step multiplies the state by $(1 - 0.4h)$, which is 0.6 for a one-year step, 0.8 for six months, 0.9 for three months and 0.95 for an eighth of a year.

2. State at two years. Raising each factor to the number of steps and multiplying by 1,000 gives the following.

3. Observed order. The error ratios are 89.33/39.73 = 2.25, 39.73/18.86 = 2.11 and 18.86/9.20 = 2.05, approaching 2 as a first-order method should. Without an exact solution, differences between successive runs (49.60, 20.87, 9.66) give ratios of 2.38 and 2.16.

4. Person-years and cost. The exact person-years in the state over two years are 1,000 × (1 − 0.44933)/0.4 = 1,376.68. With six-month steps, left-rectangle counting on the Euler path gives 0.5 × (1,000 + 800 + 640 + 512) = 1,476.00, an overstatement of 99.32, or 7.2%. The trapezoidal rule on the same path gives 0.5 × (500 + 800 + 640 + 512 + 204.80) = 1,328.40, an understatement of 48.28, because the Euler path lies below the true curve. At GBP 3,000 per person-year, the left-rectangle Euler total of GBP 4,428,000 overstates the exact GBP 4,130,033 by GBP 297,967.

5. A short-stay state. Suppose 100 people enter a hospital state with an exit rate of 52 per year, a mean stay of about one week. A monthly step gives a factor of 1 − 52/12 = −3.33, so the state holds −333.33 people after one month and 1,111.11 after two, against an exact 1.31 after one month. Stability needs a step below 2/52 of a year, about 14 days.

6. Interpretation. Six-month steps are too coarse for a cost total that should be right to within a few per cent, and the trapezoidal rule only partly offsets an error that comes from the state solver. A model with monthly cycles cannot run the short-stay state through forward Euler at all, so that part of the model needs a much shorter step or a solver built for stiff problems.

Choosing and checking the step

A step is adequate when halving it no longer changes the decision outputs by an amount that matters. Step checks form one part of model validation.

1. Match units. Step and rates share one time unit, so a yearly rate with a monthly step needs $h = 1/12$.
2. Update simultaneously. All flows for a step are computed from the same starting state before any state changes.
3. Respect the stability limit. The step satisfies $h < 2/r$ for the fastest exit rate, and $h \le 1/r$ if negative values are to be avoided.
4. Refine the step. Runs at $h$, $h/2$ and $h/4$ show errors falling by about half each time, for incremental costs, QALYs and net benefit as well as for the states.
5. Compare with another solver. The deSolve package for R, for example, offers the Euler method beside the classical fourth-order Runge-Kutta method and an lsoda routine that selects a stiff or non-stiff method automatically.

Euler's method beside neighbouring methods

The finite difference method is the wider family that replaces derivatives with differences, and forward Euler is its simplest time-stepping scheme. Numerical integration approximates the integral of a known function, whereas Euler's method advances a state whose slope depends on the state itself. The half-cycle correction changes how outcomes are counted within a cycle and leaves the state path untouched. Any numerical method still carries errors in the equations or parameters into its results.

Misreadings and limits

Three misreadings of Euler output each understate solver error.

1. Taking a smooth curve as an accurate one. The worked example shows a plausible decline that is about 9% too low at two years with six-month steps.
2. Expecting errors to cancel between arms. Elbasha and Chhatwal showed that integration errors need not cancel in incremental outcomes, so a solver error can move the ICER.
3. Ignoring solver error because parameter uncertainty is larger. Numerical bias is systematic, so it does not average out across probabilistic runs and a wide uncertainty interval does not absorb it.

Sources

• Crielaard L, Uleman JF, Châtel BDL, Epskamp S, Sloot PMA, Quax R. Refining the causal loop diagram: a tutorial for maximizing the contribution of domain expertise in computational system dynamics modeling. Psychological Methods. 2024;29(1):169-201. https://doi.org/10.1037/met0000484
• Elbasha EH, Chhatwal J. Theoretical foundations and practical applications of within-cycle correction methods. Medical Decision Making. 2016;36(1):115-131. https://doi.org/10.1177/0272989X15585121
• Pitman R, Fisman D, Zaric GS, Postma M, Kretzschmar M, Edmunds J, Brisson M. Dynamic transmission modeling: a report of the ISPOR-SMDM Modeling Good Research Practices Task Force-5. Value in Health. 2012;15(6):828-834. https://doi.org/10.1016/j.jval.2012.06.011
• Soetaert K, Petzoldt T, Setzer RW. Solving differential equations in R: package deSolve. Journal of Statistical Software. 2010;33(9):1-25. https://doi.org/10.18637/jss.v033.i09
• Süli E. Numerical Solution of Ordinary Differential Equations. Lecture notes, Mathematical Institute, University of Oxford; 2 March 2022. Sections 2.1, 2.2, 2.4, 2.5 and 4. https://people.maths.ox.ac.uk/suli/nsodes.pdf