3 Numerical Methods for Solving SIR Models

The SIR equations tell us the rate at which each compartment changes. But a rate of change is not the same as a trajectory. To get from “how fast things are moving” to “where things are at every point in time,” we need to solve the equations. This chapter introduces two ways to do that in code, starting with the simplest approach and building toward something more powerful.

3.1 Introduction to Solving SIR Models

Learning Objectives

By the end of Section 3.1, you should be able to:

Explain what “solving” a differential equation (or a system of differential equations) means and interpret it in the context of solving an SIR model.

When we talk about “solving” a differential equation (or system of differential equations like our SIR model), we are essentially answering the question:

What function changes at this rate?

In Calculus 1, we learned that to find the rate of change of a function we take its derivative. Solving a differential equation means undoing this process (i.e. finding the function whose derivative is equal to our differential equation). To fully solve a differential equation we also need some initial values, or values of the function at some specific times.

Why do we need initial values to solve a differential equation?

Remembering back to our derivative rules, we know that the derivative of any constant is always $0$. So when we are given a differential equation (which is just a derivative of some function) and are looking to “go back” to the function that produced it, we can only obtain a solution down to the precision of some unknown constant. Our initial values help us solve for that unknown constant and “narrow down” the potential solutions down to one specific one.

To illustrate this, let’s say you are given the differential equation $\frac{dy}{dx}=3x^2$. We know that the function $f(x)=y=x^3$ has a derivative equal to our differential equation (and thus is a solution to it). But the functions $f(x)=x^3+1$, $f(x)=x^3+2$,$f(x)=x^3+1003$, and so on also have derivatives equation to our differential equation and thus are solutions as well. Hence, given no other information we can only say the solution to the differential equation $\frac{dy}{dx}=3x^2$ is $f(x)=x^3+C$ where $C$ is an unknown constant. If we were given an initial value, however, we could solve for this unknown constant and know our solution exactly. Say we have the initial value $f(1)=2$, then by solving for $C$ in $f(x)=x^3+C$ we know that $C=1$ and our exact solution must be $f(x)=x^3+1$.

When we say we are “solving an SIR model” (or “solving SIR equations”), we are answering the question:

Given that we know how fast the susceptible, infected, and recovered populations are changing and we know the initial number / proportion of people in each compartment at the beginning of the epidemic, what function gives us the number / proportion of people in each compartment at various times?

As you can imagine, solving our SIR model is extremely useful since we are most interested in knowing the number / proportion of susceptible, infectious, and recovered people at various times (not their rates of change).

Two common methods of solving differential equations that are difficult to solve by hand is:

Euler’s method
Fourth-Order Runge-Kutta method (implemented using the pracma package in R)

We will go over each method separately in Section 3.2 and Section 3.3 and the differences between the two in Section 3.4.

3.2 Euler’s Method

Learning Objectives

By the end of Section 3.2, you should be able to:

Explain the step-by-step mathematical process for solving differential equation (or a system of differential equations) using Euler’s method.

Euler’s method estimates the solution to a differential equation using a series of short, connected linear approximations at a discrete set of times. Although Euler’s method is widely regarded to be a relatively inefficient method for solving otherwise “unsolvable” differential equations, it is still valued for its simplicity and accessibility.

Suppose we have approximate initial values for the number of susceptible, infectious, and recovered individuals at the beginning of our study period. We will denote these values as $S(0)$, $I(0)$, and $R(0)$, respectively.

Note

To keep it simple, we will only discuss the process of Euler’s method as it is applied to the number of individuals in each compartment, but keep in mind this process can be applied in exactly the same way to proportions as well.

Our goal is to find approximate solutions to $S(t)$, $I(t)$, and $R(t)$—the exact number of individuals in each compartment at time $t$—using our initial values ($S(0)$, $I(0)$, and $R(0)$) and system of differential equations that is our SIR model (which we found in Section 2.1). For clarity, we will be using $S(t)$, $I(t)$, and $R(t)$ to denote the exact solutions to the SIR equations and $s(t)$, $i(t)$, and $r(t)$ to denote the approximate solutions given by Euler’s method.

Note

These types of problems, one where we are given a set of intial values and corresponding differential equations, are known as initial-value problems (or IVPs).

Euler’s method involves repeatedly finding tangent line approximations to $S(t)$, $I(t)$, and $R(t)$ at a discrete set of times (often called nodes). We will use $s_t$, $i_t$, and $r_t$ to denote these approximations at a time $t$.

What is a tangent line approximation?

A tangent line approximation uses the equation of the tangent line at a point on a function to estimate function values at nearby points.

Given a function $f(x)$ who has a derivative $f'(x)$, we can find the equation of a tangent line (denoted with $y$) at the point $(a,f(a))$ using the formula:

