Step size and stability of Euler forward method

I'm trying to calculate the maximum step size that provides stability for the following nonlinear IVP using the Euler forward method:

$$u'(t) = -200tu(t)^2,\qquad u_0 = 1, \qquad t\in [0,3]$$,

with the analytical solution being $$u = (1+100t^2)^{-1}$$ (see figure below). From linear stability analysis, that is by solving $$u' = \lambda u, \qquad u_0 = 1, \qquad \lambda <0$$, one can show that the region of stability for the Euler forward is $$|1+\lambda h|\le1$$, where $$h$$ is the step size. So the step size must be $$h<2/|\lambda|$$.
I realize that this holds only for linear problems, however when solving a general IVP (including nonlinear terms):

$$u'(t) = f(t,u(t)), \qquad u(t_0)=u_0$$

calculating the eigenvalues $$\lambda_i(t)$$ of $$A=f_{u_i}(t,u(t))$$, where $$f_{u_i}$$ is the Jacobian matrix, should also lead to reasonable step sizes when using the lowest value of all $$\lambda_i$$. That is if the $$u(t)$$ is asymptotically stable.
So the next step would be to calculate the Jacobian matrix (in this case with just one entry):
$$f_{u_i}=-\frac{400t}{1+100t^2}=\lambda$$. The infimum of $$\lambda$$ on the given interval is $$\lambda=-20$$ and thus $$h = 2/20$$. The figure below shows $$\lambda(t)$$.

However the step size seems to be too large as the model blows up. Stability seems to be somewhere around $$h=2/29$$. So I have two Questions:
1: What did I do wrong while calculating the step size or is the example just pathological and if so why?
2: How can I actually calculate the maximum step size providing stability?

• Keep in mind there is no guarantee of stability for forward Euler method. As a result, the step-size that you calculated based on your stability analysis is just an approximation and for nonlinear ODEs it might be really terrible approximation. In fact, foward Euler method has only conditional stability even for linear ODEs and you should not expect to estimate step-size precisely even for linear ODEs. Nov 5 '19 at 15:49
• Just to make sure, we're on same page: by conditional stability you mean, that it is not sufficient for the stability of Euler forward that $\lambda <0$, but it also needs to be within the region of stability (since forward Euler is not A-stable)? Nov 5 '19 at 16:05
• Yes, it means for plenty of cases if you continue integrating long enough, finally you will see instability. Nov 5 '19 at 16:07

Using the exact solution to estimate $$\lambda$$ may not give a time step that works for the numerical scheme. The local estimate of $$\lambda$$ is $$\lambda = 400 t u$$ If you try this step $$h = \min(2/20, 2/(400 t u+10^{-12}))$$ it is stable, but of course accuracy is poor. Here is a code

import numpy as np
import matplotlib.pyplot as plt

t = 0.0
u = 1.0
tdata, udata = [], []
tdata.append(t); udata.append(u)
while t<3.0:
dt = np.min([2.0/20.0, 2.0/(400*t*u+1.0e-12)])
rhs = - 200.0*t*u**2
u = u + dt*rhs
t += dt
print("dt, t =", dt, t)
tdata.append(t); udata.append(u)

tdata = np.array(tdata); udata = np.array(udata)
uexact = 1.0/(1.0 + 100*tdata**2)
plt.plot(tdata,udata,tdata,uexact)
plt.legend(('Forward Euler','Exact'))
plt.xlabel('t'); plt.ylabel('u')
plt.show()


In first step, step size is $$h=0.1$$. After first step, the numerical solution is still $$u=1$$ and $$\lambda = 400*0.1*1 = 40$$ so you need to use step size $$h \le 2/40 = 0.05$$ in second step, which is less than $$2/20=0.1$$.

There is another issue with this particular problem. Starting at $$u=1$$, the solution should remain in $$[0,1]$$. If it becomes negative for some reason, then forward euler will take it to $$-\infty$$, as happens for $$h=2/20$$. The ODE itself has this behaviour. So for this problem, if exact solution is supposed to be positive, time step should be chosen so that numerical solution remains positive. $$u_{n+1} = u_n - 200 t_n u_n^2 h$$ We have $$u_n \ge 0$$ $$\implies$$ $$u_{n+1} \ge 0$$ iff $$h \le \frac{1}{200 t_n u_n}$$ and this condition also agrees with the restriction from linear stability analysis.

• So my mistake was to look at the analytical derivative instead of looking at the numerical one. Thanks Nov 6 '19 at 10:46
• I updated my answer. Nov 6 '19 at 11:29