Floating point and global error in Euler Method

Question

Inspired by this answer, I tried to understand when floating point errors come into visibility and to check it also comparing the plot of the numerical solution with Explicit Euler with the analytical one for this simple problem.

\begin{cases} y'(x)= L \sin(x) \\ y(0)= 0 \end{cases}

where $L$ is a number. The Solution of the problem is $y(x)=L(1-cos(x))$. Notice that $$|L(\sin(x)-\sin(y)| \leq L |\cos(c)||x-y| \leq L|x-y|$$

so that the Lipschitz constant of the r.h.s is $L$.

---

Solving this ODE with Explicit Euler, we have a global error that is of the order $O(e^{Lx}(h+\mu))$ where $\mu$ is the unit roundoff.

Trying to use to @LutzLehmann's argument, I'd say that floating point errors come into visibility as $e^{Lx} \mu \approx 1$, i.e. $x \approx \log(\mu)/L$. I'd expect that for such an $x$, the numerical solution will deviate from the correct one.

To see this, I tested for $L=10$. It turns out that from about $x \approx 3.67$ I should start seeing some small problems, but in my plot it seems that everything is fine also for "extremely" large times, as may be seen in this plot with $T=40$ as final time.

My question: Why can't I see a significative difference starting from that point? What's really happening?

The simple Python code is the following:

import numpy as np
import matplotlib.pyplot as plt


L = 10
def Fun(t,x):
    return L*np.sin(t)
   
def Sol(t):
    return L*(1-np.cos(t))

tf=40.0
t0=0
k = 0.01
ts = np.int(tf/k)
t=np.linspace(t0,tf,ts+1)
y=np.zeros((ts+1)) #trapz
y0=0
y[0]=y0


for i in range (0,ts):
    y[i+1] = y[i] + k*Fun(t[i],y[i])

plt.plot(t,y,'o',label='Explicit Euler')
plt.plot(t,Sol(t),'-',label='Solution')
plt.title("Integration up to T="+str(tf))
plt.show()
plt.legend()

Answer 1

In a one-step method $$ y_{n+1}=y_n+h\Phi_f(x_n,y_n,h) $$ one gets a truncation error for the exact solution $$ y(x_{n+1})=y(x_n)+h\Phi_f(x_n,y(x_n),h)+h^{p+1}\tau(x_n) $$ For the error propagation of $e_n=y_n-y(x_n)$ this gives approximately $$ e_{n+1}=e_n+h\partial_y\Phi_f(x_n,y(x_n),h)e_n-h^{p+1}\tau(x_n)+... $$ The error is then of order $p$, $e_n=c(x_n)h^p+...$, and interpreting the difference equation as discretized differential equation one gets backwards $$ c'(x)=\partial_yf(x,y(x))c(x)-\tau(x). $$ In magnitude estimates of the error one assumes that $f$ around the exact solution has a Lipschitz constant $L$ so that then $$ \|c'(x)\|\le L\|c(x)\|+\|\tau\|_\infty\implies \|c(x)\|\le \frac{e^{L|x-x_0|}-1}L\|\tau\|_\infty $$ This however might be a very pessimistic bound, for instance in the scalar case where $\partial_yf(x,y)$ is negative. In general the magnification of the error follows the magnification of the value, so that the relative error does not increase much. To get observable differences one probably needs a rapidly oscillating system, which requires state dimension larger than 2.

---

There are two points where floating point errors enter the error calculation. First during the evaluation of $\Phi_f$, and then summing up the many small steps to the numerical solution. The first error has the size proportional to $\mu h$ for a single step, so that globally it will be proportional to $\mu$, mainly independent of the step size. The second error source is proportional to $\mu$ per step, so $\mu/h$ globally. However, this source is easy to eliminate or reduce to an error of size $\mu$ independent of $h$ using compensated addition.

---

In total it looks like in general the method step errors dominate the floating point errors, so that it will be difficult to separate the floating point error from the background of the method errors.

The obvious solution is to take a problem where the exact solution has zero truncation error, so that any deviation from the exact solution only starts with floating point errors. This means that the exact solution has to be a polynomial with a degree not larger than the order of the method. That is, take some polynomial $p$ and construct the IVP $$ y'(x)=Ly(x)+[p'(x)-Lp(x)], ~~~ y(0)=p(0). $$ This is the type of problems used in your previous question and as result in the comment discussion above. The method of manufactured solutions (MMS) can be extended to more general terms, like $$ y'(x)=L\sin(y(x))+p'(x)-L\sin(p(x)), ~~~ y(0)=p(0), $$ or even higher order constructs, $y''=F(x,y,y')+p''-F(x,p,p')$ etc.

---

One other method to see the difference between truncation and floating point errors is to compute the coefficients $c(x)$ of the leading error terms. This works well in the working range, and will deviate from that when the floating point errors start to dominate. Below the 5th order method of the DoPri 45 method with fixed steps is used with another non-polynomial MMS problem. The plot is of $(y_h(x)-p(x))/h^5$ and gives without compensated summation

and with compensated summation.

Note that the error for the smallest step size is bounded by $0.006·(0.0025)^5=5.859375·10^{-16}$, which is about what is maximally "physically" possible in 64bit floating point.

In the time loop below activating the commented line restores the naive uncompensated summation.

def DoPri45integratefixed(f, t, x0):
    N = len(t)-1;
    x = np.asarray((N+1)*[x0]);
    # update accumulates the errors of the summation
    update = 0*x[0];
    for k in range(N):
        h = t[k+1]-t[k]
        v4, v5 = DoPri45Step(f,t[k],x[k],h)
        #update = 0*x[0];
        update += h*v5
        x[k+1] = x[k] + update
        # remove the step at the precision of the solution function
        update -= (x[k+1]-x[k])
    return x