Analyzing Equations Numerically

Section 1.4 Analyzing Equations Numerically

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 Newton-Raphson Algorithm is very useful for approximating

$\sqrt{2}$ as a decimal.⁷ 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.

See any calculus text for a description of the Newton-Raphson Algorithm.

🔗

Subsection 1.4.1 Euler's Method

🔗

Suppose that we wish to solve the initial value problem

$\begin{align} y' &= f(t, y) = y + t, \label{firstlook04-equation-euler-1}\tag{1.4.1}\\ y(0) &= 1. \label{firstlook04-equation-euler-2}\tag{1.4.2} \end{align}$

🔗

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

$\begin{align*} y' &= f(t, y)\\ y(t_0) &= y_0, \end{align*}$

🔗

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

$\begin{equation*} t_k = t_0 + kh, \end{equation*}$

🔗

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 slope of the solution curve at

$t_0\text{,}$ then

$\begin{equation*} y'(t_0) = f(t_0, y_0). \end{equation*}$

🔗

Making use of the fact that

$\begin{equation*} \frac{y(t_0 + h) - y(t_0)}{h} \approx y'(t_0, y(t_0)) = y'(t_0, y_0) \end{equation*}$

🔗

or equivalently

$\begin{equation*} y(t_0 + h) = y(t_1) \approx y(t_0) + h y'(t_0, y_0), \end{equation*}$

🔗

the estimate for our solution at

$t_1 = t_0 + h$ is

$\begin{equation*} Y_1 = Y_0 + h f(t_0, Y_0). \end{equation*}$

🔗

Similarly, the approximation at

$t_2 = t_0 + 2h$ will be

$\begin{equation*} Y_2 = Y_1 + h f(t_1, Y_1). \end{equation*}$

🔗

Our general algorithm is

$\begin{equation*} Y_{k+1} = Y_k + h f(t_k, Y_k). \end{equation*}$

🔗

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,

$\begin{equation*} Y_0 = y(0) = 1. \end{equation*}$

🔗

To generate the second approximation, we will construct a tangent line to the solution at

$y(0) = 1\text{.}$ If we use slope of the solution curve at

$t_0 = 0\text{,}$

$\begin{equation*} y'(0) = f(y(0), 0) = y(0) + 0 = 1 + 0 = 1, \end{equation*}$

🔗

and make use of the fact that

$\begin{equation*} \frac{y(h) - y(0)}{h} \approx y'(0, y(0)) \qquad \text{or} \qquad y(h) \approx y(0) + y'(0, y(0)), \end{equation*}$

🔗

the estimate for our solution at

$t = 0.1$ is

$\begin{align*} Y_1 & = Y_0 + h f(t_0, Y_0)\\ & = Y_0 + h[Y_0 + t_0]\\ & = 1 + (0.1) [1 + 0]\\ & = 1.1000. \end{align*}$

🔗

Similarly, the approximation at

$t = 0.2$ will be

$\begin{align*} Y_2 & = Y_1 + h f(t_1, Y_1)\\ & = Y_1 + h[Y_1 + t_1]\\ & = 1.1000 + (0.1) [1.1000 + 0.1]\\ & = 1.2200. \end{align*}$

🔗

Our general algorithm is

$\begin{equation*} Y_{k+1} = Y_k + h f(t_k, Y_k) = Y_k + h[Y_k + t_k] = (1.1) Y_k + (0.01)k. \end{equation*}$

🔗

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).

🔗

Table 1.4.1. Euler's approximation for

$y' = y + t$

$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%

three curves that increase at nearly the same rate with the exact race increasing the most quickly and the smaller step sizes increasing faster than the larger step sizes — 🔗
Figure 1.4.2.

🔗

Table 1.4.3. Step sizes for Euler's approximation

$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

$\begin{align*} y' \amp = x + xy\\ y(0) \amp = 1. \end{align*}$

🔗

(a)

Use separation of varibles 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.b) 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

$\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{.}$

🔗

Given the initial value problem

