3.1 Euler’s Method

If an initial value problem

\begin{equation} y'=f(x,y),\quad y(x_0)=y_0 \tag{3.1.1}\end{equation}

can’t be solved analytically, it’s necessary to resort to numerical methods to obtain useful approximations to a solution of (3.1.1). We’ll consider such methods in this chapter.

We’re interested in computing approximate values of the solution of (3.1.1) at equally spaced points \(x_0\), \(x_1\), …, \(x_n=b\) in an interval \([x_0,b]\). Thus,

\[ x_i=x_0+ih,\quad i=0,1, \dots,n, \]

where

\[ h={b-x_0\over n}. \]

We’ll denote the approximate values of the solution at these points by \(y_0\), \(y_1\), …, \(y_n\); thus, \(y_i\) is an approximation to \(y(x_i)\). We’ll call

\[ e_i=y(x_i)-y_i \]

the error at the \(i\)th step. Because of the initial condition \(y(x_0)=y_0\), we’ll always have \(e_0=0\). However, in general \(e_i\ne0\) if \(i>0\).

We encounter two sources of error in applying a numerical method to solve an initial value problem:

Since a careful analysis of roundoff error is beyond the scope of this book, we’ll consider only truncation errors.

Euler’s Method

The simplest numerical method for solving (3.1.1) is Euler’s method. This method is so crude that it is seldom used in practice; however, its simplicity makes it useful for illustrative purposes.

