computer scienceeuler's methodMATLABnumerical methodsrounding error
In particular, if Euler's method is implemented on a computer, what's the minimum step size that can be used before rounding errors cause the Euler approximations to become completely unreliable?
I presume it's when step size reaches the machine epsilon? E.g. if machine epsilon is e-16, then once step size is roughly e-16, the Euler approximations are unreliable.
Yes, the process of computation stays the same when $h<0$; you get a solution to the left of the initial point. An adjustment in error bound formulas may be needed, since they are written under the assumption $h>0$. So, replace $\text{error}\le Ch$ with $\text{error}\le C|h|$.
As hardmath pointed out in a comment, using $h<0$ in Euler's method is different from using Backward (aka implicit) Euler method.
Yes, there is a "generic type" limit on the size of the time step. It is related to the stability of the method. Stability means that if the solution or one of its components (or a linear combination of them) converges in the exact solution, then this should also happen in the numerical solution. This does not guarantee that the solution is good in an accuracy sense, only that it is not fundamentally wrong.
In the Euler method, to avoid the more radically wrong qualitative behavior, you need to satisfy $$ |1+\lambda h|<1 $$ for all eigenvalues of the Jacobian of the order 1 formulation with negative real part. A weakened, i.e., not sufficient, condition for not completely useless step sizes is $$ Lh<2 \qquad (\text{better use }Lh<1.5) $$ where $L$ is a Lipschitz constant of the first order ODE system.
The first order formulation of this ODE has a constant matrix with the same characteristic polynomial as the 2nd order ODE and thus eigenvalues $λ_1=-1$ and $λ_2=-2$, resulting in the condition $h<1$.
If one uses a semi-factored first order system with $v=y'+y$, then \begin{align} y'&=v-y\\ v'&=-2v \end{align} has Lipschitz constant $L=2$ in the maximum norm. Only $h=0.2$ is not disqualified from the stability condition in both versions.
Interpreting that $f=\delta$ is the Dirac delta "function", with the usual assumption of $y(t)=0$ for $t<0$, this means that $y'$ has a jump from $0$ to $1$ at $t=0$. Thus one gets the initial conditions (in the right-sided limits) $y(0^+)=0$, $y'(0^+)=1$ for the solution for $t>0$. The exact solution there is $$y(t)=e^{-t}-e^{-2t}.$$
It is easy to perform the numerical experiment using any scripting or programming language. Here I set this up in a python script to compare the exact solution to the numerical solutions with the given step sizes
# set up the equation
def odesys(u): y,v = u; return np.array([ v, -3*v-2*y]);
# plot the exact solution
t = np.arange(0,10,0.1);
plt.plot(t, np.exp(-t)-np.exp(-2*t), color='lightgray', lw=4);
# Euler method for the 4 step sizes
for h in [ 2, 1.5, 1, 0.2]:
t = np.arange(0,10,h);
# initial conditions
u = [ np.array([0.0,1.0]) ];
# Perform Euler steps according to the time array
for s in t[:-1]: u.append( u[-1]+h*odesys(u[-1]) );
# plot the solution
u = np.asarray(u);
plt.plot(t, u[:,0], '-o', ms=2+2*h);
plt.ylim(-2.5,2.5); plt.grid(); plt.show()
As can be seen, the solution for $h=1$ stays barely bounded, while $h=0.2$ slightly overshoots the exact solution, but correctly falls to zero.
Remember that the Euler methods are order $1$. As example, to get an error of magnitude $10^{-4}$ over a time interval of length $1$ would require $10000$ steps/function evaluations.
With an order $2$ method, Heun or improved Euler, one would need $100$ steps with $200$ function evaluations.
With an order $4$ method like the classical Runge-Kutta, one would need for the same error magnitude $10$ steps with $40$ function evaluations.
Thus the advantage of higher order methods. And why the question about the Euler method is purely academic.
Best Answer
It depends on the problem. Basically, you have two sources of errors in each step, the rounding errors that are some multiple of the machine constant $\mu=2\cdot10^{-16}$ (scaled by the scale of the functions involved) and the local truncation error, which depends on derivatives of the ODE function (thus also involving the function scale) and the square $h^2$ of the step size. You want the truncation error to dominate the rounding errors, which in the simplest case demands $h^2>\mu$ or $h>\sqrt\mu\sim 10^{-8}$. This changes if the ODE function is in some sense "strongly curved".