$\begin{align*} y' & = f(t, y),\\ y_0 & = y (t_0), \end{align*}$

🔗

choose

$t_1, t_2, \ldots, t_N$ to be equally spaced on

$[t_0, a]\text{,}$ we can write

$\begin{equation*} t_k = t_0 + kh, \end{equation*}$

🔗

where

$h = (a - t_0)/N$ and

$k = 1, 2, \ldots, N\text{.}$ Taylor's Theorem tells us that

$\begin{equation*} y(t_{k + 1}) = y(t_k + h) = y(t_k) + y'(t_k) h + \frac{y''(t_k)}{2!} h^2 + \cdots. \end{equation*}$

🔗

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,

$\begin{equation*} y_{k + 1} = y(t_{k+1}) \approx y(t_k) + y'(t_k) h = y(t_k) + f(t_k, y_k) h. \end{equation*}$

🔗

This gives us Euler's method,

$\begin{align*} Y_0 & = y(t_0)\\ Y_1 & = Y_0 + h f(t_0, Y_0)\\ Y_2 & = Y_1 + h f(t_1, Y_1)\\ & \vdots\\ Y_{k+1} & = Y_k + h f(t_k, Y_k). \end{align*}$

🔗

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

$\begin{align*} y' & = f(t, y),\\ y(a) & = \alpha, \end{align*}$

where $t \in [a, b]\text{.}$ Suppose that $f$ is continuous and 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}\text{.}$ Also assume that there exists an $M$ such that

$\begin{equation*} |y''(t)| \leq M \end{equation*}$

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

$\begin{equation*} |y(t_i) - Y_i| \leq \frac{hM}{2L} [ e^{L(t_i - a)} - 1]. \end{equation*}$

🔗

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. 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.

🔗

Table 1.4.6. Error bound and actual error

$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

$\begin{equation*} |f(t, y_1) - f(t, y_2)| = |(y_1 + t) - (y_2 + t)| = |y_1 - y_2|, \end{equation*}$

🔗

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

$\begin{equation*} |y(t_i) - Y_i| \leq \frac{hM}{2L} [ e^{L(t_i - a)} - 1] = 0.1e( e^{t_i } - 1) \end{equation*}$

🔗

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,

$\begin{equation*} y(t_1) \approx y_0 + f(t_0, y_0) h + \frac{y''(t_0)}{2} h^2. \end{equation*}$

🔗

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

$\begin{equation*} y'' = \frac{\partial f}{\partial t} \frac{dt}{dt} + \frac{\partial f}{\partial y} \frac{dy}{dt} = f_t + f f_y. \end{equation*}$

🔗

Thus, our approximation becomes

$\begin{equation*} y(t_1) \approx y_0 + f(t_0, y_0) h + \frac{1}{2} \left( f_t(t_0, y_0) + f(t_0, y_0) f_y (t_0, y_0) \right) h^2. \end{equation*}$

🔗

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 Runge-Kutta methods, estimate the higher-order 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

$\begin{align*} y' & = f(t, y),\\ y(t_0) & = y_0, \end{align*}$

🔗

then

$\begin{equation*} y(t_1) = y(t_0) + \int_{t_0}^{t_1} f(s, y(s)) \, ds \end{equation*}$

🔗

$\begin{equation} y_1 - y_0 = y(t_1) - y(t_0) = \int_{t_0}^{t_1} f(s, y(s)) \, ds\label{firstlook04-equation-runge-kutta-1}\tag{1.4.3} \end{equation}$

🔗

by the Fundamental Theorem of Calculus. In Euler's method, we approximate the right-hand side of (1.4.3) by

$\begin{equation*} y_1 - y_0 = f(t_0, y_0) h. \end{equation*}$

🔗

In terms of the definite integral, this is simply a left-hand sum. In the improved Euler's method or the second-order Runge-Kutta method we will estimate the right-hand side of (1.4.3) using the trapezoid rule from calculus,

