Test some numerical methods on a stiff scalar problem¶
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 numerical 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 8.3 of the textbook to illustrate the need for L-stable methods.
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$.
This exact solution is explored in the notebook ScalarStiffness.html.
%matplotlib inline
from pylab import *
from ipywidgets import interact, IntSlider, FloatSlider
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
Forward Euler¶
def forward_euler(nsteps, eta, lam):
from scipy.optimize import fsolve
f = lambda u,t: lam*(u-cos(t)) - sin(t)
t = linspace(0, tfinal, nsteps+1)
dt = t[1] - t[0]
U = empty((nsteps+1)) # array for computed solution
U[0] = eta
for n in range(nsteps):
Un = U[n]
tn = t[n]
U[n+1] = Un + dt*f(Un, tn)
figure(figsize=(10,4))
axis([-1,tfinal, -3,3])
tfine,ufine = utrue(0., eta, lam)
plot(tfine, ufine, 'k', label='true solution')
plot(t, U, 'bo-', label='Forward Euler')
legend(loc='lower right')
title('Forward Euler method with $k = %g, \lambda = %g, \quad k\lambda = %g$' \
% (dt,lam, dt*lam))
Forward Euler is absolutely stable only if $-2 \leq k\lambda \leq 0$. Here's a case where it is barely stable. It looks ok if the initial data is $u(0) = 1$ so there is no rapid transient in the true solution, and because the time step is small enough that the one-step errors introduced are not too large:
forward_euler(nsteps=95, eta=1, lam=-10)
But if we solve the equation with the initial condition $u(0) = 0$, giving a rapid transient in the true solution, then the near-instability is apparent:
forward_euler(nsteps=95, eta=0, lam=-10)
With a slightly larger timestep it goes unstable and the numerical solution grows exponentially:
forward_euler(nsteps=93, eta=0, lam=-10)
Trapezoidal method¶
Next we implement the Trapezoidal method on this same problem.
def trapezoid(nsteps, eta, lam):
from scipy.optimize import fsolve
f = lambda u,t: lam*(u-cos(t)) - sin(t)
t = linspace(0, tfinal, nsteps+1)
dt = t[1] - t[0]
U = empty((nsteps+1)) # array for computed solution
U[0] = eta
for n in range(nsteps):
Un = U[n]
tn = t[n]
tnp = t[n+1]
g = lambda u: u - Un - 0.5*dt*f(Un,tn) - 0.5*dt*f(u,tnp)
Unp = fsolve(g, Un)
U[n+1] = Unp
figure(figsize=(10,4))
axis([-1,tfinal, -3,3])
tfine,ufine = utrue(0., eta, lam)
plot(tfine, ufine, 'k', label='true solution')
plot(t, U, 'bo-', label='Trapezoid')
legend(loc='lower right')
title('Trapezoid method with $k = %g, \lambda = %g, \quad k\lambda = %g$' \
% (dt,lam, dt*lam))
This method does much better with the parameters we used above for Forward Euler:
trapezoid(nsteps=93, eta=0, lam=-10)
The Trapezoid method is A-stable and remains stable even when $k\lambda \rightarrow -\infty$. But it is not L-stable, and $R(z) \rightarrow -1$ as $|z| \rightarrow \infty$. Hence if $|k\lambda|$ is very large we expect that a rapid transient will not be damped but rather will oscillate:
trapezoid(nsteps=80, eta=0, lam=-1000)
Backward Euler or the second-order accurate TR-BDF2 method would both do better in this case, since they are L-stable.
def tr_bdf2(nsteps, eta, lam):
from scipy.optimize import fsolve
f = lambda u,t: lam*(u-cos(t)) - sin(t)
t = linspace(0, tfinal, nsteps+1)
dt = t[1] - t[0]
U = empty((nsteps+1)) # array for computed solution
U[0] = eta
for n in range(nsteps):
Un = U[n]
tn = t[n]
tnp = t[n+1]
tnph = 0.5*(t[n] + t[n+1]) # half time step
g1 = lambda u: u - Un - 0.25*dt*f(Un,tn) - 0.25*dt*f(u,tnph)
Ustar = fsolve(g1, Un)
g2 = lambda u: u - 4/3*Ustar + 1/3*Un - 1/3*dt*f(u,tnp)
Unp = fsolve(g2, Ustar)
U[n+1] = Unp
figure(figsize=(10,4))
axis([-1,tfinal, -3,3])
tfine,ufine = utrue(0., eta, lam)
plot(tfine, ufine, 'k', label='true solution')
plot(t, U, 'bo-', label='TR-BDF2')
legend(loc='lower right')
title('TR-BDF2 method with $k = %g, \lambda = %g, \quad k\lambda = %g$' \
% (dt,lam, dt*lam))
tr_bdf2(nsteps=93, eta=0, lam=-10)
tr_bdf2(nsteps=80, eta=0, lam=-1000)
