Floating point and global error in Euler Method
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 <PERSON>'s argument, I'd say that floating point errors come into visibility as $e^{Lx} \mu \approx 1$, i.e. 521
$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() 509