$\begin{align*} y(t_1) - y(t_0) & = \int_{t_0}^{t_1} f(s, y(s)) \, ds\\ & \approx \left( f(t_0, y_0) + f(t_1, y_1) \right) \frac{h}{2}. \end{align*}$

🔗

Thus, our algorithm becomes

$\begin{equation} y_1 = y_0 + \left( f(t_0, y_0) + f(t_1, y_1) \right) \frac{h}{2}.\label{firstlook04-equation-runge-kutta-2}\tag{1.4.4} \end{equation}$

🔗

However, we have a problem since

$y_1$ appears in the right-hand side of our approximation. To get around this difficulty, we will replace

$y_1$ in the right-hand side of (1.4.4) with the Euler approximation for

$y_1\text{.}$ Thus,

$\begin{equation*} y_1 = y_0 + \left( f(t_0, y_0) + f(t_1, y_0 + f(t_0, y_0) h) \right) \frac{h}{2}. \end{equation*}$

🔗

To understand that the second-order Runge-Kutta 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

$\begin{align*} f(x + h, y + k) & = f(x, y) + h \frac{\partial}{\partial x} f(x, y) + k \frac{\partial}{\partial y}f(x, y)\\ & + \frac{1}{2!} \left( h^2 \frac{\partial^2}{\partial^2 x} f(x, y) + hk \frac{\partial^2}{\partial x \partial y}f(x, y) + k^2 \frac{\partial^2}{\partial^2 y}f(x, y) \right)\\ & + \frac{1}{3!} \left( h^3 \frac{\partial^3}{\partial^3 x} f(x, y) + 3h^2k \frac{\partial^3}{\partial^2 x \partial y}f(x, y) \right.\\ & + \left. hk^2 \frac{\partial^3}{\partial x\partial^2 y}f(x, y) + k^3 \frac{\partial^3}{\partial^3 y} f(x, y) \right) + \cdots. \end{align*}$

🔗

As in the case of the single variable Taylor series, we can write a Taylor polynomial if the Taylor series is truncated,

$\begin{align*} f(x + h, y + k) & = \sum_{n = 0}^{N} \frac{1}{n!} \left( h \frac{\partial}{\partial x} + k \frac{\partial}{\partial y} \right)^n f(x, y)\\ & + \frac{1}{(N+1)!} \left( h \frac{\partial}{\partial x} + k \frac{\partial}{\partial y} \right)^{N+1} f(\overline{x}, \overline{y} ), \end{align*}$

🔗

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

$\begin{equation*} y(t + h) = y(t) + w_1 F_1 + w_2 F_2, \end{equation*}$

🔗

where

$\begin{align*} F_1 & = h f(t, y)\\ F_2 & = h f(t + \alpha h, y + \beta F_1). \end{align*}$

🔗

That is,

$\begin{equation} y(t + h) = y(t) + w_1 h f(t, y) + w_2 h f(t + \alpha h, y + \beta h f(t, y)).\label{firstlook04-equation-runge-kutta-3}\tag{1.4.5} \end{equation}$

🔗

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

$\begin{equation} y(t + h) = y(t) + h y'(t) + \frac{h^2}{2!} y''(t) + \frac{h^3}{3!} y''(t) + \cdots.\label{firstlook04-equation-runge-kutta-4}\tag{1.4.6} \end{equation}$

🔗

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

$\begin{equation*} f(t + \alpha h, y + \beta h f) = f + \alpha h f_t + \beta h f f_y + {\mathcal O}(h^2), \end{equation*}$

🔗

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),

$\begin{equation} y(t + h) = y(t) + (w_1 + w_2) hf + \alpha w_2 h^2 f_t + \beta w_2 h^2 f f_y + {\mathcal O}(h^3).\label{firstlook04-equation-runge-kutta-5}\tag{1.4.7} \end{equation}$

🔗

Since

$y'' = f_t + f_y f$ by the chain rule, we can rewrite (1.4.6) as

