Calculus 2

Euler's method

Most differential equations have no solution formula. How can the solution’s value at a point still be found, how far off is the answer, and in which direction?

You need

Slope fieldsCalculus 1 Concavity and inflection pointsCalculus 1

Where you are 7 in Differential equations

Calculus8

Differential equations7

▴ 1 above Differential equations and their solutions Slope fields

Euler's method

Separable differential equations Exponential growth and decay models ▾ Show all 7 (1 more)

Why this exists

A slope field shows the shape of the solutions of \dfrac{dy}{dx} = f(x, y), and separation of variables gives a formula for some of them. Most equations allow neither a formula nor an exact reading: the solutions of \dfrac{dy}{dx} = x - y^2, for example, cannot be written with the functions of this course. A model still has to answer “what is y at x = 2?” with a number.

Euler’s method, described by Leonhard Euler in the eighteenth century, turns “follow the segments” into arithmetic. Because it only approximates, it comes with two questions, and they are the point of this lesson: how the error shrinks as the steps shrink, and whether the estimate lies above or below the truth.

Definition

Start at the known point and walk along the tangent line for a short run h. The equation gives the slope of that tangent, so the rise is slope times run. At the new point, compute the slope again and take another step. Each step is the tangent-line estimate of Tangent lines and normal lines, restarted from wherever the last step landed.

In symbols, for the initial value problem \dfrac{dy}{dx} = f(x, y), y(x_0) = y_0, with step size h:

x_{n+1} = x_n + h, \qquad y_{n+1} = y_n + f(x_n, y_n)\,h.

The number y_n is the Euler estimate of y(x_n), and y(x_n) - y_n is its error. A negative h steps to the left. The work is usually laid out in a table with columns x_n, y_n, the slope f(x_n, y_n) and the change f(x_n, y_n)\,h.

Intuition

Figure 1 shows two runs for the same equation. Each straight piece starts tangent to a solution curve, but the curve bends away from its tangent, so at the end of the step the point is off the curve. The next slope is then computed at a point that is slightly wrong, so the method follows a neighboring solution curve for a step, and so on. Two things follow. Shorter steps leave less room for the bending, so the error falls as h falls. And the direction of the bending, the concavity, decides which side each step lands on.

Figure 1Euler’s method for \dfrac{dy}{dx} = 0.5y, y(0) = 3 (Example 2). Each straight piece follows a tangent line; the curve bends up away from it, so both broken lines run low, and the one with half the step size ends closer, 0.83 below the curve against 1.40.

Results

Euler’s method is first order. Suppose f, \dfrac{\partial f}{\partial x} and \dfrac{\partial f}{\partial y} are continuous near the solution on [x_0, X], and N steps of size h = \dfrac{X - x_0}{N} reach X. There is a constant C, depending on f and the interval but not on h, with, for every small enough h,

\lvert y(X) - y_N \rvert \le C\,h.

For small h the error is close to a constant multiple of h: halving the step roughly halves the error, at the cost of twice as many steps. This is stated without full proof; the argument below shows where the h comes from.

Concavity decides the lean. Compute \dfrac{d^2y}{dx^2} from the equation by differentiating f(x, y) with y a function of x and replacing \dfrac{dy}{dx} by f, as in Slope fields. If \dfrac{d^2y}{dx^2} > 0 throughout a region that holds the solution curve, the broken line of Euler steps and everything between them, every estimate is low: y_n < y(x_n) for n \ge 1. If \dfrac{d^2y}{dx^2} < 0 there, every estimate is high. If the concavity changes along the way, the sign of the error can change too.

Why it is true

The lean. Let u be the solution curve through (x_n, y_n). One step runs along the tangent to u at that point, and a concave up curve lies above its tangent lines (Concavity and inflection points), so y_{n+1} < u(x_{n+1}). Suppose y_n \le y(x_n), which holds at n = 0. Solution curves never cross when f and \dfrac{\partial f}{\partial y} are continuous (stated without proof in Slope fields), so u, starting at or below the true solution, stays at or below it, and y_{n+1} < u(x_{n+1}) \le y(x_{n+1}). Induction gives every n. The concave down case reverses each inequality.

First order. On one step the curve u and its tangent differ by

