Section 1.4 Analyzing Equations Numerically
Objectives
To understand that numerical algorithms such as Euler's method allow the approximation of solutions to the initial value problems and that there are more efficient algorithms than Euler's method such as those algorithms that use the RungeKutta methods.
To understand that Taylor's Theorem is a very useful tool for studying differential equations.
To understand that error analysis of the rate of convergence is very important for any numerical algorithm.
Just as numerical algorithms are useful when finding the roots of polynomials, numerical methods will prove very useful in our study of ordinary differential equations. Consider the polynomial \(f(x) = x^2  2\text{.}\) We do not need a numerical algorithm to see that the roots of this polynomial are \(x = \sqrt{2}\) and \(x =  \sqrt{2}\text{.}\) However, a numerical method such as the NewtonRaphson Algorithm is very useful for approximating \(\sqrt{2}\) as a decimal.^{â€‰1â€‰} Similarly, it may be easier to generate a numerical solution for differential equations if our goal is simply to plot a solution. In addition, there will be differential equations for which it is impossible to find a solution in terms of elementary functions such as polynomials, trigonometric functions, and exponential functions.
Subsection 1.4.1 Euler's Method
Suppose that we wish to solve the initial value problem
The equation \(y' = y + t\) is not separable, which currently is the only analytic technique at our disposal. However, we can try to find a numerical approximation for the solution. A numerical approximation is simply a table (possibly very large) of \(t\) and \(y\) values.
We will attempt to find a numerical solution for (1.4.1)â€“(1.4.2) on the interval \([0, 1]\text{.}\) Even with the use of a computer, we cannot approximate the solution at every single point on an interval. For the initial value problem
we might be able to find approximations at \(a = t_0, t_1, t_2, \ldots, t_N = b\) in \([a, b]\) at best. If we choose \(t_1, t_2, \ldots, t_N\) to be equally spaced on \([a, b]\text{,}\) we can write
where \(h = 1/N\) and \(k = 1, 2, \ldots, N\text{.}\) We say that \(h\) is the step size for our approximation.
Given an approximation \(Y_k\) for the solution \(y_k = y(t_k)\text{,}\) the question is how to find an approximate solution \(Y_{k+1}\) at \(t_{k+1}\text{.}\) To generate the second approximation, we will construct a tangent line to the solution at \(y(t_0) = y_0\text{.}\) If we use the slope of the solution curve at \(t_0\text{,}\) then
Making use of the fact that
or equivalently
the estimate for our solution at \(t_1 = t_0 + h\) is
Similarly, the approximation at \(t_2 = t_0 + 2h\) will be
Our general algorithm is
The idea is to compute tangent lines at each step and use this information to get our next approximation.
The algorithm that we have described is known as Euler's method. Let us estimate a solution to (1.4.1)â€“(1.4.2) on the interval \([0, 1]\) with step size \(h = 0.1\text{.}\) Since \(y(0) = 1\text{,}\) we can make our first approximation exact,
To generate the second approximation, we will construct a tangent line to the solution at \(y(0) = 1\text{.}\) If we use the slope of the solution curve at \(t_0 = 0\text{,}\)
and make use of the fact that
the estimate for our solution at \(t = 0.1\) is
Similarly, the approximation at \(t = 0.2\) will be
Our general algorithm is
The initial value problem (1.4.1)â€“(1.4.2) is, in fact, solvable analytically with solution \(y(t) = 2e^t  t  1\text{.}\) We can compare our approximation to the exact solution in TableÂ 1.4.1. We can also see graphs of the approximate and exact solutions in FigureÂ 1.4.2. Notice that the error grows as we get further away from our initial value. In fact, the graph of the approximation for \(h = 0.001\) is obscured by the graph of the exact solution. In addition, a smaller step size gives us a more accurate approximation (TableÂ 1.4.3).
\(k\)  \(t_k\)  \(Y_k\)  \(y_k\)  \(y_k  Y_k\)  Percent Error 
0  0.0  1.0000  1.0000  0.0000  0.00% 
1  0.1  1.1000  1.1103  0.0103  0.93% 
2  0.2  1.2200  1.2428  0.0228  1.84% 
3  0.3  1.3620  1.3997  0.0377  2.69% 
4  0.4  1.5282  1.5836  0.0554  3.50% 
5  0.5  1.7210  1.7974  0.0764  4.25% 
6  0.6  1.9431  2.0442  0.1011  4.95% 
10  1.0  3.1875  3.4366  0.2491  7.25% 
\(t_k\)  \(h = 0.1\)  \(h = 0.02\)  \(h = 0.001\)  Exact Solution 
0.1  1.1000  1.1082  1.1102  1.1103 
0.2  1.2200  1.2380  1.2426  1.2428 
0.3  1.3620  1.3917  1.3993  1.3997 
0.4  1.5282  1.5719  1.5831  1.5836 
0.5  1.7210  1.7812  1.7966  1.7974 
0.6  1.9431  2.0227  2.0431  2.0442 
0.7  2.1974  2.2998  2.3261  2.3275 
0.8  2.4872  2.6161  2.6493  2.6511 
0.9  2.8159  2.9757  3.0170  3.1092 
1.0  3.1875  3.3832  3.4238  3.4366 
Activity 1.4.1. Euler's Method and Error.
Consider the initial value problem
(a)
Use separation of variables to solve the initial value problem.
(b)
Compute \(y(x)\) for \(x = 0, 0.2, 0.4, \ldots, 1\text{.}\)
(c)
Use Euler's method to approximate solutions to the initial value problem for \(x = 0, 0.2, 0.4, \ldots, 1\text{.}\)
(d)
Compare the exact values of the solution (TaskÂ 1.4.1.b) to the approximate values of the solution (TaskÂ 1.4.1.c) and comment on what happens as \(x\) varies from \(0\) to \(1\text{.}\)
Subsection 1.4.2 Finding an Error Bound
To fully understand Euler's method, we will need to recall Taylor's theorem from calculus.
Theorem 1.4.4.
If \(x \gt x_0\text{,}\) then
where \(\xi \in (x_0, x)\text{.}\)
Given the initial value problem
choose \(t_1, t_2, \ldots, t_N\) to be equally spaced on \([t_0, a]\text{,}\) we can write
where \(h = (a  t_0)/N\) and \(k = 1, 2, \ldots, N\text{.}\) Taylor's Theorem tells us that
If we know the values of \(y\) and its derivatives at \(t_k\text{,}\) then we can determine the value of \(y\) at \(t_{k + 1}\text{.}\)
The simplest approximation can be obtained by taking the first two terms of the Taylor series. That is, we will use a linear approximation,
This gives us Euler's method,
The terms that we are omitting, all contain powers of \(h\) of at least degree two. If \(h\) is small, then \(h^n\) for \(n = 2, 3, \ldots\) will be very small and these terms will not matter much.
We can actually estimate the error incurred by Euler's method if we make use of Taylor's Theorem.
Theorem 1.4.5.
Let \(y\) be the unique solution to the initial value problem
where \(t \in [a, b]\text{.}\) Suppose that \(f\) is continuous and there exists a constant \(L \gt 0\) such that
whenever \((t, y_1)\) and \((t, y_2)\) are in \(D = [a, b] \times {\mathbb R}\text{.}\) Also assume that there exists an \(M\) such that
for all \(t \in [a, b]\text{.}\) If \(Y_0, \ldots, Y_N\) are the approximations generated by Euler's method for some positive integer \(N\text{,}\) then
The condition that there exists a constant \(L \gt 0\) such that
whenever \((t, y_1)\) and \((t, y_2)\) are in \(D = [a, b] \times {\mathbb R}\) is called a Lipschitz condition. Many of the functions that we will consider satisfy such a condition. If the condition is satisfied, we can usually say a great deal about the function.
\(k\)  \(t_k\)  \(Y_k\)  \(y_k = y(t_k)\)  \(y_k  Y_k\)  Estimated Error 
0  0.0  1.0000  1.0000  0.0000  0.0000 
1  0.1  1.1000  1.1103  0.0103  0.0286 
2  0.2  1.2200  1.2428  0.0228  0.0602 
3  0.3  1.3620  1.3997  0.0377  0.0951 
4  0.4  1.5282  1.5836  0.0554  0.1337 
5  0.5  1.7210  1.7974  0.0764  0.1763 
6  0.6  1.9431  2.0442  0.1011  0.2235 
7  0.7  2.1974  2.3275  0.1301  0.2756 
8  0.8  2.4872  2.6511  0.1639  0.3331 
9  0.9  2.8159  3.1092  0.2033  0.3968 
10  1.0  3.1875  3.4366  0.2491  0.4671 
We can now compare the estimated error from our theorem to the actual error of our example. We first need to determine \(M\) and \(L\text{.}\) Since
we can take \(L\) to be one. Since \(y'' = 2e^t\text{,}\) we can bound \(y''\) on the interval \([0, 1]\) by \(M = 2e\text{.}\) Thus, we can bound the error by
for \(h =0.1\text{.}\) Our results are in TableÂ 1.4.6.
Subsection 1.4.3 Improving Euler's Method
If we wish to improve upon Euler's method, we could add more terms of Taylor series. For example, we can obtain a more accurate approximation by using a quadratic Taylor polynomial,
However, we need to know \(y''(t_0)\) in order to use this approximation. Using the chain rule from multivariable calculus, we can differentiate both sides of \(y' = f(t, y)\) to obtain
Thus, our approximation becomes
The problem is that some preliminary analytic work must be done. That is, before we can write a program to compute our solution, we must find \(\partial f/\partial t\) and \(\partial f / \partial y\text{,}\) although this is less of a problem with the availability of computer algebra systems such as Sage.
Around 1900, two German mathematicians, Carle Runge and Martin Kutta, independently invented several numerical algorithms to solve differential equations. These methods, known as RungeKutta methods, estimate the higherorder terms of the Taylor series to find an approximation that does not depend on computing derivatives of \(f(t, y)\text{.}\)
If we consider the initial value problem
then
or
by the Fundamental Theorem of Calculus. In Euler's method, we approximate the righthand side of (1.4.3) by
In terms of the definite integral, this is simply a lefthand sum. In the improved Euler's method or the secondorder RungeKutta method we will estimate the righthand side of (1.4.3) using the trapezoid rule from calculus,
Thus, our algorithm becomes
However, we have a problem since \(y_1\) appears in the righthand side of our approximation. To get around this difficulty, we will replace \(y_1\) in the righthand side of (1.4.4) with the Euler approximation for \(y_1\text{.}\) Thus,
To understand that the secondorder RungeKutta method is actually an improvement over the traditional Euler's method, we will need to use the Taylor approximation for a function of two variables. Let us assume that \(f(x,y)\) is defined on some rectangle and that all of the derivatives of \(f\) are continuously differentiable. Then
As in the case of the single variable Taylor series, we can write a Taylor polynomial if the Taylor series is truncated,
where the second term is the remainder term and \((\overline{x}, \overline{y} )\) lies on the line segment joining \((x, y)\) and \((x + h, y + k)\text{.}\)
In the Improved Euler's Method, we adopt a formula
where
That is,
The idea is to choose the constants \(w_1\text{,}\) \(w_2\text{,}\) \(\alpha\text{,}\) and \(\beta\) as accurately as possible in order to duplicate as many terms as possible in the Taylor series
We can make equations (1.4.5) and (1.4.6) agree if we choose \(w_1 = 1\) and \(w_2 = 0\text{.}\) Since \(y' = f\text{,}\) we obtain Euler's method.
If we are more careful about choosing our parameters, we can obtain agreement up through the \(h^2\) term. If we use the two variable Taylor series to expand \(f(t + \alpha h, y + \beta h f)\text{,}\) we have
where \({\mathcal O}(h^2)\) means that of the subsequent terms have a factor of \(h^n\) with \(n \geq 2\text{.}\) Using this expression, we obtain a new form for (1.4.5),
Since \(y'' = f_t + f_y f\) by the chain rule, we can rewrite (1.4.6) as
We can make equations (1.4.7) and (1.4.8) agree up through the quadratic terms if we require that
If we choose \(\alpha = \beta = 1\) and \(w_1 = w_2 = 1/2\text{,}\) these equations are satisfied, and we obtain the improved Euler's method
The improved Euler's method or the secondorder RungeKutta method is a more sophisticated algorithm that is less prone to error due to the step size \(h\text{.}\) Euler's method is based on truncating the Taylor series after the linear term. Since
we know that the error depends on \(h\text{.}\) On the other hand, the error for the improved Euler's method depends on \(h^2\text{,}\) since
If we use Simpson's rule to estimate the integral in
we can improve our accuracy up to \(h^4\text{.}\) The idea is exactly the same, but the algebra becomes much more tedious. This method is known as the RungeKutta method of order 4 and is given by
where
Subsection 1.4.4 Important Lessons

We can use Euler's method to find an approximate solution to the initial value problem
\begin{align*} y' & = f(t, y),\\ y(a) & = \alpha \end{align*}on an interval \([a, b]\text{.}\) If we wish to find approximations at \(N\) equally spaced points \(t_1, \ldots, t_N\text{,}\) where \(h = (ba)/N\) and \(t_i = a + i h\text{,}\) our approximations should be
\begin{align*} Y_0 & = \alpha,\\ Y_1 & = Y_0 + h f(\alpha, Y_0),\\ Y_2 & = Y_1 + h f(t_1, Y_1,)\\ & \vdots\\ Y_{k+1} & = Y_k + h f(t_k, Y_k),\\ Y_N & = Y_{N1} + h f(t_{N1}, Y_{N1}). \end{align*}In practice, no one uses Euler's method. The RungeKutta methods are better algorithms.

Taylor's Theorem is a very useful tool for studying differential equations. If \(x \gt x_0\text{,}\) then
\begin{equation*} f(x) = f(x_0) + f'(x_0)(x  x_0) + \frac{f''(x_0)}{2!} (x  x_0)^2 + \cdots + \frac{f^{(n)}(\xi)}{n!}(x  x_0)^n, \end{equation*}where \(\xi \in (x_0, x)\text{.}\)

Error analysis rate of convergence is very important for any numerical algorithm. Our approximation is more accurate for smaller values of \(h\text{.}\) Under reasonable conditions we can also bound the error by
\begin{equation*} y(t_i)  Y_i \leq \frac{hM}{2L} [ e^{L(t_i  a)}  1], \end{equation*}where \(y\) is the unique solution to the initial value problem
\begin{align*} y' & = f(t, y),\\ y(a) & = \alpha. \end{align*} 
The condition that there exists a constant \(L \gt 0\) such that
\begin{equation*} f(t, y_1)  f(t, y_2) \leq L y_1  y_2, \end{equation*}whenever \((t, y_1)\) and \((t, y_2)\) are in \(D = [a, b] \times {\mathbb R}\) is called a Lipschitz condition.
Using Taylor series, we can develop better numerical algorithms to compute solutions of differential equations. The RungeKutta methods are an important class of these algorithms.

The improved Euler's method is given by
\begin{equation*} y(t + h) = y(t) + \frac{h}{2} f(t, h) + \frac{h}{2} f(t + h, y + h f(t, y)) \end{equation*}with the error bound depending on \(h^2\text{.}\)

The RungeKutta method of order 4 is given by
\begin{equation*} y(t + h) = y(t) + \frac{1}{6} (F_1 + 2 F_2 + 2 F_3 + F_4), \end{equation*}where
\begin{align*} F_1 & = h f(t, y)\\ F_2 & = hf\left( t + \frac{1}{2} h, y + \frac{1}{2} F_1 \right)\\ F_3 & = hf\left( t + \frac{1}{2} h, y + \frac{1}{2} F_2 \right)\\ F_4 & = hf(t + h, y + F_3) \end{align*}with the error bound depending on \(h^4\text{.}\)
Reading Questions 1.4.5 Reading Questions
1.
We can use Taylor polynomials to approximate a function \(f(x)\) near a point \(x_0\text{.}\) Explain why this approximation can only be expected to be accurate near \(x_0\text{.}\)
2.
Should we always use Euler's method when approximating a solution to an initial value problem? Why or why not?
Exercises 1.4.6 Exercises
Finding Solutions.
For each of the initial value problem
in Exercise GroupÂ 1.4.6.1â€“6,
Write the Euler's method iteration \(Y_{k+1} = Y_k + h f(t_k, Y_k)\) for the given problem, identifying the values \(t_0\) and \(y_0\text{.}\)
Using a step size of \(h = 0.1\text{,}\) compute the approximations for \(Y_1\text{,}\) \(Y_2\text{,}\) and \(Y_3\text{.}\)
Solve the problem analytically if possible. If it is not possible for you to find the analytic solution, use Sage.
Use the results of (c) and (d), to construct a table of errors for \(Y_i  y_i\) for \(i = 1, 2, 3\text{.}\)
1.
\(y' = 2y\text{,}\) \(y(0) = 0\)
2.
\(y' = ty\text{,}\) \(y(0) = 1\)
3.
\(y' = y^3\text{,}\) \(y(0) = 1\)
4.
\(y' = y\text{,}\) \(y'(0) = 1\)
5.
\(y' = y + t\text{,}\) \(y'(0) = 2\)
This equation is a firstorder linear equation (SectionÂ 1.5), but it is possible to find the analytic solution using Sage (SubsectionÂ 1.2.10).
6.
\(y' = 1/y\text{,}\) \(y'(0) = 2\)
7.
Consider the differential equation
with initial value \(y(0) = 2\text{.}\)
Find the exact solution of the initial value problem.
Use Euler's method with step size \(h = 0.5\) to approximate the solution to the initial value problem on the interval \([0,2]\) Your solution should include a table of approximate values of the dependent variable as well as the exact values of the dependent variable. Make sure that your approximations are accurate to four decimal places.
Sketch the graph of the approximate and exact solutions.
Use the error bound theorem (TheoremÂ 1.4.5) to estimate the error at each approximation. Your solution should include a table of approximate values of the dependent variable the exact values of the dependent variable, the error estimates, and the actual error. Make sure that your approximations are accurate to four decimal places.
8.
In this series of exercises, we will prove the error bound theorem for Euler's method (TheoremÂ 1.4.5).

Use Taylor's Theorem to show that for all \(x \geq 1\) and any positive \(m\text{,}\)
\begin{equation*} 0 \leq (1 + x)^m \leq e^{mx}. \end{equation*} 
Use part (1) and geometric series to prove the following statement: If \(s\) and \(t\) are positive real numbers, and \(\{ a_i \}_{i = 0}^k\) is a sequence satisfying
\begin{align*} a_0 & \geq \frac{t}{s}\\ a_{i+1} & \leq (1 + s) a_i + t, \text{ for }i = 0, 1, 2, \ldots, k, \end{align*}then
\begin{equation*} a_{i + 1} \leq e^{(i + 1)s} \left( \frac{t}{s} + a_0 \right)  \frac{t}{s}. \end{equation*} 
When \(i = 0\text{,}\) \(y(t_0) = Y_0 = \alpha\) and the theorem is true. Use Euler's method and Taylor's Theorem to show that
\begin{equation*} y_{i + 1}  Y_{i + 1}  \leq y_i  Y_i + h f(t_i, y_i)  f(t_i, Y_i) + \frac{h^2}{2} y''(\xi_i) \end{equation*}where \(\xi_i \in (t_i, t_{i+1})\) and \(y_k = y(t_k)\text{.}\)

Show that
\begin{equation*} y_{i + 1}  Y_{i + 1}  \leq y_i  Y_i (1 + hL) + \frac{Mh^2}{2}. \end{equation*} 
Let \(a_j = y_j  Y_j\) for \(j = 0, 1, \ldots, N\text{,}\) \(s = hL\text{,}\) and \(t = M h^2/2\) and apply part (2) to show
\begin{equation*} y_{i + 1}  Y_{i + 1}  \leq e^{(1 + i)hL} \left( y_0  Y_0 + \frac{M h^2 }{2hL} \right)  \frac{M h^2}{2hL}. \end{equation*}and derive that
\begin{equation*} y_{i + 1}  Y_{i + 1}  \leq \frac{M h }{2L} \left( e^{(t_{i+ 1}  a)L}  1\right) \end{equation*}for each \(i = 0, 1, \ldots, N1\text{.}\)
Hints for part (2):

For fixed \(i\) show that
\begin{align*} a_{i + 1} & \leq (1 + s)a_i + t\\ & \leq (1 + s)[(1 + s)a_{i  1}+ t] + t\\ & \leq (1 + s)\{ (1 + s)[(1 + s) a_{i  2} + t]+ t\} + t\\ & \vdots\\ & \leq (1 + s)^{i+1}a_0 + [1 + (1 + s) + (1+s)^2 + \cdots + (1 + s)^i]t. \end{align*} 
Now use a geometric sum to show that
\begin{equation*} a_{i + 1} \leq (1 + s)^{i + 1} a_0 + \frac{t}{s}[(1 + s)^{i + 1}  1] = (1 + s)^{i + 1} \left( \frac{t}{s} + a_0 \right)  \frac{t}{s}. \end{equation*} 
Apply part (1) to derive
\begin{equation*} a_{i + 1} \leq e^{(i + 1)s} \left( \frac{t}{s} + a_0 \right)  \frac{t}{s}. \end{equation*}
Subsection 1.4.7 Sageâ€”Numerical Routines for solving ODEs
Not all differential equations can be solved using algebra and calculus even if we are very clever. If we encounter an equation that we cannot solve or use Sage to solve, we must resort to numerical algorithms like Euler's method or one of the RungeKutta methods, which are best implemented using a computer. Fortunately, Sage has some very good numerical solvers. Sage will need to know the following to solve a differential equation:
An abstract function,
A differential equation, including an initial condition,
A Sage command to solve the equation.
Suppose we wish to solve the initial value problem \(dy/dx = x + y\text{,}\) \(y(0) = 1\text{.}\) We can use Sage to find an algebraic solution.
We can also use Euler's method to find a solution for our initial value problem.
The syntax of eulers_method
for the inital value problem
with step size \(h\) on the interval \([x_0, x_1]\) is eulers_method(f, x0, y0, h, x1)
Notice that we obtained a table of values. However, we can use the line
command from Sage to plot the values \((x, y)\text{.}\)
As we pointed out, eulers_method
is not very sophisticated. We have to use a very small step size to get good accuracy, and the method can generate errors if we are not careful. Fortunately, Sage has much better algorithms for solving initial value problems. One such algorithm is desolve_rk4
, which implements the fourth order RungeKutta method.
Again, we just get a list of points. However, desolve_rk4
has some nice builtin graphing utilities.
Write the Sage commands to compare the graphs obtained using eulers_method
and desolve
with the exact solution.
Not only is desolve_rk4
more accurate, it is much easier to use. For more information, see www.sagemath.org/doc/reference/calculus/sage/calculus/desolvers.html
.