Euler’s method is based on the assumption that the tangent line to the integral curve of (3.1.1) at \((x_i,y(x_i))\) approximates the integral curve over the interval \([x_i,x_{i+1}]\). Since the slope of the integral curve of (3.1.1) at \((x_i,y(x_i))\) is \(y'(x_i)=f(x_i,y(x_i))\), the equation of the tangent line to the integral curve at \((x_i,y(x_i))\) is

\begin{equation} y=y(x_i)+f(x_i,y(x_i))(x-x_i). \tag{3.1.2}\end{equation}

Setting \(x=x_{i+1}=x_i+h\) in (3.1.2) yields

\begin{equation} y_{i+1}=y(x_i)+hf(x_i,y(x_i)) \tag{3.1.3}\end{equation}

as an approximation to \(y(x_{i+1})\). Since \(y(x_0)=y_0\) is known, we can use (3.1.3) with \(i=0\) to compute

\[ y_1=y_0+hf(x_0,y_0). \]

However, setting \(i=1\) in (3.1.3) yields

\[ y_2=y(x_1)+hf(x_1,y(x_1)), \]

which isn’t useful, since we don’t know \(y(x_1)\). Therefore we replace \(y(x_1)\) by its approximate value \(y_1\) and redefine

\[ y_2=y_1+hf(x_1,y_1). \]

Having computed \(y_2\), we can compute

\[ y_3=y_2+hf(x_2,y_2). \]

In general, Euler’s method starts with the known value \(y(x_0)=y_0\) and computes \(y_1\), \(y_2\), …, \(y_n\) successively by with the formula

\begin{equation} y_{i+1}=y_i+hf(x_i,y_i),\quad 0\le i\le n-1. \tag{3.1.4}\end{equation}

The next example illustrates the computational procedure indicated in Euler’s method.

Example 3.1.1

Use Euler’s method with \(h=0.1\) to find approximate values for the solution of the initial value problem

\begin{equation} y'+2y=x^3e^{-2x},\quad y(0)=1 \tag{3.1.5}\end{equation}

at \(x=0.1,0.2,0.3\).

Solution We rewrite (3.1.5) as

\[ y'=-2y+x^3e^{-2x},\quad y(0)=1, \]

which is of the form (3.1.1), with

\[ f(x,y)=-2y+x^3e^{-2x},\ x_0=0,\mbox{\ and}\ y_0=1. \]

Euler’s method yields

\begin{eqnarray*} y_1 & = & y_0+hf(x_0,y_0) \\ & = & 1+(.1)f(0,1)=1+(.1)(-2)=.8, \\[12pt] y_2 & = & y_1+hf(x_1,y_1) \\ & = & .8+(.1)f(.1,.8)=.8+(.1)\left(-2(.8)+(.1)^3e^{-.2}\right)= .640081873, \\[12pt] y_3 & = & y_2+hf(x_2,y_2) \\ & = & .640081873+(.1)\left(-2(.640081873)+(.2)^3e^{-.4}\right)= .512601754.\bbox \end{eqnarray*}

We’ve written the details of these computations to ensure that you understand the procedure. However, in the rest of the examples as well as the exercises in this chapter, we’ll assume that you can use a programmable calculator or a computer to carry out the necessary computations.

Examples Illustrating The Error in Euler’s Method

Example 3.1.2

Use Euler’s method with step sizes \(h=0.1\), \(h=0.05\), and \(h=0.025\) to find approximate values of the solution of the initial value problem

\[ y'+2y=x^3e^{-2x},\quad y(0)=1 \]

at \(x=0\), \(0.1\), \(0.2\), \(0.3\), …, \(1.0\). Compare these approximate values with the values of the exact solution

\begin{equation} y={e^{-2x}\over4}(x^4+4), \tag{3.1.6}\end{equation}

which can be obtained by the method of Section 2.1. (Verify.)

Solution Table 3.1.1 shows the values of the exact solution (3.1.6) at the specified points, and the approximate values of the solution at these points obtained by Euler’s method with step sizes \(h=0.1\), \(h=0.05\), and \(h=0.025\). In examining this table, keep in mind that the approximate values in the column corresponding to \(h=.05\) are actually the results of 20 steps with Euler’s method. We haven’t listed the estimates of the solution obtained for \(x=0.05\), \(0.15\), …, since there’s nothing to compare them with in the column corresponding to \(h=0.1\). Similarly, the approximate values in the column corresponding to \(h=0.025\) are actually the results of 40 steps with Euler’s method.

Table 3.1.1. Numerical solution of \(y'+2y=x^3e^{-2x},\ y(0)=1\), by Euler’s method.
\(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)Exact
0.01.0000000001.0000000001.0000000001.000000000
0.10.8000000000.8100056550.8145183490.818751221
0.20.6400818730.6562664370.6636359530.670588174
0.30.5126017540.5322909810.5413394950.549922980
0.40.4115631950.4328870560.4427747660.452204669
0.50.3321262610.3537850150.3639155970.373627557
0.60.2702995020.2914042560.3013598850.310952904
0.70.2227453970.2427072570.2522029350.261398947
0.80.1866545930.2051057540.2139563110.222570721
0.90.1596607760.1763968830.1844924630.192412038
1.00.1397789100.1547159250.1620032930.169169104

You can see from Table 3.1.1 that decreasing the step size improves the accuracy of Euler’s method. For example,

\[ y_{\mbox{\tiny exact}}(1)-y_{\mbox{\tiny approx}}(1)\approx \left\{\begin{array}{l} .0293\mbox{ with }h=0.1,\\ .0144\mbox{ with }h=0.05,\\ .0071\mbox{ with }h=0.025. \end{array}\right. \]

Based on this scanty evidence, you might guess that the error in approximating the exact solution at a fixed value of \(x\) by Euler’s method is roughly halved when the step size is halved. You can find more evidence to support this conjecture by examining Table 3.1.2, which lists the approximate values of \(y_{\mbox{\tiny exact}}-y_{\mbox{\tiny approx}}\) at \(x=0.1\), \(0.2\), …, \(1.0\).

Table 3.1.2. Errors in approximate solutions of \(y'+2y=x^3e^{-2x},\ y(0)=1\), obtained by Euler’s method.
\(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)
0.10.01870.00870.0042
0.20.03050.01430.0069
0.30.03730.01760.0085
0.40.04060.01930.0094
0.50.04150.01980.0097
0.60.04060.01950.0095
0.70.03860.01860.0091
0.80.03590.01740.0086
0.90.03270.01600.0079
1.00.02930.01440.0071

Example 3.1.3

Tables 3.1.3 and 3.1.4 show analogous results for the nonlinear initial value problem

\begin{equation} y'=-2y^2+xy+x^2,\ y(0)=1, \tag{3.1.7}\end{equation}

except in this case we can’t solve (3.1.7) exactly. The results in the “Exact” column were obtained by using a more accurate numerical method known as the Runge-Kutta method with a small step size. They are exact to eight decimal places.

Since we think it’s important in evaluating the accuracy of the numerical methods that we’ll be studying in this chapter, we often include a column listing values of the exact solution of the initial value problem, even if the directions in the example or exercise don’t specifically call for it. If quotation marks are included in the heading, the values were obtained by applying the Runge-Kutta method in a way that’s explained in Section 3.3. If quotation marks are not included, the values were obtained from a known formula for the solution. In either case, the values are exact to eight places to the right of the decimal point.

Table 3.1.3. Numerical solution of \(y'=-2y^2+xy+x^2,\ y(0)=1\), by Euler’s method.
\(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)“Exact”
0.01.0000000001.0000000001.0000000001.000000000
0.10.8000000000.8213750000.8299770070.837584494
0.20.6810000000.7077953770.7192262530.729641890
0.30.6058678000.6337765900.6461152270.657580377
0.40.5596286760.5874545260.6000457010.611901791
0.50.5353769720.5629061690.5755563910.587575491
0.60.5298201200.5571435350.5698241710.581942225
0.70.5414674550.5687169350.5814354230.593629526
0.80.5697327760.5969519880.6096849030.621907458
0.90.6143923110.6414577290.6541108620.666250842
1.00.6751920370.7017644950.7141516260.726015790
Table 3.1.4. Errors in approximate solutions of \(y'=-2y^2+xy+x^2,\ y(0)=1\), obtained by Euler’s method.
\(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)
0.10.03760.01620.0076
0.20.04860.02180.0104
0.30.05170.02380.0115
0.40.05230.02440.0119
0.50.05220.02470.0121
0.60.05210.02480.0121
0.70.05220.02490.0122
0.80.05220.02500.0122
0.90.05190.02480.0121
1.00.05080.02430.0119

Truncation Error in Euler’s Method

Consistent with the results indicated in Tables 3.1.13.1.4, we’ll now show that under reasonable assumptions on \(f\) there’s a constant \(K\) such that the error in approximating the solution of the initial value problem

\[ y'=f(x,y),\quad y(x_0)=y_0, \]

at a given point \(b>x_0\) by Euler’s method with step size \(h=(b-x_0)/n\) satisfies the inequality

\[ |y(b)-y_n|\le Kh, \]

where \(K\) is a constant independent of \(n\).

There are two sources of error (not counting roundoff) in Euler’s method:

  1. The error committed in approximating the integral curve by the tangent line (3.1.2) over the interval \([x_i,x_{i+1}]\).

  2. The error committed in replacing \(y(x_i)\) by \(y_i\) in (3.1.2) and using (3.1.4) rather than (3.1.2) to compute \(y_{i+1}\).

Euler’s method assumes that \(y_{i+1}\) defined in (3.1.2) is an approximation to \(y(x_{i+1})\). We call the error in this approximation the local truncation error at the \(i\)th step, and denote it by \(T_i\); thus,

\begin{equation} T_i=y(x_{i+1})-y(x_i)-hf(x_i,y(x_i)). \tag{3.1.8}\end{equation}

We’ll now use Taylor’s theorem to estimate \(T_i\), assuming for simplicity that \(f\), \(f_x\), and \(f_y\) are continuous and bounded for all \((x,y)\). Then \(y''\) exists and is bounded on \([x_0,b]\). To see this, we differentiate

\[ y'(x)=f(x,y(x)) \]

to obtain

\begin{eqnarray*} y''(x) & = & f_x(x,y(x))+f_y(x,y(x))y'(x) \\ & = & f_x(x,y(x))+f_y(x,y(x))f(x,y(x)). \end{eqnarray*}

Since we assumed that \(f\), \(f_x\) and \(f_y\) are bounded, there’s a constant \(M\) such that

\[ |f_x(x,y(x))+f_y(x,y(x))f(x,y(x))|\le M,\quad x_0<x<b, \]

which implies that

\begin{equation} |y''(x)|\le M,\quad x_0<x<b. \tag{3.1.9}\end{equation}

Since \(x_{i+1}=x_i+h\), Taylor’s theorem implies that

\[ y(x_{i+1})=y(x_i)+hy'(x_i)+{h^2\over2}y''(\tilde x_i), \]

where \(\tilde x_i\) is some number between \(x_i\) and \(x_{i+1}\). Since \(y'(x_i)=f(x_i,y(x_i))\) this can be written as

\[ y(x_{i+1})=y(x_i)+hf(x_i,y(x_i))+{h^2\over2}y''(\tilde x_i), \]

or, equivalently,

\[ y(x_{i+1})-y(x_i)-hf(x_i,y(x_i))={h^2\over2}y''(\tilde x_i). \]

Comparing this with (3.1.8) shows that

\[ T_i={h^2\over2}y''(\tilde x_i). \]

Recalling (3.1.9), we can establish the bound

\begin{equation} |T_i|\le{Mh^2\over2},\quad 1\le i\le n. \tag{3.1.10}\end{equation}

Although it may be difficult to determine the constant \(M\), what is important is that there’s an \(M\) such that (3.1.10) holds. We say that the local truncation error of Euler’s method is of order \(h^2\), which we write as \(O(h^2)\).

Note that the magnitude of the local truncation error in Euler’s method is determined by the second derivative \(y''\) of the solution of the initial value problem. Therefore the local truncation error will be larger where \(|y''|\) is large, or smaller where \(|y''|\) is small.

Since the local truncation error for Euler’s method is \(O(h^2)\), it’s reasonable to expect that halving \(h\) reduces the local truncation error by a factor of 4. This is true, but halving the step size also requires twice as many steps to approximate the solution at a given point. To analyze the overall effect of truncation error in Euler’s method, it’s useful to derive an equation relating the errors

\[ e_{i+1}=y(x_{i+1})-y_{i+1}\mbox{\quad and \quad} e_i=y(x_i)-y_i. \]

To this end, recall that

\begin{equation} y(x_{i+1})=y(x_i)+hf(x_i,y(x_i))+T_i \tag{3.1.11}\end{equation}

and

\begin{equation} y_{i+1}=y_i+hf(x_i,y_i). \tag{3.1.12}\end{equation}

Subtracting (3.1.12) from (3.1.11) yields

\begin{equation} e_{i+1}=e_i+h\left[f(x_i,y(x_i))-f(x_i,y_i)\right]+T_i. \tag{3.1.13}\end{equation}

The last term on the right is the local truncation error at the \(i\)th step. The other terms reflect the way errors made at previous steps affect \(e_{i+1}\). Since \(|T_i|\le Mh^2/2\), we see from (3.1.13) that

\begin{equation} |e_{i+1}|\le |e_i|+h|f(x_i,y(x_i))-f(x_i,y_i)|+{Mh^2\over2}. \tag{3.1.14}\end{equation}

Since we assumed that \(f_y\) is continuous and bounded, the mean value theorem implies that

\[ f(x_i,y(x_i))-f(x_i,y_i)=f_y(x_i,y_i^*)(y(x_i)-y_i)=f_y(x_i,y_i^*)e_i, \]

where \(y_i^*\) is between \(y_i\) and \(y(x_i)\). Therefore

\[ |f(x_i,y(x_i))-f(x_i,y_i)|\le R|e_i| \]

for some constant \(R\). From this and (3.1.14),

\begin{equation} |e_{i+1}|\le (1+Rh)|e_i|+{Mh^2\over2},\quad 0\le i\le n-1. \tag{3.1.15}\end{equation}

For convenience, let \(C=1+Rh\). Since \(e_0=y(x_0)-y_0=0\), applying (3.1.15) repeatedly yields

\begin{eqnarray} |e_1| & \le & {Mh^2\over2} \\ |e_2| & \le & C|e_1|+{Mh^2\over2}\le(1+C){Mh^2\over2} \\ |e_3| & \le & C|e_2|+{Mh^2\over2}\le(1+C+C^2){Mh^2\over2} \\ & \vdots & \\ |e_n| & \le & C|e_{n-1}|+{Mh^2\over2}\le(1+C+\cdots+C^{n-1}){Mh^2\over2}. \tag{3.1.16}\end{eqnarray}

Recalling the formula for the sum of a geometric series, we see that

\[ 1+C+\cdots+C^{n-1}={1-C^n\over 1-C}={(1+Rh)^n-1\over Rh} \]

(since \(C=1+Rh\)). From this and (3.1.16),

\begin{equation} |y(b)-y_n|=|e_n|\le{(1+Rh)^n-1\over R}{Mh\over2}. \tag{3.1.17}\end{equation}

Since Taylor’s theorem implies that

\[ 1+Rh<e^{Rh} \]

(verify),

\[ (1+Rh)^n<e^{nRh}=e^{R(b-x_0)}\quad (\mbox{since }nh=b-x_0). \]

This and (3.1.17) imply that

\begin{equation} |y(b)-y_n|\le Kh, \tag{3.1.18}\end{equation}

with

\[ K=M{e^{R(b-x_0)}-1\over2R}. \]

Because of (3.1.18) we say that the global truncation error of Euler’s method is of order \(h\), which we write as \(O(h)\).

Semilinear Equations and Variation of Parameters

An equation that can be written in the form

\begin{equation} y'+p(x)y=h(x,y) \tag{3.1.19}\end{equation}

with \(p\not\equiv0\) is said to be semilinear. (Of course, (3.1.19) is linear if \(h\) is independent of \(y\).) One way to apply Euler’s method to an initial value problem

\begin{equation} y'+p(x)y=h(x,y),\quad y(x_0)=y_0 \tag{3.1.20}\end{equation}

for (3.1.19) is to think of it as

\[ y'=f(x,y),\quad y(x_0)=y_0, \]

where

\[ f(x,y)=-p(x)y+h(x,y). \]

However, we can also start by applying variation of parameters to (3.1.20), as in Sections 2.1 and 2.4; thus, we write the solution of (3.1.20) as \(y=uy_1\), where \(y_1\) is a nontrivial solution of the complementary equation \(y'+p(x)y=0\). Then \(y=uy_1\) is a solution of (3.1.20) if and only if \(u\) is a solution of the initial value problem

\begin{equation} u'=h(x,uy_1(x))/y_1(x),\quad u(x_0)=y(x_0)/y_1(x_0). \tag{3.1.21}\end{equation}

We can apply Euler’s method to obtain approximate values \(u_0\), \(u_1\), …, \(u_n\) of this initial value problem, and then take

\[ y_i=u_iy_1(x_i) \]

as approximate values of the solution of (3.1.20). We’ll call this procedure the Euler semilinear method.

The next two examples show that the Euler and Euler semilinear methods may yield drastically different results.

Example 3.1.4

In Example 2.1.7 we had to leave the solution of the initial value problem

\begin{equation} y'-2xy=1,\quad y(0)=3 \tag{3.1.22}\end{equation}

in the form

\begin{equation} y=e^{x^2}\left(3 +\int^x_0 e^{-t^2}dt\right) \tag{3.1.23}\end{equation}

because it was impossible to evaluate this integral exactly in terms of elementary functions. Use step sizes \(h=0.2\), \(h=0.1\), and \(h=0.05\) to find approximate values of the solution of (3.1.22) at \(x=0\), \(0.2\), \(0.4\), \(0.6\), …, \(2.0\) by (a) Euler’s method; (b) the Euler semilinear method.

Solution (a) Rewriting (3.1.22) as

\begin{equation} y'=1+2xy,\quad y(0)=3 \tag{3.1.24}\end{equation}

and applying Euler’s method with \(f(x,y)=1+2xy\) yields the results shown in Table 3.1.5. Because of the large differences between the estimates obtained for the three values of \(h\), it would be clear that these results are useless even if the “exact” values were not included in the table.

Table 3.1.5. Numerical solution of \(y'-2xy=1,\ y(0)=3\), with Euler’s method.
\(x\)\(h=0.2\)\(h=0.1\)\(h=0.05\)“Exact”
0.03.0000000003.0000000003.0000000003.000000000
0.23.2000000003.2620000003.2943485373.327851973
0.43.6560000003.8020288003.8814211033.966059348
0.64.4409600004.7268102144.8888707835.067039535
0.85.7067904006.2491912826.5707962356.936700945
1.07.7329633288.7718930269.41910562010.184923955
1.211.02614865913.06405139114.40577206716.067111677
1.416.51870001620.63727389323.52293587227.289392347
1.625.96917202434.57042375841.03344125750.000377775
1.842.78944212061.38216554376.49101824698.982969504
2.073.797840446115.440048291152.363866569211.954462214

It’s easy to see why Euler’s method yields such poor results. Recall that the constant \(M\) in (3.1.10) – which plays an important role in determining the local truncation error in Euler’s method – must be an upper bound for the values of the second derivative \(y''\) of the solution of the initial value problem (3.1.22) on \((0,2)\). The problem is that \(y''\) assumes very large values on this interval. To see this, we differentiate (3.1.24) to obtain

\[ y''(x)=2y(x)+2xy'(x)=2y(x)+2x(1+2xy(x))=2(1+2x^2)y(x)+2x, \]

where the second equality follows again from (3.1.24). Since (3.1.23) implies that \(y(x)>3e^{x^2}\) if \(x>0\),

\[ y''(x)>6(1+2x^2)e^{x^2}+2x,\quad x>0. \]

For example, letting \(x=2\) shows that \(y''(2)>2952\).

Solution (b) Since \(y_1=e^{x^2}\) is a solution of the complementary equation \(y'-2xy=0\), we can apply the Euler semilinear method to (3.1.22), with

\[ y=ue^{x^2}\mbox{\quad and \quad} u'=e^{-x^2},\quad u(0)=3. \]

The results listed in Table 3.1.6 are clearly better than those obtained by Euler’s method.

Table 3.1.6. Numerical solution of \(y'-2xy=1,\ y(0)=3\), by the Euler semilinear method.
\(x\)\(h=0.2\)\(h=0.1\)\(h=0.05\)“Exact”
0.03.0000000003.0000000003.0000000003.000000000
0.23.3305944773.3295588533.3287888893.327851973
0.43.9807341573.9740676283.9702304153.966059348
0.65.1063602315.0877052445.0776227235.067039535
0.87.0210034176.9801908916.9587795866.936700945
1.010.35007660010.26917082410.22746429910.184923955
1.216.38118009216.22614639016.14712906716.067111677
1.427.89000338027.59202608527.44129223527.289392347
1.651.18332326250.59450386350.29810665950.000377775
1.8101.424397595100.20665907699.59556276698.982969504
2.0217.301032800214.631041938213.293582978211.954462214

We can’t give a general procedure for determining in advance whether Euler’s method or the semilinear Euler method will produce better results for a given semilinear initial value problem (3.1.19). As a rule of thumb, the Euler semilinear method will yield better results than Euler’s method if \(|u''|\) is small on \([x_0,b]\), while Euler’s method yields better results if \(|u''|\) is large on \([x_0,b]\). In many cases the results obtained by the two methods don’t differ appreciably. However, we propose the an intuitive way to decide which is the better method: Try both methods with multiple step sizes, as we did in Example 3.1.4, and accept the results obtained by the method for which the approximations change less as the step size decreases.

Example 3.1.5

Applying Euler’s method with step sizes \(h=0.1\), \(h=0.05\), and \(h=0.025\) to the initial value problem

\begin{equation} y'-2y={x\over1+y^2},\quad y(1)=7 \tag{3.1.25}\end{equation}

on \([1,2]\) yields the results in Table 3.1.7. Applying the Euler semilinear method with

\[ y=ue^{2x}\mbox{\quad and \quad}u'={xe^{-2x}\over1+u^2e^{4x}},\quad u(1)=7e^{-2} \]

yields the results in Table 3.1.8. Since the latter are clearly less dependent on step size than the former, we conclude that the Euler semilinear method is better than Euler’s method for (3.1.25). This conclusion is supported by comparing the approximate results obtained by the two methods with the “exact” values of the solution.

Table 3.1.7. Numerical solution of \(y'-2y=x/(1+y^2),\ y(1)=7\), by Euler’s method.
\(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)“Exact”
1.07.0000000007.0000000007.0000000007.000000000
1.18.4020000008.4719705698.5104939558.551744786
1.210.08393645010.25257016910.34601410110.446546230
1.312.10189235412.40671938112.57672082712.760480158
1.414.52315244515.01295241615.28787210415.586440425
1.517.42844355418.16627740518.58307940619.037865752
1.620.91462447121.98163848722.58826621723.253292359
1.725.09791431026.59810518027.45647969528.401914416
1.830.11776662732.18394134033.37373894434.690375086
1.936.14151817238.94273825240.56614315842.371060528
2.043.36996715547.12083525149.30851112651.752229656
Table 3.1.8. Numerical solution of \(y'-2y=x/(1+y^2),\ y(1)=7\), by the Euler semilinear method.
\(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)“Exact”
1.07.0000000007.0000000007.0000000007.000000000
1.18.5522621138.5519939788.5518670078.551744786
1.210.44756867410.44703854710.44678764610.446546230
1.312.76201979912.76122131312.76084354312.760480158
1.415.58853514115.58744860015.58693468015.586440425
1.519.04058061419.03917224119.03850621119.037865752
1.623.25672163623.25494251723.25410125323.253292359
1.728.40618459728.40396910728.40292158128.401914416
1.834.69564922234.69291276834.69161897934.690375086
1.942.37754413842.37418009042.37258962442.371060528
2.051.76017844651.75605413351.75410426251.752229656

Example 3.1.6

Applying Euler’s method with step sizes \(h=0.1\), \(h=0.05\), and \(h=0.025\) to the initial value problem

\begin{equation} y'+3x^2y=1+y^2,\quad y(2)=2 \tag{3.1.26}\end{equation}

on \([2,3]\) yields the results in Table 3.1.9. Applying the Euler semilinear method with

\[ y=ue^{-x^3}\mbox{\quad and \quad}u'=e^{x^3}(1+u^2e^{-2x^3}),\quad u(2)=2e^8 \]

yields the results in Table 3.1.10. Noting the close agreement among the three columns of Table 3.1.9 (at least for larger values of \(x\)) and the lack of any such agreement among the columns of Table 3.1.10, we conclude that Euler’s method is better than the Euler semilinear method for (3.1.26). Comparing the results with the exact values supports this conclusion.

Table 3.1.9. Numerical solution of \(y'+3x^2y=1+y^2,\quad y(2)=2\), by Euler’s method.
\(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)“Exact”
2.02.0000000002.0000000002.0000000002.000000000
2.10.1000000000.4932312500.6096111710.701162906
2.20.0687000000.1228795860.1801134450.236986800
2.30.0694195690.0706708900.0839344590.103815729
2.40.0597326210.0613389560.0633375610.068390786
2.50.0568714510.0560023630.0562496700.057281091
2.60.0505609170.0514652560.0515175010.051711676
2.70.0482790180.0474847160.0475142020.047564141
2.80.0429258920.0439670020.0439892390.044014438
2.90.0421484580.0408396830.0408571090.040875333
3.00.0359855480.0380446920.0380585360.038072838
Table 3.1.10. Numerical solution of \(y'+3x^2y=1+y^2,\quad y(2)=2\), by the Euler semilinear method.
\(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)“Exact”
\(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)\(h=.0125\)
2.02.0000000002.0000000002.0000000002.000000000
2.10.7084262860.7025681710.7012142740.701162906
2.20.2145018520.2225994680.2289422400.236986800
2.30.0698614360.0836204940.0928528060.103815729
2.40.0324873960.0470792610.0568258050.068390786
2.50.0218955590.0360300180.0456838010.057281091
2.60.0173320580.0307501810.0401899200.051711676
2.70.0142714920.0269319110.0361346740.047564141
2.80.0118195550.0237206700.0326797670.044014438
2.90.0097767920.0209255220.0296365060.040875333
3.00.0080650200.0184723020.0269310990.038072838

In the next two sections we’ll study other numerical methods for solving initial value problems, called the improved Euler method, the midpoint method, Heun’s method and the Runge-Kutta method. If the initial value problem is semilinear as in (3.1.19), we also have the option of using variation of parameters and then applying the given numerical method to the initial value problem (3.1.21) for \(u\). By analogy with the terminology used here, we’ll call the resulting procedure the improved Euler semilinear method, the midpoint semilinear method, Heun’s semilinear method or the Runge-Kutta semilinear method, as the case may be.

3.1 Exercises

You may want to save the results of these exercises, sincewe’ll revisit in the next two sections. In Exercises 15 use Euler’s method to find approximate values of the solution of the given initial value problem at the points \(x_i=x_0+ih\), where \(x_0\) is the point wher the initial condition is imposed and \(i=1\), \(2\), \(3\). The purpose of these exercises is to familiarize you with the computational procedure of Euler’s method.

  1. C \(y'=2x^2+3y^2-2,\quad y(2)=1;\quad h=0.05\)

    Show answer

    \(y_1=1.450000000,\; y_2=2.085625000,\; y_3= 3.079099746\)

  2. C \(y'=y+\sqrt{x^2+y^2},\quad y(0)=1;\quad h=0.1\)

    Show answer

    \(y_1=1.200000000,\; y_2=1.440415946,\; y_3=1.729880994\)

  3. C \(y'+3y=x^2-3xy+y^2,\quad y(0)=2;\quad h=0.05\)

    Show answer

    \(y_1=1.900000000,\; y_2=1.781375000,\; y_3=1.646612970\)

  4. C \(y'=\dst{1+x\over1-y^2},\quad y(2)=3;\quad h=0.1\)

    Show answer

    \(y_1=2.962500000,\; y_2=2.922635828,\; y_3=2.880205639\)

  5. C \(y'+x^2y=\sin xy,\quad y(1)=\pi;\quad h=0.2\)

    Show answer

    \(y_1=2.513274123,\; y_2=1.814517822,\; y_3=1.216364496\)

  6. C Use Euler’s method with step sizes \(h=0.1\), \(h=0.05\), and \(h=0.025\) to find approximate values of the solution of the initial value problem

    \[ y'+3y=7e^{4x},\quad y(0)=2 \]

    at \(x=0\), \(0.1\), \(0.2\), \(0.3\), …, \(1.0\). Compare these approximate values with the values of the exact solution \(y=e^{4x}+e^{-3x}\), which can be obtained by the method of Section 2.1. Present your results in a table like Table 3.1.1.

    Show answer

    \(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)Exact
    1.048.29814736251.49282564353.07667368554.647937102

  7. C Use Euler’s method with step sizes \(h=0.1\), \(h=0.05\), and \(h=0.025\) to find approximate values of the solution of the initial value problem

    \[ y'+{2\over x}y={3\over x^3}+1,\quad y(1)=1 \]

    at \(x=1.0\), \(1.1\), \(1.2\), \(1.3\), …, \(2.0\). Compare these approximate values with the values of the exact solution

    \[ y={1\over3x^2}(9\ln x+x^3+2), \]

    which can be obtained by the method of Section 2.1. Present your results in a table like Table 3.1.1.

    Show answer

    \(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)Exact
    2.01.3902420091.3709967581.3619211321.353193719

  8. C Use Euler’s method with step sizes \(h=0.05\), \(h=0.025\), and \(h=0.0125\) to find approximate values of the solution of the initial value problem

    \[ y'={y^2+xy-x^2\over x^2},\quad y(1)=2 \]

    at \(x=1.0\), \(1.05\), \(1.10\), \(1.15\), …, \(1.5\). Compare these approximate values with the values of the exact solution

    \[ y={x(1+x^2/3)\over1-x^2/3} \]

    obtained in Example 2.4.3. Present your results in a table like Table 3.1.1.

    Show answer

    \(x\)\(h=0.05\)\(h=0.025\)\(h=0.0125\)Exact
    1.507.8861704378.8524637939.54803990710.500000000

  9. C In Example 2.2.3 it was shown that

    \[ y^5+y=x^2+x-4 \]

    is an implicit solution of the initial value problem

    \[ y'={2x+1\over5y^4+1},\quad y(2)=1. \tag*{\rm(A)} \]

    Use Euler’s method with step sizes \(h=0.1\), \(h=0.05\), and \(h=0.025\) to find approximate values of the solution of (A) at \(x=2.0\), \(2.1\), \(2.2\), \(2.3\), …, \(3.0\). Present your results in tabular form. To check the error in these approximate values, construct another table of values of the residual

    \[ R(x,y)=y^5+y-x^2-x+4 \]

    for each value of \((x,y)\) appearing in the first table.

    Show answer

    \(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)\(h=0.1\)\(h=0.05\)\(h=0.025\)
    3.01.4694582411.4625144861.4592170100.32100.15370.0753
    Approximate SolutionsResiduals

  10. C You can see from Example 2.5.1 that

    \[ x^4y^3+x^2y^5+2xy=4 \]

    is an implicit solution of the initial value problem

    \[ y'=-{4x^3y^3+2xy^5+2y\over3x^4y^2+5x^2y^4+2x},\quad y(1)=1. \tag*{\rm(A)} \]

    Use Euler’s method with step sizes \(h=0.1\), \(h=0.05\), and \(h=0.025\) to find approximate values of the solution of (A) at \(x=1.0\), \(1.1\), \(1.2\), \(1.3\), …, \(2.0\). Present your results in tabular form. To check the error in these approximate values, construct another table of values of the residual

    \[ R(x,y)=x^4y^3+x^2y^5+2xy-4 \]

    for each value of \((x,y)\) appearing in the first table.

    Show answer

    \(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)\(h=0.1\)\(h=0.05\)\(h=0.025\)
    2.00.4734567370.4832274700.487986391-0.3129-0.1563-0.0781
    Approximate SolutionsResiduals

  11. C Use Euler’s method with step sizes \(h=0.1\), \(h=0.05\), and \(h=0.025\) to find approximate values of the solution of the initial value problem

    \[ (3y^2+4y)y'+2x+\cos x=0, \quad y(0)=1; \mbox{\; (Exercise~2.2.~\hspace*{-3pt}\ref{exer:2.2.13})} \]

    at \(x=0\), \(0.1\), \(0.2\), \(0.3\), …, \(1.0\).

    Show answer

    \(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)“Exact”
    1.00.6910667970.6762695160.6683274710.659957689

  12. C Use Euler’s method with step sizes \(h=0.1\), \(h=0.05\), and \(h=0.025\) to find approximate values of the solution of the initial value problem

    \[ y'+{(y+1)(y-1)(y-2)\over x+1}=0, \quad y(1)=0 \mbox{\; (Exercise~2.2.~\hspace*{-3pt}\ref{exer:2.2.14})} \]

    at \(x=1.0\), \(1.1\), \(1.2\), \(1.3\), …, \(2.0\).

    Show answer

    \(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)“Exact”
    2.0-0.772381768-0.761510960-0.756179726-0.750912371

  13. C Use Euler’s method and the Euler semilinear method with step sizes \(h=0.1\), \(h=0.05\), and \(h=0.025\) to find approximate values of the solution of the initial value problem

    \[ y'+3y=7e^{-3x},\quad y(0)=6 \]

    at \(x=0\), \(0.1\), \(0.2\), \(0.3\), …, \(1.0\). Compare these approximate values with the values of the exact solution \(y=e^{-3x}(7x+6)\), which can be obtained by the method of Section 2.1. Do you notice anything special about the results? Explain.

    Show answer
    Euler’s method
    \(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)Exact
    1.00.5388711780.5930023250.6201315250.647231889

    Euler semilinear method
    \(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)Exact
    1.00.6472318890.6472318890.6472318890.647231889

    Applying variation of parameters to the given initial value problem yields

    \(y=ue^{-3x}\), where (A) \(u'=7,\quad u(0)=6\). Since \(u''=0\), Euler’s method yields the exact solution of (A). Therefore the Euler semilinear method produces the exact solution of the given problem

    .

The linear initial value problems in Exercises 1419 can’t be solved exactly in terms of known elementary functions. In each exercise, use Euler’s method and the Euler semilinear methods with the indicated step sizes to find approximate values of the solution of the given initial value problem at 11 equally spaced points (including the endpoints) in the interval.

  1. C \(y'-2y=\dst{1\over1+x^2},\quad y(2)=2\);   \(h=0.1,0.05,0.025\) on \([2,3]\)

    Show answer
    Euler’s method
    \(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)“Exact”
    3.012.80422613513.91294466214.55962305515.282004826

    Euler semilinear method
    \(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)“Exact”
    3.015.35412228715.31725770515.29942942115.282004826

  2. C \(y'+2xy=x^2,\quad y(0)=3\)  (Exercise 2.1. 38);   \(h=0.2,0.1,0.05\) on \([0,2]\)

    Show answer
    Euler’s method
    \(x\)\(h=0.2\)\(h=0.1\)\(h=0.05\)“Exact”
    2.00.8675650040.8857192630.8950247720.904276722

    Euler semilinear method
    \(x\)\(h=0.2\)\(h=0.1\)\(h=0.05\)“Exact”
    2.00.5696707890.7208618580.8084382610.904276722

  3. C \(\dst{y'+{1\over x}y={\sin x\over x^2},\quad y(1)=2}\);  (Exercise 2.1. 39);   \(h=0.2,0.1,0.05\) on \([1,3]\)

    Show answer
    Euler’s method
    \(x\)\(h=0.2\)\(h=0.1\)\(h=0.05\)“Exact”
    3.00.9220943790.9456048000.9567528680.967523153

    Euler semilinear method
    \(x\)\(h=0.2\)\(h=0.1\)\(h=0.05\)“Exact”
    3.00.9939547540.9807513070.9741403200.967523153

  4. C \(\dst{y'+y={e^{-x}\tan x\over x},\quad y(1)=0}\);  (Exercise 2.1. 40);   \(h=0.05,0.025,0.0125\) on \([1,1.5]\)

    Show answer
    Euler’s method
    \(x\)\(h=0.0500\)\(h=0.0250\)\(h=0.0125\)“Exact”
    1.500.3198921310.3307971090.3370201230.343780513

    Euler semilinear method
    \(x\)\(h=0.0500\)\(h=0.0250\)\(h=0.0125\)“Exact”
    1.500.3055969530.3233402680.3332045190.343780513

  5. C \(\dst{y'+{2x\over 1+x^2}y={e^x\over (1+x^2)^2}, \quad y(0)=1}\);   (Exercise 2.1.  41);   \(h=0.2,0.1,0.05\) on \([0,2]\)

    Show answer
    Euler’s method
    \(x\)\(h=0.2\)\(h=0.1\)\(h=0.05\)“Exact”
    2.00.7545725600.7438698780.7383039140.732638628

    Euler semilinear method
    \(x\)\(h=0.2\)\(h=0.1\)\(h=0.05\)“Exact”
    2.00.7226104540.7277429660.7302202110.732638628

  6. C \(xy'+(x+1)y=e^{x^2},\quad y(1)=2\);  (Exercise 2.1. 42);   \(h=0.05,0.025,0.0125\) on \([1,1.5]\)

    Show answer
    Euler’s method
    \(x\)\(h=0.0500\)\(h=0.0250\)\(h=0.0125\)“Exact”
    1.502.1759599702.2102595542.2272075002.244023982

    Euler semilinear method
    \(x\)\(h=0.0500\)\(h=0.0250\)\(h=0.0125\)“Exact”
    1.502.1179533422.1798445852.2116479042.244023982

In Exercises 2022, use Euler’s method and the Euler semilinear method with the indicated step sizes to find approximate values of the solution of the given initial value problem at 11 equally spaced points (including the endpoints) in the interval.

  1. C \(y'+3y=xy^2(y+1),\quad y(0)=1\);   \(h=0.1,0.05,0.025\) on \([0,1]\)

    Show answer
    Euler’s method
    \(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)“Exact”
    1.00.0321051170.0439970450.0501593100.056415515

    Euler semilinear method
    \(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)“Exact”
    1.00.0560201540.0562439800.0563364910.056415515

  2. C \(\dst{y'-4y={x\over y^2(y+1)},\quad y(0)=1}\);   \(h=0.1,0.05,0.025\) on \([0,1]\)

    Show answer
    Euler’s method
    \(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)“Exact”
    1.028.98781665638.42695751645.36726968854.729594761

    Euler semilinear method
    \(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)“Exact”
    1.054.70913494654.72415048554.72822801554.729594761

  3. C \(\dst{y'+2y={x^2\over1+y^2},\quad y(2)=1}\);   \(h=0.1,0.05,0.025\) on \([2,3]\)

    Show answer
    Euler’s method
    \(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)“Exact”
    3.01.3614279071.3613208241.3613325891.361383810

    Euler semilinear method
    \(x\)\(h=0.1\)\(h=0.05\)\(h=0.025\)“Exact”
    3.01.2913455181.3265357371.3440041021.361383810

  4. Numerical Quadrature. The fundamental theorem of calculus says that if \(f\) is continuous on a closed interval \([a,b]\) then it has an antiderivative \(F\) such that \(F'(x)=f(x)\) on \([a,b]\) and

    \[ \int_a^bf(x)\,dx=F(b)-F(a). \tag*{\rm(A)} \]

    This solves the problem of evaluating a definite integral if the integrand \(f\) has an antiderivative that can be found and evaluated easily. However, if \(f\) doesn’t have this property, (A) doesn’t provide a useful way to evaluate the definite integral. In this case we must resort to approximate methods. There’s a class of such methods called numerical quadrature, where the approximation takes the form

    \[ \int_a^bf(x)\,dx\approx \sum_{i=0}^n c_if(x_i), \tag*{\rm(B)} \]

    where \(a=x_0<x_1<\cdots<x_n=b\) are suitably chosen points and \(c_0\), \(c_1\), …, \(c_n\) are suitably chosen constants. We call (B) a quadrature formula.

    1. Derive the quadrature formula

      \[ \int_a^bf(x)\,dx\approx h\sum_{i=0}^{n-1}f(a+ih)\mbox{\quad(where $h=(b-a)/n)$} \tag*{\rm(C)} \]

      by applying Euler’s method to the initial value problem

      \[ y'=f(x),\quad y(a)=0. \]
    2. The quadrature formula (C) is sometimes called the left rectangle rule. Draw a figure that justifies this terminology.

    3. L For several choices of \(a\), \(b\), and \(A\), apply (C) to \(f(x)=A\) with \(n = 10,20,40,80,160,320\). Compare your results with the exact answers and explain what you find.

    4. L For several choices of \(a\), \(b\), \(A\), and \(B\), apply (C) to \(f(x)=A+Bx\) with \(n=10\), \(20\), \(40\), \(80\), \(160\), \(320\). Compare your results with the exact answers and explain what you find.