Step size and stability of Euler forward method

Question

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?

Answer 1

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.