$\begin{equation} y(t + h) = y(t) + h f + \frac{h^2}{2!}( f_t + f f_y) + {\mathcal O}(h^3).\label{firstlook04-equation-runge-kutta-6}\tag{1.4.8} \end{equation}$

🔗

We can make equations (1.4.7) and (1.4.8) agree up through the quadratic terms if we require that

$\begin{align*} w_1 + w_2 & = 1,\\ \alpha w_2 & = \frac{1}{2},\\ \beta w_2 & = \frac{1}{2}. \end{align*}$

🔗

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

$\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*}$

🔗

The improved Euler's method or the second-order Runge-Kutta 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

$\begin{equation*} y(t + h) = y(t) + h y'(t) + {\mathcal O}(h^2), \end{equation*}$

🔗

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

$\begin{equation*} y(t + h) = y(t) + h y'(t) + \frac{h^2}{2!} y''(t) + {\mathcal O}(h^3). \end{equation*}$

🔗

If we use Simpson's rule to estimate the integral in

$\begin{equation*} y(t_1) - y(t_0) = \int_{t_0}^{t_1} f(s, y(s)) \, ds \end{equation*}$

🔗

we can improve our accuracy up to

$h^4\text{.}$ The idea is exactly the same, but the algebra become much more tedious. This method is known as the Runge-Kutta method of order 4 and 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*}$

🔗

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 = (b-a)/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_{N-1} + h f(t_{N-1}, Y_{N-1}). \end{align*}$
In practice, no one uses Euler's method. The Runge-Kutta methods are better algorithms.
Taylor's Theorem 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 Runge-Kutta 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 Runge-Kutta 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{.}$

🔗

Exercises 1.4.5 Exercises

🔗

Finding Solutions

For each of the initial value problem

$\begin{align*} y' & = f(t, y),\\ y(t_0) & = y_0 \end{align*}$

in Exercise Group 1.4.5.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$

Hint

This equation is a first-order linear equation (Section 1.5), but it is possible to find the analytic solution using Sage (Subsection 1.2.9).

🔗

6.

$y' = 1/y\text{,}$ $y'(0) = 2$

🔗

7.

Consider the differential equation

$\begin{equation*} \frac{dy}{dt} = 3y - 1 \end{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 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, N-1\text{.}$

Hint

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.6 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 Runge-Kutta 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.

xxxxxxxxxx
 
y = function('y')(x)
de = diff(y,x) == x + y
solution = desolve(de, y, ics = [0,1])
solution.show()
plot(solution, x, 0, 5)

🔗

We can also use Euler's method to find a solution for our initial value problem.

xxxxxxxxxx
 
x,y = PolynomialRing(RR, 2, "xy").gens() #declare x, y as generators of a polynomial ring
eulers_method(x + y, 0,1, 0.5, 5)

🔗

The syntax of eulers_method for the inital value problem

$\begin{align*} y \amp = f(x,y)\\ f(x_0) \amp = y_0 \end{align*}$

🔗

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{.}$

xxxxxxxxxx
 
x,y = PolynomialRing(RR, 2, "xy").gens()
pts = eulers_method(x + y, 0, 1, 1/2, 4.5, algorithm="none")
p = list_plot(pts, color="red")
p += line(pts)
p

🔗

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 Runge-Kutta method.

xxxxxxxxxx
 
x,y = var('x y')
desolve_rk4(x + y, y, ics=[0, 1], end_points = 5, step = 0.5)

🔗

Again, we just get a list of points. However, desolve_rk4 has some nice built-in graphing utilities.

xxxxxxxxxx
 
x,y = var('x y')
p = desolve_rk4(x + y, y, ics = [0,1], ivar = x, output = 'slope_field', end_points = [0, 5], thickness = 2, color = 'red')
p

🔗

Write the Sage commands to compare the graphs obtained using eulers_method and desolve with the exact solution.

xxxxxxxxxx
 
​

🔗

Not only is desolve_rk4 more accurate, it is much easier to use. For more information, see http://www.sagemath.org/doc/reference/calculus/sage/calculus/desolvers.html.