\[ y=f(a)+f'(a)(x-a) \]

This formula might look familiar because it is just the point slope formula from Algebra 1 but with the inclusion of the derivative $f'(a)$ as the slope.

Now that we have this equation for the tangent line, we can use it to approximate function values at nearby points. Say we want to estimate the value of $a+\Delta t$ where $\Delta t$ is just some change on the x-axis to the right or left of $a$. We can do so by plugging in $a+\Delta t$ for $x$ into our tangent line equation:

\[ y=f(a)+f'(a)(\Delta t) \]

Whatever we get out for $y$ is our tangent line approximation for $a+\Delta t$ using the tangent line at $x=a$.

These approximations are found using the following recursive formulas:

\[ \begin{aligned} s_{t+1} &= s_t - \Delta t \ (\beta) \left( \frac{i_t}{N} \right) (s_t) \\[10pt] i_{t+1} &= i_t + \Delta t \left[\beta \left( \frac{i_t}{N} \right) (s_t) - \gamma (i_t) \right] \\[10pt] r_{t+1} &= r_t + \Delta t \ (\gamma) (i_t) \end{aligned} \]

$s_{t+1}$, $i_{t+1}$, $r_{t+1}$: Approximations for $S(t)$, $I(t)$, and $R(t)$ at the next node
$s_{t}$, $i_{t}$, $r_{t}$: Approximations for $S(t)$, $I(t)$, and $R(t)$ at the previous node
$\Delta t$: Step size (how far each node is spaced apart on the $t$-axis)
$N$, $\beta$, $\gamma$: Population size, transmission rate, and recovery rate, respectively

How does this connect back to the “general” way of doing Euler’s method?

This formula is just an application of the generalized Euler’s method for first-order differential equation $Y'(t)=f\big(t,Y(t)\big)$, true solution $Y(t)$, and approximate solutions $y_0,\dots,y_t,\dots,y_b$, which are calculated at evenly spaced $t_n$ as:

\[ y_{t+1}=y_t+ \Delta t \big[f\big(t_n,Y(t_n)\big)\big] \]

Where $\Delta t$ still represents how far each $t_n$ are apart.

This process can be a little confusing to grasp at first, so let’s do an example. Suppose the following SIR model equations with initial values $S(0)=99$, $I(0)=1$, $E(0)=0$, population size $N=100$, transmission rate $\beta=1.5$, and recovery rate $\gamma=0.5$:

\[ \begin{aligned} \frac{dS}{dt} &= -1.5 \tilde{I} S \\[10pt] \frac{dI}{dt} &= 1.5 \tilde{I} S - 0.5 I \\[10pt] \frac{dR}{dt} &= 0.5 I \end{aligned} \]

Reminder

Remember that $\tilde{I}=\frac{I}{N}$.

The following table illustrates the first few steps of the calculation process to solve this SIR model using Euler’s method with a step size of $\Delta t=1$:

Time	Approximations	Calculation Process
$t=0$	$\textcolor{blue}{s(0)=99.0}$ $\textcolor{red}{i(0)=1.0}$ $\textcolor{green}{r(0)=0.0}$	—
$t=1$	$\textcolor{blue}{s(1)=97.5}$ $\textcolor{red}{i(1)=2.0}$ $\textcolor{green}{r(1)=0.5}$	$s(1)=\textcolor{blue}{99.0}-(1.0)(1.5)\left(\frac{\textcolor{red}{1.0}}{100.0}\right)(\textcolor{blue}{99.0})$ $i(1)=\textcolor{red}{1.0}+(1.0)\left[(1.5)\left(\frac{\textcolor{red}{1.0}}{100.0}\right)(\textcolor{blue}{99.0})-0.5(\textcolor{red}{1.0})\right]$ $r(1)=\textcolor{green}{0.0}+(1.0)(0.5)(\textcolor{red}{1.0})$
$t=2$	$\textcolor{blue}{s(2)=94.6}$ $\textcolor{red}{i(2)=3.9}$ $\textcolor{green}{r(2)=2.0}$	$s(2)=\textcolor{blue}{97.5}-(1.0)(1.5)\left(\frac{\textcolor{red}{2.0}}{100.0}\right)(\textcolor{blue}{97.5})$ $i(2)=\textcolor{red}{2.0}+(1.0)\left[(1.5)\left(\frac{\textcolor{red}{2.0}}{100.0}\right)(\textcolor{blue}{97.5})-0.5(\textcolor{red}{2.0})\right]$ $r(2)=\textcolor{green}{0.5}+(1.0)(0.5)(\textcolor{red}{2.0})$
$t=3$	$\textcolor{blue}{s(3)=89.1}$ $\textcolor{red}{i(3)=7.5}$ $\textcolor{green}{r(3)=3.4}$	$s(3)=\textcolor{blue}{94.6}-(1.0)(1.5)\left(\frac{\textcolor{red}{3.9}}{100.0}\right)(\textcolor{blue}{94.6})$ $i(3)=\textcolor{red}{3.9}+(1.0)\left[(1.5)\left(\frac{\textcolor{red}{3.9}}{100.0}\right)(\textcolor{blue}{94.6})-0.5(\textcolor{red}{3.9})\right]$ $r(3)=\textcolor{green}{2.0}+(1.0)(0.5)(\textcolor{red}{3.9})$
$t=4$	$\textcolor{blue}{s(4)=79.1}$ $\textcolor{red}{i(4)=13.7}$ $\textcolor{green}{r(4)=7.2}$	$s(4)=\textcolor{blue}{89.1}-(1.0)(1.5)\left(\frac{\textcolor{red}{7.5}}{100.0}\right)(\textcolor{blue}{89.1})$ $i(4)=\textcolor{red}{7.5}+(1.0)\left[(1.5)\left(\frac{\textcolor{red}{7.5}}{100.0}\right)(\textcolor{blue}{89.1})-0.5(\textcolor{red}{7.5})\right]$ $r(4)=\textcolor{green}{3.4}+(1.0)(0.5)(\textcolor{red}{7.5})$
$t=5$	$\textcolor{blue}{s(5)=62.8}$ $\textcolor{red}{i(5)=23.1}$ $\textcolor{green}{r(5)=14.0}$	$s(5)=\textcolor{blue}{79.1}-(1.0)(1.5)\left(\frac{\textcolor{red}{13.7}}{100.0}\right)(\textcolor{blue}{79.1})$ $i(5)=\textcolor{red}{13.7}+(1.0)\left[(1.5)\left(\frac{\textcolor{red}{13.7}}{100.0}\right)(\textcolor{blue}{79.1})-0.5(\textcolor{red}{13.7})\right]$ $r(5)=\textcolor{green}{7.2}+(1.0)(0.5)(\textcolor{red}{13.7})$

What does $\approx$ mean?

The symbol $\approx$ means “approximately equal to” and indicates that the number has been rounded.

Important Note

The accuracy of the solution produced by Euler’s method is a direct result of what we choose the step size to be. If the step size is larger, then the curvature of the true solution will be lost since we are repeatedly taking linear approximations, resulting in an approximate solution that deviates a lot from the true solution curve. On the other hand, as the step size gets smaller, the approximate solution provided by Euler’s method converges to the true solution curve.

Read on to Section 3.2.1 to experiment with what happens when you change the step size.

You’re probably noticing by now that the calculation process above is very tedious, which is why Euler’s method is most often calculated using code rather than by hand. Read on to Section 3.2.1 for a demonstration on how to code Euler’s method.

3.2.1 Euler’s Method Code Demonstration

Learning Objectives

By the end of Section 3.2.1, you should be able to:

Translate the mathematical process of Euler’s method into code.

This section contains a step by step walkthrough of how to code an Euler’s method algorithm in R.

Step 1: Define the biological parameters.

Every model starts with its parameters: the biological quantities that characterize the disease and host system. Here, we have two: how fast the disease spreads ($\beta$) and how fast individuals recover ($\gamma$).

Step 2: Set the initial conditions.

Before the model can run, we need to specify where the system starts. In an epidemic model, this means the proportion of the population in each compartment at time zero.

This sanity check, confirming that proportions sum to 1, is a simple example of a habit that will save you debugging time later. When you add more compartments (SEIR, endemic model), this check generalizes naturally.

Step 3: Set up the time grid.

Euler’s method steps through time discretely. We need to decide how far apart those steps are ($\Delta t$, the step size) and how long to simulate.

The step size is a computational choice, not a biological one, but it affects accuracy. A smaller $\Delta t$ gives a more accurate approximation of the continuous solution, at the cost of more computation. We will explore this tradeoff shortly.

Step 4: Pre-allocate storage.

Before running the loop, we create empty vectors to hold the results at each time step. This is more efficient than growing a vector inside the loop.

Think of these vectors as the model’s memory: they store the entire trajectory of the epidemic, one row per time step.

Step 5: Run the Euler loop.

This is where the model actually runs. At each time step, we calculate the rates of change using the SIR equations and use them to update the compartment values. Each line in the loop maps directly to an equation.

The loop is doing exactly what the Euler update formulas say, mechanically, one step at a time. This is the computational equivalent of doing the hand calculations from the table shown earlier, just thousands of times faster.

The loop is the model

Notice that the three lines computing dS, dI, and dR are not arbitrary calculations. They are the SIR equations, written in R. Every biological assumption we made in the previous chapter is present in those three lines. The loop is not running the model alongside the equations, the loop is the equations.

Step 6: Visualize the results.

The vectors S, I, and R now hold the full epidemic trajectory. We can plot them against time.

Try this

Go back to Step 3 and try the following changes, re-running each step in order after each change:

Change dt <- 0.1 to dt <- 2.0. What happens to the curves? This is the cost of a large step size.
Change the step size to dt <- 0.001. Did the code output results faster or slower than when using dt <- 2.0? (NOTE: How fast the code runs depends on the computer you are using right now, but in general note that the smaller our step size gets, the more accurate our results become at the cost of being increasingly computationally demanding.)

3.3 Solving Models Using the `pracma` Package in R

Learning Objectives

By the end of Section 3.3, you should be able to:

Explain what values the function rk4sys in pracma takes in.
Translate an SIR model to a function and solve it using rk4sys.

Euler’s method is honest: every step is visible, every assumption is explicit. Its weakness is that those small linear steps accumulate error, especially when the epidemic curve is changing quickly. The next method, RK4, is more sophisticated. It looks further ahead before committing to each step, which makes it more accurate. The tradeoff is that it is less transparent.

The package pracma in R, short for “Practical Numerical Math Functions,” provides a large repertoire of functions used in numerical analysis, linear algebra, numerical optimization, and differential equations (Borchers 2025).

Is there a specific name for the type of differential equations we are solving?

The SIR models this module covers deal with systems of first-order ordinary differential equations.

Breaking up each word in this description, we have:

First-order: This just means we are only dealing with equations involving the first derivative and not subsequent derivatives (the second, third, fourth, etc. derivative)
Ordinary: The unknown function which we are trying to solve for depends only one independent variable.

Although pracma contains many mathematical functions, we are primarily interested in rk4sys, a function to solve systems of differential equations using an advanced method known as the Classic Fourth-Order Runge-Kutta Method (or RK4 for short). This method for solving differential equations is much more complex than Euler’s method but outputs solutions to a much greater accuracy. Like Euler’s method, RK4’s accuracy depends on the step size. However, in general, a solution to a model found using RK4 will be more accurate than one found using Euler’s method, even when the same step size is used.

How exactly does RK4 work?

In simple terms, RK4 estimates the value of a solution at the next point by taking a weighted average of four slope calculations, each computed at different points within the step interval. The Wikipedia article on Runge-Kutta methods goes more into the specific mathematical formulation of RK4 (Wikipedia contributors 2026).

The following table summarizes the values rk4sys takes in:

Variable Name	What Does It Take In?
`f`	A user-defined function that defines the system of differential equations to be solved.
`a`	A number specifying the starting point of the interval we want solutions for.
`b`	A number specifying the endpoint of the interval we want solutions for.
`y0`	A vector specifying the initial values of the system. A system of $m$ differential equations must have a `y0` of length $m$.
`n`	A number specifying the step size (i.e. the number of steps that should be taken from `a` to `b`.

The function rk4sys returns a matrix. When that matrix is converted to a data frame using the function as.data.frame, its first column provides the time points at which the SIR model was solved and subsequent columns are the solutions to the model’s state variables at those times. The state variables for an SIR model are the number or proportion of individuals susceptible, infectious, and recovered at a given time.

Before we can call rk4sys, let’s go over how to code the function to be passed into the parameter f. The function below is specifying a classic SIR model on the proportion of individuals in each compartment (see Section 2.1).

Step 1: Code the model function.

Step 2: Specify the parameters needed to call rk4sys.

We now have everything we need to call rk4sys.

Step 3: Call rk4sys using our defined values.

Now that we have a data frame of the solved SIR model equations at specified times, let’s plot it.

Step 4: Plot the results.

Try this

Go back to Step 2 and try the following changes, re-running each step in order after each change:

Change step_size <- 10 to step_size <- 50. What changes about the final plot? Is this similar to what happened when we changed the step size in our code for Euler’s method (see Section 3.2.1)?

3.4 Evaluating and Comparing Solving Methods

Learning Objectives

By the end of Section 3.4, you should be able to:

Explain the benefits and drawbacks of Euler’s method and rk4sys in relation to solving SIR models.
Identify which solving method is more suitable for a particular SIR model or application.
Visualize the difference between solutions produced using Euler’s method and rk4sys as step size changes.

The following table summarizes some potential benefits and drawbacks to using Euler’s method or rk4sys in pracma:

Solving Method	Potential Benefits	Potential Drawbacks
Euler’s Method	• Simple and intuitive code implementation • Easy to understand and visualize step-by-step • No external packages required • Can be implemented in Excel or by hand	• Low accuracy unless a very small step size is used • Computationally inefficient for small step sizes or complex models • Can be tedious or error-prone to calculate by hand for many steps • Can accumulate large numerical errors quickly over long periods
`rk4sys`	• Uses a sophisticated solving method (Fourth-Order Runge-Kutta) that has a high accuracy • Handles large models with many system equations well • Preferred for complex models / advanced applications • Accuracy scales better with step size compared to Euler’s method	• Less transparent: not much understanding of what is going on “under the hood” • Requires the external package `pracma` • Output accuracy depends on step size • Each step is more computationally expensive than Euler’s method since four derivative approximations are required

Tip

Although the accuracy of the solutions output by Euler’s method and the Fourth-Order Runge-Kutta method rely on step size, Runge-Kutta typically outputs more accurate results than Euler’s method even with the same step size.

This means we can “get away” with using a smaller step size when using rk4sys over Euler’s method (thus saving us some computational cost / time).

The following simulation shows the difference between the solution curves produced by Euler’s method and rk4sys for different values of the model parameters and step size:

Try this

Try changing the step size from $0.25$ to $0.1$. Do the curves solved using Euler’s method and Runge-Kutta method become closer or further apart? Can you explain why mathematically?
Try changing the step size from $0.5$ to $0.1$. Which set of curves (the ones solved using Euler’s method or Runge-Kutta method) change more? Can you explain why mathematically (think about which solving method scales better with step size in terms of accuracy)?

This isn’t the only way!

There are many packages in R that can solve SIR equations, and the package deSolve is a really common choice among them. deSolve uses very advanced and versatile solving methods that work well with large or complex SIR models. See the CRAN documentation for more details (Soetaert and Petzoldt 2026).

We now have two tools for solving SIR models. The next question is: how do we know which values of $\beta$ and $\gamma$ to use? That is the problem of fitting a model to data, and it is the topic of the next chapter.

Chapter References

Borchers, Hans W. 2025. Pracma: Practical Numerical Math Functions. https://doi.org/10.32614/CRAN.package.pracma.

Soetaert, Karline, and Thomas Petzoldt. 2026. deSolve: Solvers for Initial Value Problems of Differential Equations. https://doi.org/https://doi.org/10.32614/CRAN.package.deSolve.

Wikipedia contributors. 2026. “Runge–Kutta Methods — Wikipedia, the Free Encyclopedia.” https://en.wikipedia.org/w/index.php?title=Runge%E2%80%93Kutta_methods&oldid=1339475672.

--- title: "Numerical Methods for Solving SIR Models" engine: knitr webr: packages: - pracma --- The SIR equations tell us the rate at which each compartment changes. But a rate of change is not the same as a trajectory. To get from "how fast things are moving" to "where things are at every point in time," we need to solve the equations. This chapter introduces two ways to do that in code, starting with the simplest approach and building toward something more powerful. ## Introduction to Solving SIR Models {#sec-solving-sir-equations} ::: {.callout-tip title="Learning Objectives"} By the end of @sec-solving-sir-equations, you should be able to: - Explain what "solving" a differential equation (or a system of differential equations) means and interpret it in the context of solving an SIR model. ::: When we talk about "solving" a differential equation (or system of differential equations like our SIR model), we are essentially answering the question: > What function changes at this rate? In Calculus 1, we learned that to find the rate of change of a function we take its *derivative*. Solving a differential equation means undoing this process (i.e. finding the function whose derivative is equal to our differential equation). To fully solve a differential equation we also need some *initial values*, or values of the function at some specific times. :::{.callout-note collapse="true"} ## Why do we need initial values to solve a differential equation? Remembering back to our derivative rules, we know that the derivative of any constant is always $0$. So when we are given a differential equation (which is just a derivative of some function) and are looking to "go back" to the function that produced it, we can only obtain a solution down to the precision of some unknown constant. Our initial values help us solve for that unknown constant and "narrow down" the potential solutions down to one specific one. To illustrate this, let's say you are given the differential equation $\frac{dy}{dx}=3x^2$. We know that the function $f(x)=y=x^3$ has a derivative equal to our differential equation (and thus is a solution to it). But the functions $f(x)=x^3+1$, $f(x)=x^3+2$,$f(x)=x^3+1003$, and so on also have derivatives equation to our differential equation and thus are solutions as well. Hence, given no other information we can only say the solution to the differential equation $\frac{dy}{dx}=3x^2$ is $f(x)=x^3+C$ where $C$ is an unknown constant. If we were given an initial value, however, we could solve for this unknown constant and know our solution exactly. Say we have the initial value $f(1)=2$, then by solving for $C$ in $f(x)=x^3+C$ we know that $C=1$ and our exact solution must be $f(x)=x^3+1$. ::: When we say we are "solving an SIR model" (or "solving SIR equations"), we are answering the question: > Given that we know how fast the susceptible, infected, and recovered populations are changing and we know the initial number / proportion of people in each compartment at the beginning of the epidemic, what function gives us the number / proportion of people in each compartment at various times? As you can imagine, solving our SIR model is extremely useful since we are most interested in knowing the number / proportion of susceptible, infectious, and recovered people at various times (not their rates of change). Two common methods of solving differential equations that are difficult to solve by hand is: 1. **Euler's method** 2. **Fourth-Order Runge-Kutta method** (implemented using the `pracma` package in R) We will go over each method separately in [@sec-euler-method] and [@sec-pracma] and the differences between the two in [Section @sec-evaluating-and-comparing-both-methods]. ## Euler's Method {#sec-euler-method} ::: {.callout-tip title="Learning Objectives"} By the end of @sec-euler-method, you should be able to: - Explain the step-by-step mathematical process for solving differential equation (or a system of differential equations) using Euler's method. ::: Euler's method estimates the solution to a differential equation using a series of short, connected linear approximations at a discrete set of times. Although Euler's method is widely regarded to be a relatively inefficient method for solving otherwise "unsolvable" differential equations, it is still valued for its simplicity and accessibility. Suppose we have approximate initial values for the number of susceptible, infectious, and recovered individuals at the beginning of our study period. We will denote these values as $S(0)$, $I(0)$, and $R(0)$, respectively. ::: {.callout-note} To keep it simple, we will only discuss the process of Euler's method as it is applied to the number of individuals in each compartment, but keep in mind this process can be applied in exactly the same way to proportions as well. ::: Our goal is to find approximate solutions to $S(t)$, $I(t)$, and $R(t)$---the exact number of individuals in each compartment at time $t$---using our initial values ($S(0)$, $I(0)$, and $R(0)$) and system of differential equations that is our SIR model (which we found in [Section @sec-the-classic-sir-epidemic-model]). For clarity, we will be using $S(t)$, $I(t)$, and $R(t)$ to denote the **exact** solutions to the SIR equations and $s(t)$, $i(t)$, and $r(t)$ to denote the *approximate* solutions given by Euler's method. ::: {.callout-note} These types of problems, one where we are given a set of intial values and corresponding differential equations, are known as initial-value problems (or IVPs). ::: Euler's method involves repeatedly finding *tangent line approximations* to $S(t)$, $I(t)$, and $R(t)$ at a discrete set of times (often called nodes). We will use $s_t$, $i_t$, and $r_t$ to denote these approximations at a time $t$. :::{.callout-note collapse="true"} ## What is a tangent line approximation? A tangent line approximation uses the equation of the tangent line at a point on a function to estimate function values at nearby points. Given a function $f(x)$ who has a derivative $f'(x)$, we can find the equation of a tangent line (denoted with $y$) at the point $(a,f(a))$ using the formula: $$ y=f(a)+f'(a)(x-a) $$ This formula might look familiar because it is just the *point slope formula* from Algebra 1 but with the inclusion of the derivative $f'(a)$ as the slope. Now that we have this equation for the tangent line, we can use it to approximate function values at nearby points. Say we want to estimate the value of $a+\Delta t$ where $\Delta t$ is just some change on the x-axis to the right or left of $a$. We can do so by plugging in $a+\Delta t$ for $x$ into our tangent line equation: $$ y=f(a)+f'(a)(\Delta t) $$ Whatever we get out for $y$ is our tangent line approximation for $a+\Delta t$ using the tangent line at $x=a$. ::: These approximations are found using the following recursive formulas: $$ \begin{aligned} s_{t+1} &= s_t - \Delta t \ (\beta) \left( \frac{i_t}{N} \right) (s_t) \\[10pt] i_{t+1} &= i_t + \Delta t \left[\beta \left( \frac{i_t}{N} \right) (s_t) - \gamma (i_t) \right] \\[10pt] r_{t+1} &= r_t + \Delta t \ (\gamma) (i_t) \end{aligned} $$ - $s_{t+1}$, $i_{t+1}$, $r_{t+1}$: Approximations for $S(t)$, $I(t)$, and $R(t)$ at the next node - $s_{t}$, $i_{t}$, $r_{t}$: Approximations for $S(t)$, $I(t)$, and $R(t)$ at the previous node - $\Delta t$: Step size (how far each node is spaced apart on the $t$-axis) - $N$, $\beta$, $\gamma$: Population size, transmission rate, and recovery rate, respectively ::: {.callout-note collapse="true"} ## How does this connect back to the "general" way of doing Euler's method? This formula is just an application of the generalized Euler's method for first-order differential equation $Y'(t)=f\big(t,Y(t)\big)$, true solution $Y(t)$, and approximate solutions $y_0,\dots,y_t,\dots,y_b$, which are calculated at evenly spaced $t_n$ as: $$ y_{t+1}=y_t+ \Delta t \big[f\big(t_n,Y(t_n)\big)\big] $$ Where $\Delta t$ still represents how far each $t_n$ are apart. ::: This process can be a little confusing to grasp at first, so let's do an example. Suppose the following SIR model equations with initial values $S(0)=99$, $I(0)=1$, $E(0)=0$, population size $N=100$, transmission rate $\beta=1.5$, and recovery rate $\gamma=0.5$: $$ \begin{aligned} \frac{dS}{dt} &= -1.5 \tilde{I} S \\[10pt] \frac{dI}{dt} &= 1.5 \tilde{I} S - 0.5 I \\[10pt] \frac{dR}{dt} &= 0.5 I \end{aligned} $$ :::{.callout-tip title="Reminder"} Remember that $\tilde{I}=\frac{I}{N}$. ::: The following table illustrates the first few steps of the calculation process to solve this SIR model using Euler's method with a step size of $\Delta t=1$: | **Time** | **Approximations** | **Calculation Process** | |-----------------|---------------------|----------------------------------| | $t=0$ | $\textcolor{blue}{s(0)=99.0}$ $\textcolor{red}{i(0)=1.0}$ $\textcolor{green}{r(0)=0.0}$ | — | | $t=1$ | $\textcolor{blue}{s(1)=97.5}$ $\textcolor{red}{i(1)=2.0}$ $\textcolor{green}{r(1)=0.5}$ | $s(1)=\textcolor{blue}{99.0}-(1.0)(1.5)\left(\frac{\textcolor{red}{1.0}}{100.0}\right)(\textcolor{blue}{99.0})$ $i(1)=\textcolor{red}{1.0}+(1.0)\left[(1.5)\left(\frac{\textcolor{red}{1.0}}{100.0}\right)(\textcolor{blue}{99.0})-0.5(\textcolor{red}{1.0})\right]$ $r(1)=\textcolor{green}{0.0}+(1.0)(0.5)(\textcolor{red}{1.0})$ | | $t=2$ | $\textcolor{blue}{s(2)=94.6}$ $\textcolor{red}{i(2)=3.9}$ $\textcolor{green}{r(2)=2.0}$ | $s(2)=\textcolor{blue}{97.5}-(1.0)(1.5)\left(\frac{\textcolor{red}{2.0}}{100.0}\right)(\textcolor{blue}{97.5})$ $i(2)=\textcolor{red}{2.0}+(1.0)\left[(1.5)\left(\frac{\textcolor{red}{2.0}}{100.0}\right)(\textcolor{blue}{97.5})-0.5(\textcolor{red}{2.0})\right]$ $r(2)=\textcolor{green}{0.5}+(1.0)(0.5)(\textcolor{red}{2.0})$ | | $t=3$ | $\textcolor{blue}{s(3)=89.1}$ $\textcolor{red}{i(3)=7.5}$ $\textcolor{green}{r(3)=3.4}$ | $s(3)=\textcolor{blue}{94.6}-(1.0)(1.5)\left(\frac{\textcolor{red}{3.9}}{100.0}\right)(\textcolor{blue}{94.6})$ $i(3)=\textcolor{red}{3.9}+(1.0)\left[(1.5)\left(\frac{\textcolor{red}{3.9}}{100.0}\right)(\textcolor{blue}{94.6})-0.5(\textcolor{red}{3.9})\right]$ $r(3)=\textcolor{green}{2.0}+(1.0)(0.5)(\textcolor{red}{3.9})$ | | $t=4$ | $\textcolor{blue}{s(4)=79.1}$ $\textcolor{red}{i(4)=13.7}$ $\textcolor{green}{r(4)=7.2}$ | $s(4)=\textcolor{blue}{89.1}-(1.0)(1.5)\left(\frac{\textcolor{red}{7.5}}{100.0}\right)(\textcolor{blue}{89.1})$ $i(4)=\textcolor{red}{7.5}+(1.0)\left[(1.5)\left(\frac{\textcolor{red}{7.5}}{100.0}\right)(\textcolor{blue}{89.1})-0.5(\textcolor{red}{7.5})\right]$ $r(4)=\textcolor{green}{3.4}+(1.0)(0.5)(\textcolor{red}{7.5})$ | | $t=5$ | $\textcolor{blue}{s(5)=62.8}$ $\textcolor{red}{i(5)=23.1}$ $\textcolor{green}{r(5)=14.0}$ | $s(5)=\textcolor{blue}{79.1}-(1.0)(1.5)\left(\frac{\textcolor{red}{13.7}}{100.0}\right)(\textcolor{blue}{79.1})$ $i(5)=\textcolor{red}{13.7}+(1.0)\left[(1.5)\left(\frac{\textcolor{red}{13.7}}{100.0}\right)(\textcolor{blue}{79.1})-0.5(\textcolor{red}{13.7})\right]$ $r(5)=\textcolor{green}{7.2}+(1.0)(0.5)(\textcolor{red}{13.7})$ | :::{.callout-note collapse="true"} ## What does $\approx$ mean? The symbol $\approx$ means "approximately equal to" and indicates that the number has been rounded. ::: ::: {.callout-important title="Important Note"} The accuracy of the solution produced by Euler's method is a direct result of what we choose the step size to be. If the step size is larger, then the curvature of the true solution will be lost since we are repeatedly taking linear approximations, resulting in an approximate solution that deviates a lot from the true solution curve. On the other hand, as the step size gets smaller, the approximate solution provided by Euler's method converges to the true solution curve. Read on to [Section @sec-euler-code] to experiment with what happens when you change the step size. ::: You're probably noticing by now that the calculation process above is very tedious, which is why Euler's method is most often calculated using code rather than by hand. Read on to [Section @sec-euler-code] for a demonstration on how to code Euler's method. ### Euler's Method Code Demonstration {#sec-euler-code} ::: {.callout-tip title="Learning Objectives"} By the end of @sec-euler-code, you should be able to: - Translate the mathematical process of Euler's method into code. ::: This section contains a step by step walkthrough of how to code an Euler's method algorithm in R. **Step 1: Define the biological parameters.** Every model starts with its parameters: the biological quantities that characterize the disease and host system. Here, we have two: how fast the disease spreads ($\beta$) and how fast individuals recover ($\gamma$). ```{webr-r} # Transmission rate: contacts × probability of transmission per unit time beta <- 1.5 # Recovery rate: 1/gamma gives the average duration of infection gamma <- 0.5 # With these values, the average infectious period is: cat("Average infectious period:", 1 / gamma, "time units") ``` --- **Step 2: Set the initial conditions.** Before the model can run, we need to specify where the system starts. In an epidemic model, this means the proportion of the population in each compartment at time zero. ```{webr-r} # Starting proportions — must sum to 1 S0 <- 0.99 # 99% susceptible: nearly everyone is susceptible at the start I0 <- 0.01 # 1% infectious: small "spark" of infection enters population R0_start <- 0.00 # 0% recovered: no prior immunity # Verify they sum to 1 (a good sanity check to build into your code) cat("S0 + I0 + R0 =", S0 + I0 + R0_start) ``` This sanity check, confirming that proportions sum to 1, is a simple example of a habit that will save you debugging time later. When you add more compartments (SEIR, endemic model), this check generalizes naturally. --- **Step 3: Set up the time grid.** Euler's method steps through time discretely. We need to decide how far apart those steps are ($\Delta t$, the step size) and how long to simulate. ```{webr-r} dt <- 0.1 # step size — try making this larger later to see what breaks t_end <- 50 # total simulation time # Create a vector of all time points times <- seq(0, t_end, by = dt) n_steps <- length(times) cat("Number of time steps:", n_steps, "\n") cat("Time points run from", times[1], "to", times[n_steps]) ``` The step size is a computational choice, not a biological one, but it affects accuracy. A smaller $\Delta t$ gives a more accurate approximation of the continuous solution, at the cost of more computation. We will explore this tradeoff shortly. --- **Step 4: Pre-allocate storage.** Before running the loop, we create empty vectors to hold the results at each time step. This is more efficient than growing a vector inside the loop. ```{webr-r} # Empty vectors: one slot per time step, filled with zeros for now S <- numeric(n_steps) I <- numeric(n_steps) R <- numeric(n_steps) # Fill in the starting values S[1] <- S0 I[1] <- I0 R[1] <- R0_start cat("Storage ready. First values: S =", S[1], " I =", I[1], " R =", R[1]) ``` Think of these vectors as the model's memory: they store the entire trajectory of the epidemic, one row per time step. --- **Step 5: Run the Euler loop.** This is where the model actually runs. At each time step, we calculate the rates of change using the SIR equations and use them to update the compartment values. Each line in the loop maps directly to an equation. ```{webr-r} for (t in 1:(n_steps - 1)) { # Rate of change: these ARE the SIR equations, written as code dS <- -beta * I[t] * S[t] # new infections leaving S dI <- beta * I[t] * S[t] - gamma * I[t] # gain infections, lose recoveries dR <- gamma * I[t] # recoveries entering R # Step forward: current value + (step size × rate of change) S[t + 1] <- S[t] + dt * dS I[t + 1] <- I[t] + dt * dI R[t + 1] <- R[t] + dt * dR } # Quick check: do the compartments still sum to 1 at the end? cat( "Final values: S =", round(S[n_steps], 4), " I =", round(I[n_steps], 4), " R =", round(R[n_steps], 4), "\n" ) cat("Sum:", round(S[n_steps] + I[n_steps] + R[n_steps], 6), "\n") ``` The loop is doing exactly what the Euler update formulas say, mechanically, one step at a time. This is the computational equivalent of doing the hand calculations from the table shown earlier, just thousands of times faster. ::: {.callout-note title="The loop is the model"} Notice that the three lines computing `dS`, `dI`, and `dR` are not arbitrary calculations. They are the SIR equations, written in R. Every biological assumption we made in the previous chapter is present in those three lines. The loop is not running the model alongside the equations, the loop is the equations. ::: --- **Step 6: Visualize the results.** The vectors `S`, `I`, and `R` now hold the full epidemic trajectory. We can plot them against time. ```{webr-r} # Plotting susceptible curve plot(times, S, type = "l", col = "steelblue", lwd = 2, ylim = c(0, 1), xlab = "Time", ylab = "Proportion", main = paste0( "SIR Epidemic (Euler) | β=", beta, " γ=", gamma, " R₀=", round(beta / gamma, 2) ) ) # Adding infectious and recovered curves lines(times, I, col = "firebrick", lwd = 2) lines(times, R, col = "seagreen", lwd = 2) # Adding a legend legend("right", legend = c("Susceptible", "Infectious", "Recovered"), col = c("steelblue", "firebrick", "seagreen"), lwd = 2, bty = "n" ) ``` ::: {.callout-tip title="Try this"} Go back to **Step 3** and try the following changes, re-running each step in order after each change: 1. Change `dt <- 0.1` to `dt <- 2.0`. What happens to the curves? This is the cost of a large step size. 2. Change the step size to `dt <- 0.001`. Did the code output results faster or slower than when using `dt <- 2.0`? (**NOTE:** How fast the code runs depends on the computer you are using right now, but in general note that the smaller our step size gets, the more accurate our results become at the cost of being increasingly computationally demanding.) ::: ## Solving Models Using the `pracma` Package in R {#sec-pracma} ::: {.callout-tip title="Learning Objectives"} By the end of @sec-pracma, you should be able to: - Explain what values the function `rk4sys` in `pracma` takes in. - Translate an SIR model to a function and solve it using `rk4sys`. ::: Euler's method is honest: every step is visible, every assumption is explicit. Its weakness is that those small linear steps accumulate error, especially when the epidemic curve is changing quickly. The next method, RK4, is more sophisticated. It looks further ahead before committing to each step, which makes it more accurate. The tradeoff is that it is less transparent. The package `pracma` in R, short for "Practical Numerical Math Functions," provides a large repertoire of functions used in numerical analysis, linear algebra, numerical optimization, and differential equations [@pracma_package]. ::: {.callout-note} ## Is there a specific name for the type of differential equations we are solving? The SIR models this module covers deal with systems of *first-order ordinary differential equations*. Breaking up each word in this description, we have: - **First-order:** This just means we are only dealing with equations involving the first derivative and not subsequent derivatives (the second, third, fourth, etc. derivative) - **Ordinary:** The unknown function which we are trying to solve for depends only one independent variable. ::: Although `pracma` contains many mathematical functions, we are primarily interested in `rk4sys`, a function to solve systems of differential equations using an advanced method known as the *Classic Fourth-Order Runge-Kutta Method* (or RK4 for short). This method for solving differential equations is much more complex than Euler's method but outputs solutions to a much greater accuracy. Like Euler's method, RK4's accuracy depends on the step size. However, in general, a solution to a model found using RK4 will be more accurate than one found using Euler's method, even when the same step size is used. :::{.callout-note collapse="true"} ## How exactly does RK4 work? In simple terms, RK4 estimates the value of a solution at the next point by taking a weighted average of four slope calculations, each computed at different points within the step interval. The Wikipedia article on Runge-Kutta methods goes more into the specific mathematical formulation of RK4 [@wiki_rk4]. ::: The following table summarizes the values `rk4sys` takes in: | **Variable Name** | **What Does It Take In?** | |-------------|-----------------------------------------------------------| | `f` | A user-defined function that defines the system of differential equations to be solved. | | `a` | A number specifying the starting point of the interval we want solutions for. | | `b` | A number specifying the endpoint of the interval we want solutions for. | | `y0` | A vector specifying the initial values of the system. A system of $m$ differential equations must have a `y0` of length $m$. | | `n` | A number specifying the step size (i.e. the number of steps that should be taken from `a` to `b`. | The function `rk4sys` returns a matrix. When that matrix is converted to a data frame using the function `as.data.frame`, its first column provides the time points at which the SIR model was solved and subsequent columns are the solutions to the model's state variables at those times. The state variables for an SIR model are the number or proportion of individuals susceptible, infectious, and recovered at a given time. Before we can call `rk4sys`, let's go over how to code the function to be passed into the parameter `f`. The function below is specifying a classic SIR model on the proportion of individuals in each compartment (see [Section @sec-the-classic-sir-epidemic-model]). **Step 1: Code the model function.** ```{webr-r} # Function Inputs: A time and vector of the model proportions at that time # Function Outputs: The rates of change for S, I, R at the input time ModelFunc <- function(time, state) { # Extracting state variables: susceptible, infected, and recovered S <- state[1] I <- state[2] R <- state[3] # Specifying the model parameters # We only have two parameters - the transmission and recovery rate # These parameters will change depending on the model beta <- 1.5 gamma <- 0.5 # Constructing our SIR model equations dS <- -(beta * I * S) dI <- (beta * I * S) - (gamma * I) dR <- (gamma * I) # Now we will return those model equations as a vector return(c(dS, dI, dR)) } # Example call ModelFunc(1, c(S = 0.99, I = 0.01, R = 0)) ``` --- **Step 2: Specify the parameters needed to call `rk4sys`.** ```{webr-r} library(pracma) # The start and ending points of the solution interval start <- 0 end <- 20 # Initial values to be passed into the parameter y # Make sure values line up with the state variables in the model function initial_vals <- c(S = 0.99, I = 0.01, R = 0) # How many steps to take from a to b step_size <- 10 # Printing the times at which solutions will be found print(seq(start, end, by = (end - start) / step_size)) ``` We now have everything we need to call `rk4sys`. --- **Step 3: Call `rk4sys` using our defined values.** ```{webr-r} library(pracma) # Make sure the values we defined are passed into the correct arguments output <- rk4sys( f = ModelFunc, a = start, b = end, y0 = initial_vals, n = step_size ) # Converting the output matrix to a data frame and printing it output <- as.data.frame(output) output ``` Now that we have a data frame of the solved SIR model equations at specified times, let's plot it. --- **Step 4: Plot the results.** ```{webr-r} times <- output$x S <- output$y.S I <- output$y.I R <- output$y.R beta <- 1.5 gamma <- 0.5 plot(times, S, type = "l", col = "steelblue", lwd = 2, ylim = c(0, 1), xlab = "Time", ylab = "Proportion", main = paste0( "SIR Epidemic (Runge-Kutta) | β=", beta, " γ=", gamma, " R₀=", round(beta / gamma, 2) ) ) lines(times, I, col = "firebrick", lwd = 2) lines(times, R, col = "seagreen", lwd = 2) legend("right", legend = c("Susceptible", "Infectious", "Recovered"), col = c("steelblue", "firebrick", "seagreen"), lwd = 2, bty = "n" ) ``` ::: {.callout-tip title="Try this"} Go back to **Step 2** and try the following changes, re-running each step in order after each change: 1. Change `step_size <- 10` to `step_size <- 50`. What changes about the final plot? Is this similar to what happened when we changed the step size in our code for Euler's method (see @sec-euler-code)? ::: ## Evaluating and Comparing Solving Methods {#sec-evaluating-and-comparing-both-methods} ::: {.callout-tip title="Learning Objectives"} By the end of @sec-evaluating-and-comparing-both-methods, you should be able to: - Explain the benefits and drawbacks of Euler's method and `rk4sys` in relation to solving SIR models. - Identify which solving method is more suitable for a particular SIR model or application. - Visualize the difference between solutions produced using Euler's method and `rk4sys` as step size changes. ::: The following table summarizes some potential benefits and drawbacks to using Euler's method or `rk4sys` in `pracma`: | **Solving Method** | **Potential Benefits** | **Potential Drawbacks** | |--------------|-------------------------|------------------------------| | Euler's Method | • Simple and intuitive code implementation • Easy to understand and visualize step-by-step • No external packages required • Can be implemented in Excel or by hand | • Low accuracy unless a very small step size is used • Computationally inefficient for small step sizes or complex models • Can be tedious or error-prone to calculate by hand for many steps • Can accumulate large numerical errors quickly over long periods| | `rk4sys` | • Uses a sophisticated solving method (Fourth-Order Runge-Kutta) that has a high accuracy • Handles large models with many system equations well • Preferred for complex models / advanced applications • Accuracy scales better with step size compared to Euler's method | • Less transparent: not much understanding of what is going on "under the hood" • Requires the external package `pracma` • Output accuracy depends on step size • Each step is more computationally expensive than Euler's method since four derivative approximations are required| :::{.callout-tip title="Tip"} Although the accuracy of the solutions output by Euler's method and the Fourth-Order Runge-Kutta method rely on step size, Runge-Kutta typically outputs more accurate results than Euler's method even with the same step size. This means we can "get away" with using a smaller step size when using `rk4sys` over Euler's method (thus saving us some computational cost / time). ::: The following simulation shows the difference between the solution curves produced by Euler's method and `rk4sys` for different values of the model parameters and step size: <iframe src="https://a-tom.shinyapps.io/euler_vs_rk_simulation/" width="100%" height="700" style="border:none;"></iframe> :::{.callout-tip title="Try this"} 1. Try changing the step size from $0.25$ to $0.1$. Do the curves solved using Euler's method and Runge-Kutta method become closer or further apart? Can you explain why mathematically? 2. Try changing the step size from $0.5$ to $0.1$. Which set of curves (the ones solved using Euler's method or Runge-Kutta method) change more? Can you explain why mathematically (think about which solving method scales better with step size in terms of accuracy)? ::: :::{.callout-important title="This isn't the only way!"} There are many packages in R that can solve SIR equations, and the package `deSolve` is a really common choice among them. `deSolve` uses very advanced and versatile solving methods that work well with large or complex SIR models. See the CRAN documentation for more details [@desolve_package]. ::: We now have two tools for solving SIR models. The next question is: how do we know which values of $\beta$ and $\gamma$ to use? That is the problem of fitting a model to data, and it is the topic of the next chapter. ## Chapter References {.unnumbered .unlisted}

Time	Approximations	Calculation Process
\(t=0\)	\(\textcolor{blue}{s(0)=99.0}\) \(\textcolor{red}{i(0)=1.0}\) \(\textcolor{green}{r(0)=0.0}\)	—
\(t=1\)	\(\textcolor{blue}{s(1)=97.5}\) \(\textcolor{red}{i(1)=2.0}\) \(\textcolor{green}{r(1)=0.5}\)	\(s(1)=\textcolor{blue}{99.0}-(1.0)(1.5)\left(\frac{\textcolor{red}{1.0}}{100.0}\right)(\textcolor{blue}{99.0})\) \(i(1)=\textcolor{red}{1.0}+(1.0)\left[(1.5)\left(\frac{\textcolor{red}{1.0}}{100.0}\right)(\textcolor{blue}{99.0})-0.5(\textcolor{red}{1.0})\right]\) \(r(1)=\textcolor{green}{0.0}+(1.0)(0.5)(\textcolor{red}{1.0})\)
\(t=2\)	\(\textcolor{blue}{s(2)=94.6}\) \(\textcolor{red}{i(2)=3.9}\) \(\textcolor{green}{r(2)=2.0}\)	\(s(2)=\textcolor{blue}{97.5}-(1.0)(1.5)\left(\frac{\textcolor{red}{2.0}}{100.0}\right)(\textcolor{blue}{97.5})\) \(i(2)=\textcolor{red}{2.0}+(1.0)\left[(1.5)\left(\frac{\textcolor{red}{2.0}}{100.0}\right)(\textcolor{blue}{97.5})-0.5(\textcolor{red}{2.0})\right]\) \(r(2)=\textcolor{green}{0.5}+(1.0)(0.5)(\textcolor{red}{2.0})\)
\(t=3\)	\(\textcolor{blue}{s(3)=89.1}\) \(\textcolor{red}{i(3)=7.5}\) \(\textcolor{green}{r(3)=3.4}\)	\(s(3)=\textcolor{blue}{94.6}-(1.0)(1.5)\left(\frac{\textcolor{red}{3.9}}{100.0}\right)(\textcolor{blue}{94.6})\) \(i(3)=\textcolor{red}{3.9}+(1.0)\left[(1.5)\left(\frac{\textcolor{red}{3.9}}{100.0}\right)(\textcolor{blue}{94.6})-0.5(\textcolor{red}{3.9})\right]\) \(r(3)=\textcolor{green}{2.0}+(1.0)(0.5)(\textcolor{red}{3.9})\)
\(t=4\)	\(\textcolor{blue}{s(4)=79.1}\) \(\textcolor{red}{i(4)=13.7}\) \(\textcolor{green}{r(4)=7.2}\)	\(s(4)=\textcolor{blue}{89.1}-(1.0)(1.5)\left(\frac{\textcolor{red}{7.5}}{100.0}\right)(\textcolor{blue}{89.1})\) \(i(4)=\textcolor{red}{7.5}+(1.0)\left[(1.5)\left(\frac{\textcolor{red}{7.5}}{100.0}\right)(\textcolor{blue}{89.1})-0.5(\textcolor{red}{7.5})\right]\) \(r(4)=\textcolor{green}{3.4}+(1.0)(0.5)(\textcolor{red}{7.5})\)
\(t=5\)	\(\textcolor{blue}{s(5)=62.8}\) \(\textcolor{red}{i(5)=23.1}\) \(\textcolor{green}{r(5)=14.0}\)	\(s(5)=\textcolor{blue}{79.1}-(1.0)(1.5)\left(\frac{\textcolor{red}{13.7}}{100.0}\right)(\textcolor{blue}{79.1})\) \(i(5)=\textcolor{red}{13.7}+(1.0)\left[(1.5)\left(\frac{\textcolor{red}{13.7}}{100.0}\right)(\textcolor{blue}{79.1})-0.5(\textcolor{red}{13.7})\right]\) \(r(5)=\textcolor{green}{7.2}+(1.0)(0.5)(\textcolor{red}{13.7})\)

3.1 Introduction to Solving SIR Models

3.2 Euler’s Method

3.2.1 Euler’s Method Code Demonstration

3.3 Solving Models Using the pracma Package in R

3.4 Evaluating and Comparing Solving Methods

Chapter References

3.3 Solving Models Using the `pracma` Package in R