2
$\begingroup$

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).enter image description here 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)$.
enter image description here
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?

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Alone Programmer Nov 5 '19 at 15:49
  • $\begingroup$ 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)? $\endgroup$ – J. Fregin Nov 5 '19 at 16:05
  • $\begingroup$ Yes, it means for plenty of cases if you continue integrating long enough, finally you will see instability. $\endgroup$ – Alone Programmer Nov 5 '19 at 16:07
1
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ So my mistake was to look at the analytical derivative instead of looking at the numerical one. Thanks $\endgroup$ – J. Fregin Nov 6 '19 at 10:46
  • $\begingroup$ I updated my answer. $\endgroup$ – cfdlab Nov 6 '19 at 11:29

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.