f[x_, y_] := y
exactSolution[x_] := Exp[x]
x0 = 0; y0 = 1; xEnd = 2;
eulerStep[fx_, x_, y_, h_] := y + h*fx[x, y]
eulerMethod[fx_, x0_, y0_, xEnd_, h_] :=
Module[{xc = x0, yc = y0, pts = {{x0, y0}}},
While[xc < xEnd - h/2,
yc = eulerStep[fx, xc, yc, h];
xc = xc + h;
AppendTo[pts, {xc, yc}]];
pts]
pts025 = eulerMethod[f, x0, y0, xEnd, 0.25];
exactCurve = Table[{t, Exp[t]}, {t, 0, 2, 0.01}];
geomFig = ListLinePlot[{pts025, exactCurve},
PlotStyle -> {{Red, Thick}, {Blue, Thick}},
PlotLegends -> {"Euler (h = 0.25)", "Exact: e^x"},
Frame -> True, FrameLabel -> {"x", "y"},
PlotLabel -> Style["Euler method vs. exact solution", 13]];
Export["04_RK1_Eulers_Method/euler-vs-exact.pdf", geomFig];
Export["04_RK1_Eulers_Method/euler-vs-exact.svg", geomFig];4 Euler’s Method — Runge-Kutta Order 1
Every Runge–Kutta method, from the simplest to the most elaborate, is a recipe of the same form: take a weighted sum of slope estimates, multiply by the step size, and add the result to the current value. Euler’s method is the degenerate case of that recipe — one slope estimate, weight one. Deriving it carefully from the Taylor series accomplishes three things at once: we see exactly where the formula comes from, we measure how much error one step introduces, and we establish the pattern that every subsequent chapter will refine.
4.1 The Unifying Structure: A Weighted Sum of Slopes
Before deriving anything, we name the pattern that connects Euler’s method to every Runge–Kutta method in this book. Every RK method advances the solution from \(x_n\) to \(x_{n+1} = x_n + h\) by exactly one formula:
\[y_{n+1} = y_n + h\, \sum_{i=1}^{s} b_i\, k_i\]
Here \(s\) is the number of stages — the number of times the right-hand side \(f\) is evaluated per step. The \(b_i\) are weights (they must sum to 1). Each \(k_i\) is a slope estimate: a value of \(f\) computed at some point near \((x_n, y_n)\). Different choices of \(s\), the weights \(b_i\), and where the \(k_i\) are sampled produce different RK methods with different orders of accuracy. The conditions the coefficients must satisfy are what each chapter derives from the Taylor series.
For Euler’s method the answer is immediate: \(s = 1\), \(b_1 = 1\), and \(k_1 = f(x_n, y_n)\) — the slope at the current point. Substituting into the unifying formula gives the Euler step:
\[y_{n+1} = y_n + h\, f(x_n, y_n)\]
This formula will reappear in Section 4.3 as the direct output of a Taylor expansion — derived, not assumed. Seeing it here first means the algebra in Section 4.3 produces a result you already recognize.
4.2 The Geometric Idea: Following the Slope
The solution curve \(y(x)\) passes through the initial point \((x_0, y_0)\). We do not know the curve, but the ODE tells us its slope at that point: \(y'(x_0) = f(x_0, y_0)\). The tangent line at \((x_0, y_0)\) has that slope. Euler’s method says: walk along the tangent line for one step of width \(h\), then stop. The endpoint of that walk is \(y_1\).
At \((x_1, y_1)\) the ODE supplies a new slope — \(f(x_1, y_1)\) — and we repeat. Each step follows the tangent line at the current approximate point. Because a tangent line is only a first-order approximation to the curve, the solution drifts away from the truth, but for small \(h\) the drift is slow. The plot below shows what this looks like on our benchmark problem \(y' = y\), \(y(0) = 1\), whose exact solution is \(e^x\).
The red polyline is the Euler approximation; the blue curve is the exact solution. Each segment of the polyline is one tangent-line step. Notice that the approximation consistently falls below the exact curve: the exact solution is concave up (\(y'' = e^x > 0\)), so the tangent line always undershoots. Section 4.6 will show, side by side, how the error shrinks as \(h\) decreases.
4.3 Deriving Euler’s Method from the Taylor Series in Wolfram Language
Chapter 3 established that every RK derivation is a Taylor expansion in the step size \(h\). We apply that idea directly: expand \(y(x + h)\) as a power series around \(h = 0\), use the ODE to eliminate \(y'(x)\), and truncate. The result is the Euler formula. The full Taylor series is:
\[y(x + h) = y(x) + h\, y'(x) + \frac{h^2}{2}\, y''(x) + O(h^3)\]
We ask Series[] to produce this expansion. We treat \(y\) as an unspecified smooth function — Wolfram Language will keep \(y'[x]\) as a formal derivative because it does not yet know the ODE:
Series[y[x + h], {h, 0, 1}]
The output is y[x] + h y'[x] + O[h]^2. Now we supply the ODE. The rule \(y' = f(x, y)\) tells us \(y'[x]\) should be replaced by \(f[x, y[x]]\). We apply that substitution and drop the remainder with Normal[]:
Normal[Series[y[x + h], {h, 0, 1}]] /. y'[x] -> f[x, y[x]]
The output is y[x] + h f[x, y[x]]. Translating the continuous notation to the discrete stepping notation — \(x\) becomes \(x_n\), \(y[x]\) becomes \(y_n\), and \(y[x + h]\) becomes \(y_{n+1}\) — gives exactly the Euler formula from Section 4.1:
\[y_{n+1} = y_n + h\, f(x_n, y_n)\]
This formula was not assumed — it fell out of the Taylor series when we kept only the \(O(h)\) term and used the ODE to replace \(y'(x_n)\). The term we discarded, starting at \(\tfrac{h^2}{2}\, y''(x_n)\), is the local truncation error. Section 4.4 measures it precisely.
4.4 Local Truncation Error: How Much We Lose Each Step
The local truncation error (LTE) is the gap between the exact value \(y(x_n + h)\) and the Euler approximation \(y_{n+1}\), measured in one isolated step that starts from an exact value of \(y_n\). It quantifies how much a single step departs from the truth before errors begin to interact.
We already have both quantities. The exact value is the full Taylor series; the Euler value keeps only the first two terms. The LTE is everything we dropped. We ask Wolfram Language to expand to one order higher to make the dropped term visible:
Series[y[x + h], {h, 0, 2}]
The output includes the term \((1/2)\, y''[x]\, h^2\). That is the leading term of the LTE. Since \(y''(x_n)\) does not depend on \(h\) (it is determined by the ODE and the exact solution at \(x_n\)), the \(h^2\) factor tells us the LTE is \(O(h^2)\):
\[\text{LTE} = \frac{h^2}{2}\, y''(x_n) + O(h^3)\]
For the benchmark \(y' = y\) we have \(y''(x) = e^x\), which is positive on the whole interval \([0, 2]\). That is why the Euler approximation always undershoots: the LTE is always positive. Halving \(h\) reduces the LTE by a factor of four. This is a local statement about one isolated step. The next section asks what happens when we chain many steps together.
4.5 Global Error: How Mistakes Accumulate
The global error is the total deviation from the exact solution after integrating from \(x_0\) to \(x_\text{end}\). We take \(N = (x_\text{end} - x_0) / h\) steps, so \(N\) grows as \(h\) shrinks. A rough but instructive argument: if each step introduces an LTE of order \(h^2\) and the errors accumulate additively, the total is proportional to \(N\) times the LTE per step.
\[\text{Global Error} \approx N \cdot \text{LTE} = \frac{x_\text{end} - x_0}{h} \cdot O(h^2) = O(h)\]
One power of \(h\) in the numerator of LTE cancels the \(1/h\) from the step count, and the interval length \((x_\text{end} - x_0)\) is a fixed constant. The result is a global error proportional to \(h\). Halving \(h\) halves the global error — but it also doubles the number of steps, so the computation doubles. Euler’s method is therefore a first-order method.
Critical distinction: LTE is \(O(h^2)\) per step, but global error is \(O(h)\). A method’s order is defined by its global error, not its LTE. A reader who sees the \(h^2\) in the LTE and concludes the method is second order is making a common mistake. The order drops by one because the error accumulation over \(N = O(1/h)\) steps costs us one power of \(h\).
4.6 A Worked Example: Numerical vs. Exact Solution
We apply Euler’s method to the benchmark IVP from Chapter 2: \(y' = y\), \(y(0) = 1\), exact solution \(y(x) = e^x\), integrated over \([0, 2]\). We reload the definitions so this notebook runs independently.
f[x_, y_] := y
exactSolution[x_] := Exp[x]
x0 = 0; y0 = 1; xEnd = 2;eulerStep translates the Euler formula directly. It takes the slope function fx, the current position x, the current value y, and the step size h, and returns the next value:
eulerStep[fx_, x_, y_, h_] := y + h * fx[x, y]eulerMethod chains eulerStep from \(x_0\) to \(x_\text{end}\) and collects all solution points. The internal variables xc and yc (current \(x\), current \(y\)) are localized in Module so they do not pollute the global namespace. The While loop stops when xc reaches \(x_\text{end}\) to within half a step, which avoids floating-point overshoot:
eulerMethod[fx_, x0_, y0_, xEnd_, h_] :=
Module[{xc = x0, yc = y0, pts = {{x0, y0}}},
While[xc < xEnd - h/2,
yc = eulerStep[fx, xc, yc, h];
xc = xc + h;
AppendTo[pts, {xc, yc}]];
pts]We run the method with step size \(h = 0.25\) and compare the Euler value at \(x = 2\) to the exact answer \(e^2 \approx 7.389\):
pts = eulerMethod[f, x0, y0, xEnd, 0.25];
yEuler = Last[pts][[2]];
TableForm[
{{"Euler y(2)", "Exact y(2)", "Absolute error"},
{yEuler, N[Exp[2], 6], N[Abs[yEuler - Exp[2]], 4]}},
TableHeadings -> None]
4.7 Visualizing the Error as Step Size Changes
The figure below shows Euler approximations at four step sizes overlaid on the exact solution. Two things to watch: how the red polyline approaches the blue exact curve as \(h\) decreases, and how the absolute error at \(x = 2\) shrinks. If global error is truly \(O(h)\), halving \(h\) should halve the error.
hPanels = {0.5, 0.25, 0.125, 0.0625};
errFigs = Table[
Module[{pts, yApprox, err, exactCurve},
pts = eulerMethod[f, 0, 1, 2, h];
yApprox = Last[pts][[2]];
err = Abs[yApprox - Exp[2]];
exactCurve = Table[{t, Exp[t]}, {t, 0, 2, 0.005}];
ListLinePlot[{pts, exactCurve},
PlotStyle -> {{Red, Thick}, {Blue, Thick}},
Frame -> True, FrameLabel -> {"x", "y"},
PlotRange -> {{0, 2.1}, {0.5, 8.5}},
PlotLabel -> Style[
"h = " <> ToString[h] <> ", err at x=2: " <>
ToString[NumberForm[err, {5, 4}]], 11],
ImageSize -> 280]],
{h, hPanels}];
errGrid = GraphicsGrid[Partition[errFigs, 2], Spacings -> Scaled[0.04]];
Export["04_RK1_Eulers_Method/euler-h-grid.pdf", errGrid];
Export["04_RK1_Eulers_Method/euler-h-grid.svg", errGrid];
To see the \(O(h)\) convergence numerically, we compute the error at \(x = 2\) for a sequence of step sizes, each half the previous. A first-order method produces error ratios near \(0.5\) each time \(h\) is halved:
hValues = {0.5, 0.25, 0.125, 0.0625, 0.03125};
errors = Table[
Abs[Last[eulerMethod[f, 0, 1, 2, hv]][[2]] - Exp[2]],
{hv, hValues}];
ratios = Prepend[Rest[errors]/Most[errors], "---"];
TableForm[
Transpose[{hValues, N[errors, 4], ratios}],
TableHeadings -> {None, {"h", "Error at x=2", "Ratio"}}]
The Ratio column should be close to \(0.5\) at every row, confirming that the global error is proportional to \(h\). The values will not be exactly \(0.5\) because the constant hidden in the \(O(h)\) notation depends on the specific problem, but the convergence is clear. Chapter 8 will plot the same data on a log-log scale, where first-order convergence appears as a straight line of slope 1.
4.8 The Limitations That Motivate Everything That Follows
Euler’s method is first-order: halving \(h\) halves the global error and doubles the work. To reduce the error by a factor of 100 we need 100 times as many steps. On smooth problems over short intervals this is merely slow; on stiff problems or long-time integrations it can be completely impractical.
The deeper issue is structural: Euler uses only the slope at the left endpoint of each interval. That slope is the best first-order estimate of the average slope over \([x_n, x_{n+1}]\), but it is only that. The exact average slope over the interval is what we would need to advance the solution without error. Euler approximates the average with a single sample; every higher-order method improves that approximation by sampling at more points inside the interval.
This is the idea behind Runge–Kutta order 2. Instead of one slope estimate at the left endpoint, RK2 uses two estimates — one at the left endpoint and one near the midpoint or right endpoint — and takes a weighted average. The weights must satisfy conditions derived from the same Taylor expansion we just used. Chapter 5 derives those conditions, finds that they admit a whole family of solutions, and implements two of them: the Midpoint Method and Heun’s Method.