\begin{aligned} &u(x_n + h) - \big(u(x_n) + u'(x_n)h\big) \\ &\quad = \int_{x_n}^{x_n + h} \big(u'(s) - u'(x_n)\big)\,ds \\ &\quad \approx \int_{x_n}^{x_n + h} u''(x_n)(s - x_n)\,ds \\ &\quad = \tfrac{1}{2}u''(x_n)\,h^2. \end{aligned}

If \lvert y'' \rvert \le M near the path, each step adds at most about \tfrac{1}{2}Mh^2, and there are \dfrac{L}{h} steps, L = X - x_0, so the total is about \tfrac{1}{2}MLh: proportional to h. What this leaves out is that an error made early is carried forward and grows or shrinks on the way, by a factor that depends on \dfrac{\partial f}{\partial y}. The full proof bounds that factor and is done in courses on numerical analysis.

Worked examples

Example 1

Liko’s soup, step by step

Liko turns off the burner under a pot of soup at 94^\circ\text{C} and leaves it on the stove top in a kitchen at 22^\circ\text{C}. Newton’s law of cooling (Exponential growth and decay models) with his constant gives \dfrac{dT}{dt} = -0.04(T - 22), with t in minutes. Use Euler’s method with h = 5 to estimate T(15).

Solution

At the start the soup is cooling at -0.04(94 - 22) = -2.88 degrees per minute. Following that rate for 5 minutes changes the temperature by 5(-2.88) = -14.4 degrees, so T(5) \approx 79.6. The rate is now recomputed at the new point, -0.04(79.6 - 22) = -2.304: cooler soup cools more slowly, and the method notices only at the start of each step.

t_n	T_n	slope	change
0	94	-2.88	-14.4
5	79.6	-2.304	-11.52
10	68.08	-1.8432	-9.216
15	58.864

So T(15) \approx 58.9^\circ\text{C}. Units check: degrees per minute times minutes is degrees. The exact value 22 + 72e^{-0.6} \approx 61.5 is higher; Example 4 explains why.

Example 2

Ilyas halves the step

Ilyas estimates y(2) for \dfrac{dy}{dx} = 0.5y, y(0) = 3, with h = 1, 0.5, 0.25 and 0.125, and compares with the exact solution y = 3e^{0.5x}.

Solution

With h = 1: y_1 = 3 + 1(1.5) = 4.5 and y_2 = 4.5 + 1(2.25) = 6.75. For this equation every step is y_{n+1} = y_n + 0.5y_n h = (1 + 0.5h)y_n, so N steps multiply 3 by (1 + 0.5h)^N. With h = 0.5 that is 3(1.25)^4 = 7.3242; the other rows are found the same way. The exact value is 3e \approx 8.1548.

h	estimate	error	\text{error} \div h
1	6.75	1.4048	1.40
0.5	7.3242	0.8306	1.66
0.25	7.6974	0.4575	1.83
0.125	7.9138	0.2411	1.93

Each halving divides the error by 1.69, then 1.82, then 1.90, closer to 2 each time, and the last column levels off: for small h the error is about 0.75e \cdot h \approx 2.04h. All four estimates are low, as Figure 1 shows. That is the concavity rule: \dfrac{d^2y}{dx^2} = 0.5\dfrac{dy}{dx} = 0.25y > 0 while y > 0.

Example 3

Amina predicts the lean on a circle

Amina takes \dfrac{dy}{dx} = -\dfrac{x}{y}, y(0) = 2. Before computing, decide whether Euler estimates will be high or low; then estimate y(1.5) with h = 0.5 and compare with the exact solution \displaystyle y = \sqrt{4 - x^2}, the top of the circle x^2 + y^2 = 4.

Solution

By the quotient rule, then replacing y' by -\dfrac{x}{y},

\frac{d^2y}{dx^2} = -\frac{y - xy'}{y^2} = -\frac{y + \frac{x^2}{y}}{y^2} = -\frac{x^2 + y^2}{y^3}.

This is negative wherever y > 0, so every solution curve in the upper half plane is concave down and the estimates will be high.

x_n	y_n	slope	change
0	2	0	0
0.5	2	-0.25	-0.125
1	1.875	-0.5333	-0.2667
1.5	1.6083

The exact value is \displaystyle \sqrt{1.75} \approx 1.3229, so the estimate is high by about 0.285, as predicted (Figure 2). The error is large because \lvert y'' \rvert grows quickly as y falls toward 0: the circle turns faster than steps of 0.5 can follow.

Figure 2Euler’s method for \dfrac{dy}{dx} = -\dfrac{x}{y}, y(0) = 2 with h = 0.5 (Example 3). The solution curve is concave down and every estimate is high. The dashed step past x = 2 goes where no solution exists (Boundaries).

Example 4

Liko combines two runs

Estimate the soup’s temperature after 10 minutes in Example 1 with h = 2.5, decide from concavity whether Euler estimates run high or low, and combine the runs with h = 5 and h = 2.5 into a better estimate.

Solution

Each step is T_{n+1} - 22 = (1 - 0.04h)(T_n - 22), so with h = 2.5 each step multiplies the gap above room temperature by 0.9. Four steps give 72(0.9)^4 = 47.2392, so T(10) \approx 69.2392, against 68.08 with h = 5.

Differentiating the equation, \dfrac{d^2T}{dt^2} = -0.04\dfrac{dT}{dt} = 0.0016(T - 22), positive while the soup is above room temperature. The solution is concave up, so both estimates are low and the true value is above 69.24.

If the error is about Ch, the true value is about 68.08 + 5C and also about 69.2392 + 2.5C. Subtracting, 2.5C \approx 69.2392 - 68.08 = 1.1592, so

T(10) \approx 69.2392 + 1.1592 = 70.3984.

The exact value 22 + 72e^{-0.4} \approx 70.2630 is off from this by 0.135, against 1.024 for the better Euler run. Combining two runs as 2E(\tfrac{h}{2}) - E(h) cancels most of the first-order error; it is called Richardson extrapolation.

Example 5

Ilyas’s equation with no formula

Ilyas needs y(2) for \dfrac{dy}{dx} = x - y^2, y(0) = 0, which is not separable and has no solution formula in the functions of this course. Estimate it with h = 0.5, and decide whether the lean can be predicted.

Solution

x_n	y_n	slope	change
0	0	0	0
0.5	0	0.5	0.25
1	0.25	0.9375	0.46875
1.5	0.71875	0.98340	0.49170
2	1.21045

So y(2) \approx 1.2104. Concavity: \dfrac{d^2y}{dx^2} = 1 - 2y\dfrac{dy}{dx} = 1 - 2y(x - y^2), which is 1 at the start and about -0.37 near (2, 1.19). The solution changes from concave up to concave down near x \approx 1.2 (Figure 3), so no single lean is guaranteed.

The numbers bear this out. With h = 0.25 the estimates are y(1) \approx 0.3653 and y(2) \approx 1.1978; a far more accurate method gives y(1) \approx 0.4555 and y(2) \approx 1.1936. At x = 1, still in the concave up part, every estimate is low. At x = 2 every estimate is high: the high-leaning steps after the inflection point outweigh the earlier low-leaning ones.

Figure 3Euler’s method for \dfrac{dy}{dx} = x - y^2, y(0) = 0 with h = 0.5 (Example 5). The estimates run low while the solution is concave up and end high after the inflection point, the dashed line.

Boundaries

The method needs a solution to follow. With f and \dfrac{\partial f}{\partial y} continuous near the path, a unique solution passes through each point. The solution of Example 3 ends at x = 2, where y = 0, the tangent is vertical and -\dfrac{x}{y} is undefined. The method does not notice: from (2, 1.1420) it steps to (2.5, 0.2664), the next slope is about -9.39, and y(3) comes out near -4.43. These numbers have no solution behind them.
First order describes small steps. For Liko’s soup each step multiplies the gap T - 22 by 1 - 0.04h. With h = 25 the gap becomes 0 in one step, and with h = 50 the factor is -1: the estimate after 50 minutes is 22 - 72 = -50^\circ\text{C}, though real soup never cools below the room. Only steps short compared with how fast the slope changes behave like Ch.
The concavity rule needs one sign of \dfrac{d^2y}{dx^2} over the whole region holding the solution, the Euler steps and the space between them, not only at the start. Where the sign changes, the error can change sign (Example 5).
A computer cannot shrink h forever. Each step rounds, and many tiny steps pile up rounding error, so beyond some point a smaller h stops helping.

Where this goes next

When f depends on x alone, Euler’s method is a left Riemann sum (The definite integral as a limit of Riemann sums, Exercise 15). Better methods use more than one slope per step: averaging the slopes at both ends, as the trapezoidal rule averages heights (The trapezoidal rule and Simpson's rule), gives error proportional to h^2, and the Runge-Kutta method, used to draw the true curve in Figure 3, gives h^4. The local error \tfrac{1}{2}y''h^2 is the first omitted term of a Taylor and Maclaurin series. Logistic growth changes concavity at half its ceiling, a natural place to watch an Euler estimate switch from low to high.

Exercises

Running the method

Use Euler’s method with h = 0.5 to estimate y(1) for \dfrac{dy}{dx} = x + y, y(0) = 1.
Use Euler’s method with h = 0.5 to estimate y(1.5) for \dfrac{dy}{dx} = y - x, y(0) = 2.
Amina’s soup is at 90^\circ\text{C} when her burner goes off, in a kitchen at 24^\circ\text{C}, and \dfrac{dT}{dt} = -0.03(T - 24) with t in minutes. Estimate T(20) with h = 10.
For \dfrac{dy}{dx} = y^2, y(1) = 1, take two steps with h = -0.5 to estimate y(0). Compare with the exact solution y = \dfrac{1}{2 - x}.
Use h = 0.5 to estimate y(2) for \dfrac{dy}{dx} = x^2, y(0) = 0. Which sum for \int_0^2 x^2\,dx have you computed, and what is the exact value?

Step size and error

For \dfrac{dy}{dx} = -y, y(0) = 1, estimate y(1) with h = 0.5 and with h = 0.25. Find each error against e^{-1} and the factor by which halving cut it.
A run with h = 0.1 on [0, 2] has error 0.06 at x = 2. Predict the error with h = 0.025, and about how many steps would bring the error below 0.001.
Euler’s method gives 3.412 with h = 0.2 and 3.506 with h = 0.1 for the same value. Estimate the true value by combining them, and say what the two runs suggest about the concavity of the solution.
For \dfrac{dy}{dx} = 2x, y(0) = 0, take n steps of h = \dfrac{1}{n} to x = 1. Show the estimate is 1 - h, and show that each step falls short by exactly h^2.

Over or under

Without solving, decide whether Euler estimates for \dfrac{dy}{dx} = 1 + y^2, y(0) = 1 are high or low. Check with one step of h = 0.5 against the exact solution y = \tan\left(x + \dfrac{\pi}{4}\right).
For \dfrac{dy}{dx} = 3 - y, y(0) = 0, estimate y(2) with h = 1, predict the lean from concavity, and compare with the exact y = 3 - 3e^{-x}.
For \dfrac{dy}{dx} = xy, y(0) = -1, estimate y(1) with h = 0.5, predict the lean, and compare with the exact y = -e^{x^2/2}.
For \dfrac{dy}{dx} = \cos x, y(0) = 0, use h = \dfrac{\pi}{2} to estimate y(\pi) and y(2\pi). Compare with \sin x and explain the two errors from concavity.

Reasoning

Explain why one Euler step from a point where the solution curve through it is concave up lands below that curve. What more is needed to conclude that every estimate in a long run is low?
Show that for \dfrac{dy}{dx} = f(x), y(a) = y_0, Euler’s method gives y_0 plus a left Riemann sum. For increasing f, show that the concavity rule and the behavior of left sums agree.
For \dfrac{dT}{dt} = k(T - A) with k < 0 and T(0) = T_0 > A, show that T_n = A + (T_0 - A)(1 + kh)^n, and that the estimates are low whenever 0 < h < -\dfrac{1}{k}.
Estimate y(1) for \dfrac{dy}{dx} = -\dfrac{x}{y}, y(0) = 1 with h = 0.5 and compare with the exact solution. Why is the error so large?
In Example 5, halving h from 0.5 to 0.25 cut the error at x = 2 by a factor of about 4, not 2. Explain why, and why this does not contradict the first-order result.

Answers

2.5. At (0, 1) the slope is 1, so y_1 = 1.5; at (0.5, 1.5) it is 2, so y_2 = 1.5 + 1 = 2.5.
5.875. Slopes 2, 2.5, 3.25 at (0, 2), (0.5, 3), (1, 4.25); each change is half the slope, giving 3, 4.25, 5.875.
T(20) \approx 56.34^\circ\text{C}. Slope -0.03(66) = -1.98, change -19.8, T_1 = 70.2; slope -0.03(46.2) = -1.386, change -13.86, T_2 = 56.34.
y(0) \approx 0.375. At (1, 1) the slope is 1 and the change -0.5, so y(0.5) \approx 0.5; there the slope is 0.25 and the change -0.125. The exact value is \dfrac{1}{2} = 0.5.
1.75. Each step adds 0.5x_n^2: 0.5(0 + 0.25 + 1 + 2.25) = 1.75, the left Riemann sum with four pieces. The exact value is \dfrac{8}{3} \approx 2.667.
Each step multiplies by 1 - h: (0.5)^2 = 0.25 and (0.75)^4 \approx 0.3164. With e^{-1} \approx 0.3679 the errors are 0.1179 and 0.0515, cut by a factor of about 2.29.
About 0.015: two halvings divide the error by about 4. Error 0.001 needs h \approx 0.1 \cdot \dfrac{0.001}{0.06} \approx 0.00167, about 1{,}200 steps.
2(3.506) - 3.412 = 3.600. The estimates rise as h shrinks, so they approach from below, which is what a concave up solution produces.
y_n = \sum_{k=0}^{n-1} 2(kh)h, which is 2h^2 \cdot \dfrac{n(n-1)}{2} = \dfrac{n - 1}{n} = 1 - h, and y(1) = 1, so the error is exactly h. On one step the true change is \int_{x_k}^{x_k + h} 2x\,dx = 2x_kh + h^2 and Euler’s is 2x_kh; n shortfalls of h^2 total nh^2 = h.
Low. y'' = 2yy' = 2y(1 + y^2) > 0 while y > 0, and y' \ge 1 keeps y positive. One step: 1 + 0.5(2) = 2, against \tan\left(0.5 + \dfrac{\pi}{4}\right) \approx 3.408.
y_1 = 3, then the slope is 0, so y(2) \approx 3. y'' = -y' = -(3 - y) < 0 for y < 3: high. Indeed 3 - 3e^{-2} \approx 2.594. Euler lands on the equilibrium y = 3 in one step; the solution only approaches it.
y(1) \approx -1.25: the slope at (0, -1) is 0, then at (0.5, -1) it is -0.5, change -0.25. y'' = y + xy' = y(1 + x^2) < 0 for y < 0: high. Indeed -e^{0.5} \approx -1.6487 < -1.25.
y(\pi) \approx \dfrac{\pi}{2}(1 + 0) \approx 1.571, against \sin \pi = 0; y(2\pi) \approx \dfrac{\pi}{2}(1 + 0 - 1 + 0) = 0 = \sin 2\pi. The sine curve is concave down on (0, \pi), so the first two steps lean high; it is concave up on (\pi, 2\pi), so the last two lean low and here cancel the excess exactly.
The step follows the tangent line at its start, and a concave up curve lies above its tangents. For a long run, the solution curves through every Euler point must be concave up too, since each step follows a different one, and solution curves must not cross, so a curve starting below the true solution stays below it.
y_{n+1} = y_n + f(x_n)h gives y_N = y_0 + \sum_{n=0}^{N-1} f(x_n)h, while y(X) = y_0 + \int_a^X f(x)\,dx. For increasing f a left sum is too small; and y'' = f'(x) > 0, concave up, so the rule also says low.
T_{n+1} - A is (T_n - A) + kh(T_n - A), that is, (1 + kh)(T_n - A), which gives the formula. The exact value at t = nh is A + (T_0 - A)(e^{kh})^n. With u = kh between -1 and 0, 0 < 1 + u < e^u (the tangent to the concave up e^u at 0 lies below it), so (1 + u)^n < e^{un} and T_n is below the true value.
y_1 = 1 (slope 0), then the slope at (0.5, 1) is -0.5, so y(1) \approx 0.75. The exact solution \displaystyle \sqrt{1 - x^2} is 0 there: high, as concavity predicts, but by 0.75, because y'' = -\dfrac{1}{y^3} grows without bound as the curve turns vertical at x = 1.
The error at x = 2 is a low-leaning part from before x \approx 1.2 and a high-leaning part after it, of opposite signs. Each part shrinks roughly like h, but at x = 2 their first-order terms nearly cancel: for small h the error is only about 0.007h. At h = 0.5 that term is about 0.004, while the error is 0.017, so most of the error comes from terms like h^2, which the first-order statement ignores and which shrink by 4 on halving. The first-order result bounds the error by Ch and describes small h; with h = 0.125 the error is 0.0014, a factor of only about 3, and the factor approaches 2 as h shrinks further.