4
$\begingroup$

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?

enter image description here 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()
$\endgroup$
8
  • $\begingroup$ The Lipschitz constant relative $y$ is zero, as the right side does not depend on $y$. Try perhaps with $y'=F(x,y)=-Ly+\sin(x)$. $\endgroup$ – Lutz Lehmann Jan 7 at 14:57
  • $\begingroup$ @LutzLehmann I'll try this out now! What a mistake I did $\endgroup$ – Vefhug Jan 7 at 15:40
  • $\begingroup$ The solution of your equation, under the conditon $y(0)=0$ is $y(x)=\frac{1}{L^2+1} (e^{-Lx}+L \sin(x) - \cos(x))$, and the lipschitz constant of $F(x,y)$ is exactly $L$. However, for $L=10$ and step $h=0.05$ I don't see any strange behaviour starting from point $x \approx 3.67$. Also for larger time, everything seems fine @LutzLehmann What am I missing? $\endgroup$ – Vefhug Jan 7 at 16:19
  • $\begingroup$ Yes, and any deviation from the periodic part of the solution gets reduced exponentially, which also gives no striking error picture. In autonomous systems that somewhat stay bounded the error growth often happens as phase error, that is, with a time dilation the error becomes much smaller. Then of course you also have chaotic systems like for the Lorenz attractor where such a correction is not possible. $\endgroup$ – Lutz Lehmann Jan 7 at 16:26
  • $\begingroup$ I'm sorry but the global error is in any case $O(e^{Lx}(h+μ))$, so how is it possible that it decays? I'm missing this point @LutzLehmann Ho can I see that "floating point errors do not come into visibility" here? $\endgroup$ – Vefhug Jan 7 at 16:45
3
$\begingroup$

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

enter image description here

and with compensated summation.

enter image description here

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
$\endgroup$
6
  • $\begingroup$ I'm gonna read your answer carefully! Just two things: i) You plotted $e(x)/h^5$, which correspond to the leading coefficient $c(x)$. Then you noticed that for the smallest time step $h=0.0025$, the term $e(x)/h^5$ is bounded by roughly $0.006$, so to get the bound on the error $e(x)$ you just multiply for $h^5$ and you get $$0.006·(0.0025)^5=5.859375·10^{-16}$$, which tell us that the error is of the order of $\mu$. ii) You're seeing the difference between floating point and truncation, just by showing the compensated vs uncompensated plots. Is this correct? @LutzLehmann $\endgroup$ – Vefhug Jan 9 at 13:53
  • $\begingroup$ The deviation behavior for $h=0.0025$ in both cases is due to floating point errors, the error coefficient resulting from truncation converges visually up to that point. The difference between the graphs is just that one source of floating-point errors was radically reduced. $\endgroup$ – Lutz Lehmann Jan 9 at 13:55
  • $\begingroup$ Okay, I see the difference now, thanks. Also, in this example $p(x)$ (the analytical solution) is a poly of degree $p$ (let's say $p=4$) so that the truncation error is $0$ (because the method has order $5$): so when you plot those graphs, they must be the result of the floating point noise, right? @LutzLehmann $\endgroup$ – Vefhug Jan 9 at 14:08
  • $\begingroup$ No, that example had the non-polynomial $p(x)=\sin x$ with the second order ODE $y''(x)=p''(x)+\sin(p(x))-\sin(y(x))$, $y(0)=p(0)$, $y'(0)=p'(0)$. Thus the presence of truncation errors in an oscillating pattern. $\endgroup$ – Lutz Lehmann Jan 9 at 14:20
  • $\begingroup$ So when we can observe a deviation behaviour from the truncation error pattern, then the source is the floating point noise. Another way to see the floating point noise would be to choose a problem with exact solution with $0$ truncation error ( i.e. a polynomial with degree less than order $p$) and show the error graph. In your example, even if we can see that for $h=0.0025$ we have a "deviation", it does not affect the goodness of the approximation, right @LutzLehmann? $\endgroup$ – Vefhug Jan 9 at 14:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.