A scalar problem illustrating stiffness

AMath 586, Spring Quarter 2019 at the University of Washington. For other notebooks, see Index.html or the Index of all notebooks on Github.

This notebook shows the solution to the ODE

$$ u'(t) = \lambda(u(t)-\cos(t)) - \sin(t), \qquad u(t_0) = \eta $$

as $t_0, \eta,$ and $\lambda$ are varied. This ODE is used in Examples 7.1 and 8.1 of the textbook.

The exact solution is

$$ u(t) = \cos(t) + \exp(\lambda(t-t_0))(\eta - \cos(t_0)). $$

Note in particular that if $u(0)=1$ is specified then the solution is just $u(t)=\cos(t)$ and if $\lambda=0$ then other solutions remain parallel to this, since in this case the ODE $u'(t)=-\sin(t)$ has solution $u(t) = \cos(t) + (\eta - \cos(t_0))$. Note that in this case $f(u,t)$ is independent of $u$ and the ODE can be solved by simple integration.

If $\lambda<0$ then solutions decay towards this "slow solution". This equation exhibits stiffness when $\lambda$ is very negative and we wish to compute the solution for over times that are long relative to $-1/\lambda$.

In [1]:
%matplotlib inline
In [2]:
from pylab import *
from ipywidgets import interact, IntSlider, FloatSlider
In [3]:
tfinal = 6*pi
def utrue(t0, eta, lam):
    t = linspace(t0,tfinal,1000)
    u = cos(t) + exp(lam*(t-t0))*(eta - cos(t0))
    return t,u
In [4]:
def plot_soln(t0, eta, lam):
    t,u = utrue(t0, eta, lam)
    figure(figsize=(10,4))
    axis([0,tfinal, -3,3])
    tequil = linspace(0, tfinal, 1000)
    plot(tequil, cos(tequil), 'k', label='solution with $u(0)=1$')
    plot([t[0]], [u[0]], 'ro', label='initial data')
    plot(t, u, 'r', label='solution with specified data')
    legend(loc='upper right')
    grid(True)
    title('With $\lambda = %.1f,$' % lam)

Here's a typical solution, when $u(3) = 2$ is specified as initial data, with $\lambda=-2$:

In [5]:
plot_soln(3, 2, -2)

The interact below allows you to vary the parameters. Note that changing $\eta$ with $\lambda=0$ gives solution curves that are "parallel" to the slow solution, and that decreasing $\lambda$ gives decay towards this solution.

In [6]:
interact(plot_soln, t0=IntSlider(min=0,max=10,value=0, description='$t_0$'),
         eta=FloatSlider(min=-2.5,max=2.5,step=0.5,value=1., description='$\eta$'),
         lam=FloatSlider(min=-20,max=0,step=0.5,value=0,description='$\lambda$'));
In [ ]: