Math 555: Differential Equations I

Lecture Companion

§2.7: Numerical Methods

(see Numerical Project)

Consider a first order initial value problem $$ \begin{cases} y' = f(t,y), \\ y(t_0) = y_0. \end{cases} $$ We wish to approximate the solution curve using Euler's method, Improved Euler's method, and the Runge-Kutta method.

The general idea of each of these methods is to partition the $t$-axis into intervals of length $h > 0$, then use the slope field to approximate the $y$-values of the solution curve above each $t$-value. The output of each method is thus a list of points in the $ty$-plane whose $t$-values are evenly spaced. We interpret the piecewise linear graph that connects these points as an approximation of the solution curve.

We begin by stating the algorithms. The required inputs for each of the methods are:

${\tt f}$	function of $({\tt t, y})$ representing the right-hand side of the DE
${\tt t_0}$	initial $t$-value
${\tt y_0}$	initial $y$-value
${\tt h}$	step-size in $t$-direction
${\tt N}$	number of steps

The $t$-values for each method are the same. For $1 \leq {\tt k} \leq {\tt N}$, they are given by:

${\tt t_k = t_0 + k\cdot h}$

The $y$-values for each method are updated according to the following rules. Again, ${\tt k}$ varies from $1$ to ${\tt N}$.

Euler's Method:

${\tt m = f(t_{k-1},y_{k-1})}$

${\tt y_k = y_{k-1} + m\cdot h}$

Improved Euler's Method:

${\tt m_1 = f(t_{k-1},y_{k-1})}$
${\tt m_2 = f(t_k,y_{k-1} + m_1\cdot h)}$
${\tt m = \tfrac{1}{2}(m_1 + m_2)}$

${\tt y_k = y_{k-1} + m\cdot h}$

Runge-Kutta Method:

${\tt m_1 = f(t_{k-1},y_{k-1})}$
${\tt m_2 = f(t_{k-1} + \tfrac{1}{2}\cdot h,y_{k-1} + \tfrac{1}{2}\cdot m_1\cdot h)}$
${\tt m_3 = f(t_{k-1} + \tfrac{1}{2}\cdot h,y_{k-1} + \tfrac{1}{2}\cdot m_2\cdot h)}$
${\tt m_4 = f(t_k,y_{k-1} + m_3\cdot h)}$
${\tt m = \tfrac{1}{6}(m_1 + 2\cdot m_2 + 2\cdot m_3 + m_4)}$

${\tt y_k = y_{k-1} + m\cdot h}$

Notice that each of the numerical methods approximates the "next" $y$-value as ${\tt y_k = y_{k-1} + m\cdot h}$. The differences in the algorithms lie solely in how they compute the slope, ${\tt m}$.

In a later class you will study the derivations of these methods in detail. These methods have close connections with the approximate integration techniques of Calculus I and II. Euler's method corresponds to the Left Endpoint Rule, Improved Euler's method to the Trapezoidal Rule, and the Runge-Kutta method to Simpson's Rule.

Below, we consider one example done using three different choices of software: SageMath, Octave (or MATLAB), and LibreOffice Calc (or Microsoft Excel). You can run the SageMath and Octave code directly on this page; you'll need to download the XLSX-file and open it with an appropriate program on your computer to run that one.

Example. Consider the initial value problem, $$ \begin{cases} y' = -2ty + 2t, \\ y(0) = 10. \end{cases} $$ We wish to use the numerical methods described above to approximate the solution to this IVP, and compare them to the exact solution.

First, we must use the Fundamental Existence and Uniqueness Theorem to check that a unique solution does indeed exist, and that we are not wasting our time trying to implement these methods.

This equation is linear, and can be written in the form $$ y' + 2ty = 2t. $$ The coefficient functions are $p(t) = 2t$ and $g(t) = 2t$. Both are continuous at all values of $t$, hence at our initial value $t_0 = 0$. Thus a unique solution to this problem does exist.

Since this equation is linear, we can write down the solution. In integral form, it is $$ y(t) = e^{-t^2}\Big( \int_0^t e^{u^2}\, 2u\, du + 10\Big). $$ The integral can be computed without much trouble to give the solution, $$ y(t) = 1 + 9 e^{-t^2}. $$ We will use this exact solution as a reference to test the performance of the numerical methods.

SageMath

We begin the numerical approximation portion of these notes using SageMath. SageMath (or just Sage) is free, open-source software that is essentially just a collection of Python packages designed for doing computational mathematics.

In the code block below we define Python functions to carry out each numerical method. Each method's function takes in the information outlined above $({\tt f, t_0, y_0, h, N})$, and outputs a list of points and a corresponding piecewise linear graph connecting those points.

You must click the "Run" button below the code box below to initialize these functions.

Now we can use these functions to approximate the solution to our IVP. First, we can use Sage to solve the IVP exactly. We already did this, but it will be nice to have confirmation.

Now we apply the numerical methods with ${\tt h} = 0.1$ and ${\tt N} = 20$, and plot their graphs on the same set of axes together with the exact solution. Sometimes it is useful to plot all of these curves on the slope field for the differential equation, especially when you don't know the exact solution. This allows you to check if the numerical solutions "look" generally correct.

Octave

The next software that we will utilize is called Octave. Octave is free, open-source software that is designed to replicate MATLAB. All of the code in the cells below should run equally well in MATLAB. Any necessary tweaks in MATLAB should be minor, but please notify me of them if you find any.

Again, we begin by defining the methods as functions to be called later. Octave is not a symbolic language like Sage, so we will have to save and output our data slightly differently. Also, as of right now, all Octave code must be run in a single cell. I am working on fixing this, but for now all of the code is in the next cell.

LibreOffice

Finally, here is an XLSX-file that you can download and use on your computer. It is a rather simple spread sheet, so it should work properly in Excel, LibreOffice, and Google Sheets.

Back to main page

Your use of Wichita State University content and this material is subject to our Creative Common License.