Euler Method

The Euler method is the simplest, as well as historically the earliest, numerical method for differential equations. Consider a first-order initial value problem

\begin{displaymath} \frac{dy}{dx} = f(x,y), y(x_0) = y_0 \end{displaymath}

Let $y = \phi(x)$ be the solution of this problem. But we don't know $\phi(x)$, so we want to approximate it. More specifically, we want to compute a sequence $y_1$, $y_2$, ..., such that each $y_i$ is an approximation to $\phi(x_i)$, where $x_i = x_0 + i h$ for some positive constant $h$ (the step size). Let's start with the computation of $y_1$. The Euler method is based on the tangent line approximation from calculus:

\begin{displaymath} \phi(x_1) = \phi(x_0 + h) \approx \phi(x_0) + h \phi'(x_0) \end{displaymath}

Since (according to the differential equation) $\phi'(x_0) = f(x_0, y_0)$, we therefore take $y_1 = y_0 + h f(x_0, y_0)$. In terms of the direction field, starting at the initial point $[x_0, y_0]$ we go in a straight line in the direction of the ``road sign'' at that point until we reach $[x_1, y_1]$. Having $y_1$, we next want $y_2$. The same principle we used for $y_1$ would say $\phi(x_2) \approx \phi(x_1) + h \phi'(x_1)$. Now we don't know $\phi(x_1)$, we only have its approximation $y_1$, so this is what we use: $y_2 = y_1 + h f(x_1, y_1)$. And we repeat the process for each $i$:

\begin{displaymath}y_{i+1} = y_i + h f(x_i, y_i) \end{displaymath}

Let's look at an example. We'll consider the initial value problem $ y' = 2 (y^2 + 1)/(x^2+4)$, $y(0) = 1$. This equation is separable, and the solution turns out to be $y = (2+x)/(2-x)$, so we'll be able to compare the Euler approximation to the true solution. With step size $h = .1$, we have the following results: $y_1 = 1 + .1 f(0,1) = 1.1$, $y_2 = 1.1 + .1 f(.1, 1.1) = 1.1 + .2 (1.1^2+1)/ (.1^2+4) = 1.21022444$, etc. By comparison, the true solution has $\phi(.1) = 2.1/1.9 = 1.10526316$, $\phi(.2) = 2.2/1.8 = 1.22222222$, etc. Here is a table of the results of the first 10 steps, which takes us up to $x=1$. (I'm showing 9 digits of the results, but they are actually being computed in double precision, so roundoff error won't be a concern)


$ x_i$ $y_i$ $\phi(x_i)$ Error
0.0 1.00000000 1.00000000 0.0
0.1 1.10000000 1.10526316 0.00526316
0.2 1.21022444 1.22222222 0.01199778
0.3 1.33223648 1.35294118 0.02070470
0.4 1.46792616 1.50000000 0.03207384
0.5 1.61959959 1.66666667 0.04706708
0.6 1.79009854 1.85714286 0.06704432
0.7 1.98296335 2.07692308 0.09395973
0.8 2.20265794 2.33333333 0.13067539
0.9 2.45488648 2.63636364 0.18147716
1.0 2.74704729 3.00000000 0.25295271


In an effort to improve our accuracy, we might try a smaller step size. With $h = .05$, we get $y_1 = 1 + .05 f(0,1) = 1.05$. Note that now $x_1 = .05$, so we should compare this to $\phi(.05) = 2.05/1.95 = 1.05128205$. Then $y_2 = 1.05 + .05 f(.05, 1.05) = 1.10252967$ is the approximate value at $x_2 = .1$. The next table shows the values computed with $h = .05$ for the same $x$ values as before (thus leaving out the odd $ x_i$'s).


$ x_i$ $y_i$ $\phi(x_i)$ $Error$
0.0 1.00000000 1.00000000 0.0
0.1 1.10252967 1.10526316 0.00273349
0.2 1.21596496 1.22222222 0.00625726
0.3 1.34209198 1.35294118 0.01084920
0.4 1.48310373 1.50000000 0.01689627
0.5 1.64172213 1.66666667 0.02494454
0.6 1.82136643 1.85714286 0.03577643
0.7 2.02638978 2.07692308 0.05053330
0.8 2.26241822 2.33333333 0.07091511
0.9 2.53684738 2.63636364 0.09951625
1.0 2.85958887 3.00000000 0.14041113


If we compare the errors with $h = .05$ to the errors at the same points with $h = .1$, we see that we have approximately halved each of these errors (the actual ratio ranges from about $.52$ to $.555$). In the graph below we show the error in the Euler approximation at $x=1$ as a function of the step size $h$. We use 14 different values for $h$: $1/10$, $1/20$, $1/30$, $1/40$, $1/50$, $1/60$, $1/70$, $1/80$, $1/90$, $1/100$, $1/120$, $1/160$, $1/250$, $1/500$ (note that we need $1/h$ to be an integer $n$ so that $x_n$ will be 1). The data points are shown by circles. Especially for small $h$, they are very close to a straight line through the origin, representing a constant times $h$.

We could try other examples, and the net result would be similar: the error in the Euler approximation for a given $\phi(x)$ is approximately proportional to $h$. We'll see why later. But the consequence of this is that the Euler method is not suitable if you want reasonable accuracy. Suppose, in our example, we wanted the error at $x=1$ to be less than $10^{-10}$. Approximately what stepsize, and how many steps, would be needed? We had an error of $0.14041113$ with step size $.05$, and the error is approximately proportional to the step size. To divide the error by $1,404,111,300$ we would have to divide the step size by the same amount, so the step size would have to be about $.05/1404111300 \approx 3.561 \times 10^{-11}$. The number of steps would be about $1404111300/.05 \approx 2.808 \times 10^{10}$ (more than 28 billion). This would take a significant amount of time, even on a fast computer. Moreover, there would be severe problems with roundoff error, because each step would only add a number on the order of $10^{-10}$ to $y_i$. This means about 10 digits of accuracy is lost. Even in double precision, not enough digits will remain to attain the desired accuracy.


Robert
2002-01-27