Euler's method is an algorithm for approximating solutions to first-order IVPs \begin{equation*} \frac{dy}{dt} = f(t,y) \quad y(t_0) = y_0. \end{equation*} The key idea of Euler's method and all numerical techniques is discretization. If a variable can take on any value between two specified values, it is called a continuous variable; otherwise, it is called a discrete variable. When the independent variable is time, \(t\text{,}\) we discretize it by imagining it to correspond with the ticks of a clock, just like how the second hand on a clock moves jerkily as opposed to smoothly from one second to the next.

Our discrete ticks are determined by a start time, \(t_0\text{,}\) and a tick size, or step size \(h = \Delta t\text{,}\) for example \begin{equation*} t_0, \underbrace{t_0+h}_{t_1}, \underbrace{t_0+2h}_{t_2}, \underbrace{t_0+3h}_{t_3}, \ldots \end{equation*}

The next step is to pick a starting point or initial condition: \((t_0, y_0)\text{,}\) then we use \(f(t,y)\) to generate a list of coordinate pairs which can be plotted, i.e. \((t_0, y_0), (t_1,y_1), (t_2,y_2), (t_3,y_3), \ldots\)

The above list of pairs is generated via the following recursive algorithm.

Since Sage is really just Python with lots of extra built-in math related functions we can implement the above algorithm as a Python function. Be aware that Sage has two kinds of functions, mathematical functions and Python functions. A Python function is a way to store an algorithm or procedure that can be used repeatedly. You will need to evaluate the euler Python function below once and all the other cells in your worksheet can use it.

The Python version of the algorithm should be fairly easy to understand, but perhaps the last line needs explanation. It uses the Python function zip. This function takes two lists and combines them into one. The following Sage cell has an example of how zip works.

Note: You will need to define a mathematical function f(t,y) in Sage before using the Python euler function. This f(t,y) function is the right hand side of the differential equation \(y' = f(t,y)\text{.}\) This is done in Sage by typing, for example: f(t,y) = t*y

Evaluating the following Sage cell which contains the Python implementation of Euler's method will not result in any output. It does create a new function that all of the other Sage cells on this page depend upon, so you must evaluate it before evaluating the other cells on this page.

Let's test our Python implementation of Euler's method with the population model. Sage's table function displays the output in an easy to read table.

Of course, once we have a list of coordinate pairs, it is natural to plot them. If you want points and lines then you will need to use the list_plot function to plot points, and the line function to plot lines.

We can use Sage to plot the output from Euler's method and an exact solution on the same set of axes for comparison. Below we do this for the IVP: \begin{equation*} \frac{dy}{dt} = \sin t, \quad y(0)=0. \end{equation*} which has particular solution, \(y(t) = -\cos(t)+1\text{.}\)

It is instructive to see how the choice of step size affects the approximation. We can see this by plotting multiple approximation curves on the same set of axes. Below we do this for the linear IVP: \begin{equation*} \frac{dy}{dt} + y = \cos t, \quad y(0)=1. \end{equation*} Which has particular solution: \begin{equation*} y(t) = \frac{1}{2} \left( \cos t + \sin t + e^{-t} \right). \end{equation*}

We can combine appoximate solutions with slope field plots. Below we do so for the IVP: \begin{equation*} \dfrac{dy}{dt} = -y + \cos t, \quad y(-\pi)=5, \end{equation*} which has solution \begin{equation*} y(t) = \frac{1}{2} \left( \cos t + \sin t + 11e^{-(t+\pi)} \right). \end{equation*}

Euler's method uses the left side endpoint of the interval \([t_n, t_{n+1}]\) to predict the slope for the whole interval. But that is rather unsymmetric, because it ignores the slope at the right side of the interval. What if we compute the slopes at both endpoints and then average them?

Below is an implementation of the Heun algorithm in Python.

Next we compare Euler and Heun on the same set of axes for the population model IVP, \begin{equation*} \frac{dy}{dt} = y, \quad y(0)=1, \end{equation*} with solution, \(y(t) = e^t\text{.}\)