The truncation error does not satisfy that equation, it's just its definition.
Consider two following problems:
- The first is an ODE.
$$
y'(t) = f(t, y(t))\\
y(0) = a.
$$
Its solution is some smooth function $y(t)$.
- The second is a difference equation
$$
\frac{z_{i+1} - z_i}{h} = f(t_i, z_i)\\
z_0 = a.
$$
Its solution is some discrete function $z_i$.
I've intentionally used different letters to denote those two solutions. They are quite different, the former is a smooth function while the latter is a discrete one. One needs to be careful even to compare those two. Usually the third function is introduced. It is defined as a restriction of the smooth $y(t)$ to the grid $t_i$, where the discrete function $z_i$ is defined. Let's denote the restriction as $w_i$:
$$
w_i \equiv y(t_i).
$$
The function $w_i$ is discrete just like $z_i$ and $w_i$ coincide with $y(t)$ at grid points. Since now $w_i$ and $z_i$ are functions of the same class we can easily compare them:
$$
e_i = w_i - z_i \equiv y(t_i) - z_i.
$$
So, roughly speaking, the global error shows how close are $y(t)$ and $z_i$ (by restricting the former to the grid). When someone is solving some problem numerically the global error is what he is interesting in. Anyway, direct computation of global error is almost impossible, since we often simply do not have the exact values of $w_i = y(t_i)$ (
in contradistinction to $z_i$, which we can compute easily).
And the local truncation error concept comes to the rescue. Note that previously we've compared the solutions. Now we're going to compare problems. Take $z_i$. It is the solution to the second problem. Plugging $z_i$ into it makes it a valid identity
$$
\frac{z_{i+1} - z_i}{h} = f(t_i, z_i)\\
z_0 = a.
$$
But if we now take $w_i$ and try to plug it into the difference scheme we wont get an identity. Instead we'll get a residual:
$$
\frac{w_{i+1} - w_i}{h} = f(t_i, w_i) \color{red}{{}+ d_i}\\
w_0 = a \color{red}{{} + d_0}.
$$
If we are very lucky, some residuals may vanish, like $d_0$, but often it is not the case.
So why is $d_i$ interesting while it also is defined in terms of $w_i$ (the unknown solution to the original problem)? It turns out that we can estimate the $d_i$ without knowing the exact values of $w_i$ by just knowing the original problem.
$$
d_i = \frac{w_{i+1} - w_i}{h} - f(t_i, w_i) \equiv
\frac{y(t_{i+1}) - y(t_i)}{h} - f(t_i, y(t_i)) = \\ =
y'(t_i) + h \frac{y''(t_i)}{2} + O(h^2) - f(t_i, y(t_i)) = \\ =
\color{blue}{\left[y'(t_i) - f(t_i, y(t_i))\right]} + \color{red}{h \frac{y''(t_i)}{2} + O(h^2)}
$$
The blue term in braces is exactly the original ODE, and $y(t)$ is exactly its solution. So the term is equal to zero.
$$
d_i = h \frac{y''(t_i)}{2} + O(h^2).
$$
Similar result may be obtained if using different form of Taylor's formula:
$$
d_i = h \frac{y''(\xi_i)}{2}, \qquad \xi_i \in [t_{i}, t_{i+1}].
$$
So now we can estimate the local truncation error, but we're interested in the global error.
To relate them we need to introduce another concept of stability. Consider the two discrete problems
$$
\begin{aligned}
&\frac{z_{i+1} - z_i}{h} = f(t_i, z_i)\\
&z_0 = a
\end{aligned}
\qquad\text{and}\qquad
\begin{aligned}
&\frac{w_{i+1} - w_i}{h} = f(t_i, w_i) \color{green}{{} + d_i}\\
&w_0 = a \color{green}{{} + d_0}
\end{aligned}.
$$
Pretend that we know $d_i$. Let's view the second problem as a perturbation of the first one. That's reasonable, since $d_i$ is a small value of $O(h)$ magnitude. A difference problem is called stable if such small perturbations result in small changes of the solution. For this case this means that the difference $z_i - w_i$ will also be small. Precisely
$$
\max_i |z_i - w_i| \leq C \max_i |d_i|
$$
where $C$ is called the stability constant of the method. For the explicit Euler method it can be shown that for Lipschitz-continuous $f$
$$
C \leq e^{LT}
$$
with $L$ being the Lipschitz constant of $f$ and $T$ is the total integration time $T = \max_i t_i$.
Finally we can relate the global error and the local truncation error by
$$
|e_i| \leq C \max_i |d_i|
$$
If the local truncation error tends to zero when the discrete mesh is refined the numerical method is called consistent. The Lax theorem states that a stable consistent method converges, in sense that $e_i \to 0$ when the mesh is refined.
Let's look at the half axis $y=0$, $t>0$. There the right side is $f(t,0)=t>0$ so that no solution may cross from the upper to the lower quadrant.
Close to zero one gets $y(t)=\frac12t^2+O(t^9)$ so that the solution will indeed enter the upper quadrant from the start. Furthermore, from $y'(t)\le t$ we get $y(t)\le\frac12t^2$, so that we also know an upper bound for the solution. We can restrict the region for the estimates of the Euler method to $(t,x)\in[0,1]\times[0,1]$, or, if you want to be cautious, $(t,x)\in[0,1]\times[-1,1]$.
On that region, $$|f(t,y)|\le 1=M_1$$ is a bound for the first derivative of any solution, and $$|f_t+f_yf|=|1-4y^3(t-y^4)|\le 5=M_2$$ a bound for the second derivative. Also, for the $y$-Lipschitz constant one gets similarly $$|f_y|=|-4y^3|\le 4=L. $$
Now insert into the error estimate
$$
e(t,h)\le \frac{M_2}{2L}(e^{Lt}-1)h=\frac{5}{8}(e^{4t}-1)h.
$$
Now let's see how that bound stands up to the actual error of the numerical method
def EulerIntegrate(f,y0,t):
y = np.asarray(len(t)*[y0], dtype=float);
for k in range(len(t)-1):
y[k+1] = y[k]+(t[k+1]-t[k])*f(t[k],y[k])
return y
def f(t,y): return t-y**4
t = np.linspace(0,1,64*15+1)
N = [30*2**k for k in range(4)]
yEuler = [ EulerIntegrate(f,0,t[::n]) for n in N];
yGold = odeint(f,0,t,tfirst=True, atol=1e-12, rtol=1e-13)[:,0]
fig,ax = plt.subplots(1,2,figsize=(12,5))
for k in range(len(N)):
h = t[N[k]]-t[0];
ax[0].plot(t[::N[k]], yEuler[k],'-o',label="h=%.4f"%(1./N[k]))
ax[1].plot(t[::N[k]], (yEuler[k]-yGold[::N[k]])/h,'-o',label="h=%.4f"%(1./N[k]))
gives the plots
where the second plot shows the error profile, the estimated leading coefficient $c(t)$ of the global error $e(t,h)=c(t)h+O(h^2)$ over time. The rapidly falling gray line is the error bound, safely below the actual error. The left plot of the actual solutions against the backdrop of a much more precise numerical solution clearly shows the linear convergence of the Euler method.
Best Answer
It is a misprint, it should be $$ e^{L(1-0)} $$ which is an upper bound for $$ \left(e^{L(1-0)}-1\right). $$