% SIMIODE-TeacherStudent-Template.tex and simiode.cls for LaTeX2e
% Version 22 November 2016.
%
% This template is to be used with the  class file "simiode.cls".
%
% The latest version of this class and template files
% can be found at
%
%       <http://www.SIMIODE.org> in the Author section.
%
% A "Quick Start Guide" appears below.
%
% This file must be processed with version LaTeX2e or higher.
% Appreciation to Mark Frantz, IUPUI; John Thoo, Yuba College; and Eric Sullivan, Carroll College.
%
\documentclass{simiode}
%
% With this version of the tex file and the associated SIMIODE.cls file you can produce BOTH a Student Version and a Teacher Version of your modeling scenario or technique narrative.
%
% Student Version contains the title STUDENT VERSION along with the NAME OF SCENARIO, with no abstract, no keywords, no tags, no author identifiers, and no COMMENTS section. The Student Version contains only the STATEMENT of the modeling scenario only with appropriate figures, data, instructions, and perhaps a bibliography appropriate for student resourcing, but nothing in the bibliography which would "give away" the nature of the modeling scenario, etc.
%
%The Teacher Version contains the material in the Student Version, namely the title NAME OF SCENARIO, now headed with   TEACHER VERSION, but also full author  identification, abstract, key words, tag words, and, following the STATEMENT, a section of COMMENTS which can be discussion of presentation  ideas, "solutions," techniques, etc.
%
% Authors can prepare a full Teacher Version using the notions immediately below  to create a pdf for review, but in review please eliminate all information about the author as we use a double-blind referee system. Such information will be restored when accepted and published in SIMIODE at www.simiode.org.
%
%***********   Toggle the following line to turn the teacher version on (1) or off (0). ***********
%
\def\TeacherEdition{0} % Toggle this line's value.
\newcommand{\teacher}[1]{\ifnum\TeacherEdition=1 #1 \fi}
%
% The \teacher{ } command wraps any text that will appear only in the teacher version of
% the document. See this file for examples. This has been done for you in the document below.
% Rename the pdf after you compile each final version.
%
\usepackage{amsmath,amssymb,amsthm}
\usepackage{color}
\usepackage{graphicx}
%\usepackage{hyperref}
\usepackage{url}
\usepackage{exercise}
\def\sv{\left [ \begin{array}{c}}
\def\ev{\end{array} \right ]}
\def\smtwo{\left [ \begin{array}{cc}}
\def\smthree{\left [ \begin{array}{ccc}}
\def\smfour{\left [ \begin{array}{cccc}}
\def\smfive{\left [ \begin{array}{ccccc}}
\def\em{\end{array} \right ]}
%
\renewcommand{\baselinestretch}{1.3}
%
% NOTE regarding macros %%
%
% Any macros defined for your paper should be contained in
% the top matter. Likewise, any environment definitions such
% as \newtheorem or \newenvironment should also go in the
% top matter.
%
% Do not \input your macros as separate files within your
% final version because more files make it harder for the
% editors and users to keep track of and/or modify and
% append your paper.
%
% Top matter is above and the beginning of the document below.
%
\theoremstyle{definition}
\newtheorem{example}{Example}[section]
\begin{document}

\newpage
%
% You only need to change the NAME OF SCENARIO in the next line and delete the footnote.
% See above for instructions on toggling from Teacher Version to Student Version.
%
\title{\ifnum\TeacherEdition=1 TEACHER VERSION
%\footnote{DELETE THIS FOOTNOTE BEFORE SUBMISSION.  In the \TeX \, source code there is a toggle/option to produce a Teacher Version or a Student Version by ``carving'' out with  $\backslash$teacher\{\} that material which should only appear in the Teacher version. See the source code.}
\else STUDENT VERSION \fi\\Stiff Differential Equations}



\markboth{Stiff Differential Equations}{Stiff Differential Equations}

\author{Kurt Bryan\\
Department of Mathematics\\
Rose-Hulman Institute of Technology\\
Terre Haute IN  47840 U.S.\\
kurt.bryan@rose-hulman.edu
}

\makeStitlePDFLaTex


\begin{abstract}

This material introduces the topic of ``stiffness'' for a system of ordinary differential equations (ODE's), through a series of examples.
Stiffness is a property that a system of ODE's
may posses that make it difficult to solve numerically with standard methods, and it is a surprisingly common phenomenon. Many standard numerical methods applied
 to stiff systems of ODE's will either yield an unstable iteration that grows without limit, or just grind to a halt. Stiffness is frequently a mathematical manifestation of a physical problem that has two or more very different scales for the independent variable. There are numerical ODE methods well-suited to handling stiff systems and
 we illustrate these methods with examples. The material is suitable for an introductory ODE course in which students have encountered systems of ODE's (at least, linear systems)
 and basic numerical methods, for example, Euler's method.

\end{abstract}


\section{INTRODUCTION}

Most real-world systems of ordinary differential equations (ODE's) must be solved numerically, either because no closed-form solution exists,
or the solution is too unwieldy to manipulate symbolically.
``Stiffness'' is a particular property that some systems of ODE's have that makes them challenging to
solve numerically. Stiffness often reflects an essential disparity of scale in a physical system---one facet of a physical system may evolve in time
on a relatively slow scale, e.g., seconds, while another facet of the same system evolves much more rapidly, e.g., microseconds. It might be the slow
scale/large time evolution of the system that interests us, but the rapidly changing aspect of the system's behavior demands the attention of the numerical method
and makes long-term simulations computationally expensive. Stiff systems are surprisingly common, but rarely talked about in introductory ODE
or numerical analysis courses.

We consider here a variety of straightforward examples that illustrate the essential idea behind stiffness.
The goal here is not a treatise on the theory of numerical methods for ODE's and stiff systems, but rather an intuitive and practical guide
to what stiff ODE's are and some ideas on how to handle them. There are also some physical examples and exercises to illustrate how stiff systems come about.

\section{NUMERICAL METHODS FOR ODE'S REVIEW}

Consider a first order initial value problem
\begin{equation}
\label{ode1}
y'(t) = f(t,y(t))
\end{equation}
for an unknown scalar function $y(t)$,  with initial condition $y(t_0)=a$. You can consider $t$ as time. It is frequently
the case that we can't solve (\ref{ode1}) explicitly and so must resort to numerical approximation.
Recall that Euler's method is one of the simplest schemes for this computation, and is really just repeated
tangent line extrapolation. In the examples that follow we'll often use Euler's method to illustrate, but the conclusions
generally hold for any ``explicit one-step'' numerical method (defined below).

In Euler's method we choose a step size ``$h$'' with $h>0$ for the independent variable $t$ and then construct a sequence
of approximations $y_k \approx y(t_k)$ where $t_k=t_0+kh$ with $k=1,2,\ldots$, as follows. If $h$ is sufficiently small the approximation
\begin{equation}
\label{fordiff}
y'(t) \approx \frac{y(t+h)-y(t)}{h}
\end{equation}
should be accurate (in the limit $h\rightarrow 0$ the right side converges to $y'(t)$, assuming $y$ is differentiable). The approximation (\ref{fordiff}) can be rearranged to
$$
y(t+h) \approx y(t) + h y'(t)
$$
and if $y(t)$ satisfies (\ref{ode1}) we then have
\begin{equation}
\label{findiff}
y(t+h) \approx y(t) + h f(t,y(t)).
\end{equation}
Equation (\ref{findiff}) tells us how to extrapolate the solution forward from time $t$ to time $t+h$, at least approximately.

We can use (\ref{findiff}) to approximate $y(t)$ for $t_0\leq t\leq T$ for some final time $T=t_0+Mh$ with the following algorithm, Euler's Method:
\begin{enumerate}
\item Set $y_0=a$ and counter $k=0$. Set $t_k=t_0+kh$.
\item Use (\ref{findiff}) to generate $y_{k+1}$, an approximation to $y(t_{k+1})$, as
\begin{equation}
\label{euler}
y_{k+1} = y_k + hf(t_k,y_k).
\end{equation}
\item If $t_k=T$ (equivalently, $k=M$) stop and return $y_0,y_1,\ldots,y_M$. Otherwise, increment $k\leftarrow k+1$
and return to step 2.
\end{enumerate}
The smaller the step size $h$, the better this algorithm ``tracks'' the true solution.

\begin{example}
Consider the classic ``exponential decay'' ODE
\begin{equation}
\label{decay}
y' = \lambda y
\end{equation}
where $\lambda<0$, with initial condition $y(0)=1$. The true solution is $y(t)=e^{\lambda t}$.

For $\lambda=-1$ we have $y(t)=e^{-t}$. If we approximate $y(t)$ with Euler's method on the interval $0\leq t\leq 5$ with step size $h=0.5$, the situation
is as shown in the left panel in Figure \ref{fig1}. Also shown is the direction field for the ODE (\ref{decay}). When $h$ is smaller, the Euler procedure
tracks the direction field or solution more accurately, as shown in the right panel of Figure \ref{fig1}.
\begin{figure}[ht]
\centerline{\includegraphics[width=0.4\textwidth]{example1.jpg}\ \includegraphics[width=0.4\textwidth]{example1b.jpg}}
\caption{Euler approximation to $y'=-y$ solution, step sizes $0.5$ (left) and $0.1$ (right).}
\label{fig1}
\end{figure}
\end{example}

But now consider what happens if we solve (\ref{decay}) with $\lambda=-3$ and the same step size $h=0.5$, as shown in the left panel in Figure \ref{fig2}. The Euler approximation
is very poor, and qualitatively incorrect---the true solution is always positive, but the approximation is actually negative at times. When $\lambda=5$ the situation is even worse,
as illustrated in the right panel of Figure \ref{fig2},
for although the true solution decays smoothly to zero, Euler's method is unstable---the iterates grow without bound!
\begin{figure}[ht]
\centerline{\includegraphics[width=0.4\textwidth]{example2.jpg}\ \includegraphics[width=0.4\textwidth]{example3.jpg}}
\caption{Euler approximation to $y'=-\lambda y$ solution, step size $0.5$, $\lambda=3$ (left) and $\lambda=5$ (right).}
\label{fig2}
\end{figure}

If you examine Euler's method geometrically in relation to the direction field, it's easy to see what goes wrong. Although the true solution decays smoothly to zero for any $\lambda<0$, when $\lambda$ is ``very'' negative the solution changes rapidly, with a rather steep rate of decay. Even if the current iterate $y_k$ is exact, Euler steps that are too large end up extrapolating to a wildly incorrect value (very negative). For this ODE, a large step size in Euler's method yields an iteration that isn't simply inaccurate, but actually unstable.
\subsection{Exercises}
\begin{enumerate}
\item For the ODE (\ref{decay}) with $y(0)=1$, show that Euler's method produces iterates that satisfy $y_{k+1}=(1+\lambda h)y_k$, and that (given that $y_0=1$)
\begin{equation}
\label{eqn1}
y_k = (1+\lambda h)^k.
\end{equation}

\item Based on (\ref{eqn1}), for a fixed $\lambda<0$ show that the condition $|1+\lambda h|<1$ requires that $h<2/|\lambda|$ if we want
$y_k\rightarrow 0$ as $k\rightarrow \infty$ (so the iterates decay, just like the solution). What happens if this condition is not met ($h$ is too large?) What condition on $h$ assures that the $y_k$ iterates actually remain positive?

\item The Improved Euler Method replaces the formula (\ref{euler}) in Euler's Method with the two-step procedure
\begin{align}
u & = y_k + hf(t_k,y_k)\nonumber\\
y_{k+1} & = y_k + \frac{h}{2}\left (f(t_k,y_k)+f(t_{k+1},u)\right )\label{impeuler}
\end{align}
and is generally more accurate than Euler's Method.

Solve the ODE (\ref{decay}) with $\lambda=-1$ and $y(0)=1$ using the Improved Euler method, with step size $h=0.5$ out to time $t=5$. Look at the iterates $y_0,\ldots,y_{10}$. Do they remain positive? Decay to zero? Are they accurate? Repeat with $\lambda=-3$ and $\lambda=-5$ (step size $h=0.5$). Does the Improved Euler method suffer from the same
stability problems as Euler's method?

\end{enumerate}

The moral here is that for many numerical methods and ODE's, step sizes that are too large result in more than just inaccurate iterates. The
numerical solution actually becomes unstable and ``blows up'', even if the true solution decays. This is starting to get at what stiffness is about, but it's not the whole story. Stiffness is really a property possessed by a system of differential equations. Let's look at some more examples, now involving systems of ODE's.

\section{STIFF SYSTEMS OF ODE'S}

A system of ODE's can be written in a form similar to (\ref{ode1}), namely
\begin{equation}
\label{ode2}
\mathbf{y}'(t) = \mathbf{f}(t,\mathbf{y}(t))
\end{equation}
where $\mathbf{y}(t)=\langle y_1(t),\ldots,y_n(t)\rangle$ is a vector-valued function taking values in $\mathbb{R}^n$, whose components are the functions of interest. The function $\mathbf{f}(t,\mathbf{y}(t))$
is a vector-valued function taking values in $\mathbb{R}^n$ that indicates how to compute the derivative of each function $y_k(t)$ in terms of $t$ and $\mathbf{y}(t)$. Euler's Method for such a system takes precisely the same form as the scalar version and is simply repeated extrapolation on each component of $\mathbf{y}(t)$.

\begin{example}
Let's start by considering a system of ODE's of the form
\begin{align}
y_1'(t) & = \lambda_1 y_1(t)\label{x1eqn}\\
y_2'(t) & = \lambda_2 y_2(t)\label{x2eqn}
\end{align}
for some constants $\lambda_1,\lambda_2<0$. In the framework of (\ref{ode2}) we have $\mathbf{f}(t,\mathbf{y}) = \langle\lambda_1 y_1,\lambda_2 y_2\rangle $.
Of course, this ``system'' is in fact decoupled---each equation can be solved independently of the other. Nonetheless, it will serve to illustrate the essential
issue with stiff systems. The exact solution is $y_1(t) = y_1(0)e^{\lambda_1 t}$ and $y_2(t) = y_2(0)e^{\lambda_2 t}$.

Let's use Euler's Method to solve the system starting at time $t_0=0$ with $y_1(0)=0.5$ and $y_2(0)=0.5$. Let $h$ be the step size, $t_k=kh$, and let's use $y_1^k$ and $y_2^k$ to denote the iterates
produced by Euler's Method, as approximations to $y_1(t_k)$ and $y_2(t_k)$. Euler's Method (\ref{ode2}) for the system (\ref{x1eqn})-(\ref{x2eqn}) takes the form
\begin{align}
y_1^{k+1} & = (1+\lambda_1 h)y_1^k\label{eul1}\\
y_2^{k+1} & = (1+\lambda_2 h)y_2^k\label{eul2}.
\end{align}
It's easy to see then that $y_1^k = y_1(0) (1+\lambda_1 h)^k$ and $y_1^k = y_2(0) (1+\lambda_2 h)^k$ from Exercise 1 above.

Suppose that $\lambda_1$ and $\lambda_2$ are both negative, and recall Exercise 2 above. Both true solutions decay to zero, but $|1+\lambda_j h|>1$
if $h>2/|\lambda_j|$, and so $y_j^k$ (here $j=1$ or $j=2$) will grow without bound as the iteration counter $k$ increases, which does not reflect a qualitatively (or quantitatively!) accurate
solution; Euler's method will be unstable. This means the maximum step size we can take is dictated by the smaller of the quantities $2/|\lambda_1|,2/|\lambda_2|$, that is, the larger of $|\lambda_1|,|\lambda_2|$. If, for example, $\lambda_2<\lambda_1<0$, it's $\lambda_2$ that dictates the largest permissible step size. But it may be the case that we're interested
in the evolution of the system related to $\lambda_1$, which occurs on a longer time scale, long after the $e^{\lambda_2 t}$ part of the solution has decayed. Yet we're forced to take tiny steps with size less than $2/|\lambda_2|$ to track this solution.

To illustrate, suppose $\lambda_1=-1$ and $\lambda_2=-100$. Then $y_1(t)=y_1(0) e^{-t}$, while $y_2(t)=y_2(0) e^{-100 t}$. Here $y_2(t)$ decays much more rapidly than
$y_1(t)$ and very quickly will be insignificant in magnitude. If we are interested in tracking the decay of $y_1(t)$, however, we have to continue to iterate
(\ref{eul1}) and (\ref{eul2}) in a stable manner, which requires $h<2/100$, even though (\ref{eul1}) itself would only require $h<2$. This means taking unnecessarily many more iterations to get out to a large enough time to see the evolution of $\mathbf{y}(t)$. This system is (moderately)``stiff.''
\end{example}

\begin{example}
Of course in the above example it would make sense to solve (\ref{x1eqn}) and (\ref{x2eqn}) independently, each with an appropriate step size, but most systems are not decoupled
like this one. Consider instead a slightly more complicated (but still linear) system of the form
\begin{equation}
\label{stiffsys}
\sv y'_1(t)\\y'_2(t) \ev = \smtwo -56 & 55\\
44 & -45 \em \sv y_1(t)\\y_2(t) \ev.
\end{equation}
The eigenvalues and eigenvectors of the matrix on the right in (\ref{stiffsys}) are $\lambda_1=-1,\lambda_2=-100$ with corresponding eigenvectors $\mathbf{v}_1=<\langle 1,1\rangle $ and $\mathbf{v}_2=\langle1,-0.8\rangle $.
As such, the general solution is
\begin{equation}
\label{gensol}
\mathbf{y}(t) = c_1 e^{-t} \sv 1\\1\ev + c_2 e^{-100t}\sv 1\\-0.8\ev.
\end{equation}
Take initial conditions $y_1(0)=0.5$ and $y_2(0)=0.5$; these are slightly ``rigged'' initial conditions, to illustrate the vital point, but we'll consider
something more typical below. We find that $c_1=0.5$ and $c_2=0.0$. The solution is
$$
\mathbf{y}(t) = 0.5 e^{-t} \sv 1\\1\ev.
$$
The $e^{-100t}$ portion isn't even present in the true solution (and if it were, it would decay toward zero long before the solution above would).

It seems that a step size of $h=0.1$ or smaller in Euler's method would do a pretty good job of tracking the solution, but this is not the case!
A plot of $y_1(t)$ (which equals $y_2(t)$) is shown in Figure \ref{fig5} in the left panel for $0\leq t\leq 5$, along with the Euler iterates.
\begin{figure}[ht]
\centerline{\includegraphics[width=0.4\textwidth]{example4.jpg}\ \includegraphics[width=0.4\textwidth]{example5.jpg}}
\caption{True solution and Euler's method with $h=0.022$ (left) and RK4 method with $h=0.03$ (right).}
\label{fig5}
\end{figure}
However, this step size results in disaster---the iterates quickly overflow floating point arithmetic. A bit of experimentation shows that no step size larger than $h=0.02$ is stable
in the long run. Switching to a more sophisticated numerical solver doesn't help much. In the right panel a Runge-Kutta 4th order method still falls apart when $h=0.03$. Again, as with the decoupled
system (\ref{x1eqn})-(\ref{x2eqn}), it's the fastest decaying exponentials in the solution that dictate the largest stable step size. These fast decaying exponentials correspond to the largest magnitude eigenvalues in a linear system.

We can get a little more geometric insight into what's going wrong in this example by looking at the direction field for the system and a more typical solution. Suppose that $y_1(0)=0.55$ and $y_2(0)=0.57$. The solution to (\ref{stiffsys}) is, to good a approximation,
\begin{equation}
\label{trutru}
\mathbf{y}(t) \approx 0.5611 e^{-t}\sv 1\\1\ev - 0.0111e^{-t/100}\sv 1\\-0.8 \ev.
\end{equation}
The direction field for the system (\ref{stiffsys}) is shown in Figure \ref{fig12}, along with a parametric plot of the true solution $\mathbf{y}(t)$ (the black curve).
Bear in mind that for aesthetic reasons, the arrows that delineate the direction field are not drawn to scale---they are much longer than indicated in Figure \ref{fig12}. In particular, the arrows that lie very far off the diagonal line $y_2=y_1$ are rather long and point almost orthogonally toward this line, a geometric manifestation of the
fact that solutions like (\ref{gensol}) decay rapidly to a multiple of the vector $\langle 1,1\rangle $, and then slowly toward the origin.

Of course the true solution in this case does exactly this, starting at $\langle 0.55,0.57\rangle $ and quickly decaying to the line $y_2=y_1$. At all points the true solution curve follows the direction field (that's what it means to be a solution!)
However, the Euler iterates (shown in blue/dashed) with step size $h=0.022$ do not behave like this. The first Euler step starts at $\langle0.55,0.57\rangle$; at this point
the vector that dictates the direction the solution should go is $\langle 0.55,-1.45\rangle$, but a step of $h=0.022$ puts the next iterate at $\langle0.55,0.57\rangle+0.022\langle 0.55,1.45\rangle \approx \langle 0.5621,0.5381\rangle$. This overshoots the line $y_2=y_1$ and actually ends up farther off this line than the starting point $\langle0.55,0.57\rangle$! This makes matters even worse, for the direction field has an even larger magnitude here, and the next iterate overshoots the other way. The process is repeated, with obvious effect. In every case, however, Euler's method is following the direction field toward the true solution, but due to the large magnitude of the direction field (that stems from the large eigenvalue of $-100$) Euler overshoots
and oscillates in an unbounded manner.
\begin{figure}[ht]
\centerline{\includegraphics[width=0.8\textwidth]{example12.jpg}}
\caption{Euler iterates versus true solution and the direction field.}
\label{fig12}
\end{figure}
\end{example}

The heart of the problem in this example is the large disparity in the magnitude of the eigenvalues. The portion of the solution (\ref{gensol}) corresponding to
the $-1$ eigenvalue or $e^{-t}$ term dominates in most cases, at least after a very short time, and the solution decays smoothly and gradually to zero. But the portion of the solution corresponding to the eigenvalue of $-100$ is what limits the size of the steps we can take, no matter how small that portion of the solution might be. In many stiff systems one might
encounter eigenvalues that vary over a much larger range, e.g., five or six orders of magnitude, or more! We are then forced to take time steps perhaps one million
times smaller than the actual rate of change of the solution dictates.  See Exercises 2 and 3 below.

This is the essence of stiffness, which does not have a precise or universally accepted definition. Nonetheless, the general idea is
\begin{quote}
A system of ODE's is \textit{stiff} if the step size required to maintain stability in a numerical ODE solver is small in relation to the scale on which the solution changes
with respect to the independent variable.
\end{quote}
This can occur in a linear system $\mathbf{x}'=\mathbf{A}\mathbf{x}$ if the eigenvalues of $\mathbf{A}$ (real or complex) are of widely varying magnitude. See Exercise 3 at the end for an example with complex eigenvalues. In a nonlinear system the phenomena can be even more complicated; the solution may be stiff in some intervals and not in others.

Modern numerical ODE solvers also incorporate ``adaptive stepsizing.'' That is, the step size $h$ is adjusted at each time step in order to control the (estimated) error of the solver at each step. When things are going poorly, $h$ is decreased. This can combat stiffness to some extent, but in practice what ultimately happens is that $h$ is reduced almost to $0$ and the solver grinds to a halt. Another approach to overcoming stiffness is needed.

\section{THE IMPLICIT EULER METHOD}

If a system is stiff, what can be done?  One common solution is to use an ``implicit'' ODE solver. Euler's method and most other numerical methods encountered in an introductory ODE course are ``explicit, one-step'' methods. What this means is that $\mathbf{y}^{k+1}$ is computed from the previous iterate $\mathbf{y}^{k}$, the time $t_k$, and the step size $h$ (and the ODE itself, of course), in the form
\begin{equation}
\label{explicit}
\mathbf{y}^{k+1} = \mathbf{y}^{k} + h\phi(t_k,\mathbf{y}^{k},h)
\end{equation}
where $\phi(t,\mathbf{y},h)$ is a function of the indicated variables. For example, for an ODE (\ref{ode1}) with Euler's method we have $\phi(t,\mathbf{y},h) = \mathbf{f}(t,\mathbf{y})$. For the Improved Euler method we have (using (\ref{impeuler})) $\phi(t,\mathbf{y},h) = \frac{1}{2}(\mathbf{f}(t,\mathbf{y}) + \mathbf{f}(t+h,\mathbf{y} +h\mathbf{f}(t,\mathbf{y})))$. The ``explicit'' refers to the fact that $\mathbf{y}^{k+1}$ is given explicitly in (\ref{explicit}). The ``one-step'' means that only $\mathbf{y}^k$ (and not previous iterates $\mathbf{y}^{k-1},\mathbf{y}^{k-2},$ etc.) are involved in the right side of (\ref{explicit}).

Explicit methods are convenient because we can ``march'' the iterates forward in time with a minimum of fuss. Their drawback, for stiff systems, is that the step size
$h$ has to be small enough to maintain stability of the method, even if the solution doesn't change rapidly. An alternative is to use an ``implicit'' method. These frequently
remain stable for large (even arbitrarily large) step sizes.

The simplest example of such a method is the \textit{implicit Euler} method. The implicit version of
Euler's method for a scalar ODE (\ref{ode1}) replaces the approximation (\ref{fordiff}) with a so-called ``backwards difference''
\begin{equation}
\label{backdiff}
y'(t)\approx \frac{y(t)-y(t-h)}{h}
\end{equation}
(again, $h>0$ here). The right side of (\ref{backdiff}) converges to $y'(t)$ as $h\rightarrow 0$, so if $h$ is small enough the approximation of
(\ref{backdiff}) should be accurate. With $t=t_k$ (\ref{backdiff}) becomes $y'(t_k) \approx (y(t_k)-y(t_{k-1}))/h$, which can be rearranged to $y(t_k) \approx y(t_{k-1}) + hy'(t_k)$. Finally,
if $y$ satisfies (\ref{ode1}) then we have $y(t_k)\approx y(t_{k-1})+hf(t_k,y(t_k))$. If $y_k$ denotes our approximation to $y(t_k)$ this becomes a recipe for how
to compute $y_k$, as a solution to $y_k = y_{k-1} + h f(t_k,y_k)$. Actually, we can shift indices $k\rightarrow k+1$ and write the implicit Euler method as
\begin{equation}
\label{impliciteul}
y_{k+1} = y_{k} + h f(t_{k+1},y_{k+1}),
\end{equation}
a prescription for computing $y_{k+1}$ from $y_k$ by solving (\ref{impliciteul}). If $f$ is sufficiently simple we might be able to solve for $y_{k+1}$ in closed-form, but in most cases we actually have to use a numerical method (e.g., Newton's method) just to obtain $y_{k+1}$. Thus the implicit Euler's method probably requires more work at each iteration. The payback is that the method remains stable for large step sizes $h$. Equation (\ref{impliciteul}) generalizes to systems in an obvious way.

\begin{example}
To illustrate, let's consider the ODE (\ref{decay}) again, with $y(0)=1$, this time with the implicit Euler method. In this case equation (\ref{impliciteul}) yields
$
y_{k+1} = y_k - h \lambda y_{k+1}.
$
Because $f$ is so simple here, we can actually solve for $y_{k+1}$ in closed form as
$$
y_{k+1} = \frac{y_{k}}{1+\lambda h}.
$$
With $y_0=1$ its easy to see that in general
\begin{equation}
\label{decimp}
y_k = \frac{1}{(1+\lambda h)^k}.
\end{equation}
Contrast this to (\ref{eqn1}). In the case that $\lambda=-5$ and $h$ is reasonably small, the implicit Euler method does a good job of approximating the true solution. But
 even if $h$ is large---say $h=0.5$---the implicit Euler method remains stable, and even positive, albeit a bit inaccurate. See Figure \ref{fig6}.
\begin{figure}[ht]
\centerline{\includegraphics[width=0.4\textwidth]{example6.jpg}}
\caption{True solution and implicit Euler's method with $h=0.5$.}
\label{fig6}
\end{figure}
In fact (\ref{decimp}) makes it clear that any step size $h>0$ will yield such an approximation.
\end{example}

\begin{example}
To further illustrate, let's apply the implicit Euler Method to the system (\ref{stiffsys}). As with standard Euler's method, the implicit approach extends to a system of $n$ ODE's
(\ref{ode2}) in $n$ unknowns and leads to
\begin{equation}
\label{impliciteulsys}
\mathbf{y}^{k+1}= \mathbf{y}^k + h\mathbf{f}(t_{k+1},\mathbf{y}^{k+1}),
\end{equation}
a system of $n$ (possibly nonlinear) equations in $n$ unknowns, the components of the vector $\mathbf{y}^{k+1}$. For the system $\mathbf{y}'=\mathbf{A}\mathbf{y}$ of (\ref{stiffsys}) this means we must solve
\begin{equation}
\label{impeul2}
\mathbf{y}^{k+1}= \mathbf{y}^k + h\mathbf{A}\mathbf{y}^{k+1}
\end{equation}
to compute $\mathbf{y}^{k+1}$ from $\mathbf{y}^k$ at each iteration. This is equivalent to solving $(\mathbf{I}-h\mathbf{A})\mathbf{y}^{k+1}=\mathbf{y}^k$ for $\mathbf{y}^{k+1}$, a linear system of equations that can be handled by Gaussian elimination or any linear solver. This is a fortunate consequence of the fact that our ODE's here are linear. In general the resulting equations for $\mathbf{y}^{k+1}$ could be nonlinear and something like Newton's method might be needed.

%, so we have $\mathbf{y}^{k+1} = (\mathbf{I}-h\mathbf{A})^{-1}\mathbf{y}^k$.
In fact in this case since $\mathbf{A}$ has eigenvalues $-1$ and $-100$, the matrix $\mathbf{I}-h\mathbf{A}$
has nonzero positive eigenvalues $1+h$ and $1+100h$ if $h>0$. Thus $\mathbf{I}-h\mathbf{A}$ must be invertible and so (\ref{impeul2}) will have a unique solution for $\mathbf{y}^{k+1}$.

Figure \ref{fig7} shows the result of this iteration, a plot of the true solution $y_1(t)=e^{-t}$ and the implicit Euler iterates for $y_1$ with step size $h=0.1$. Contrast this with the results of Figure \ref{fig5}. We can get away with a much larger step size here and still maintain stability. An excessively large step size will still result in a poor approximation though, just one that doesn't blow up!
\begin{figure}[ht]
\centerline{\includegraphics[width=0.4\textwidth]{example7.jpg}}
\caption{True solution and implicit Euler's method with $h=0.1$.}
\label{fig7}
\end{figure}
\end{example}

\section*{SUMMARY}


There are many other implicit methods for systems of ODE's. These methods often take the form
\begin{equation}
\label{imp}
\psi(\mathbf{y}^{k+1},\mathbf{y}^k,t_k,t_{k+1}) = 0,
\end{equation}
for some function $\psi$,
at least in a one-step version. Here $\mathbf{y}^{k+1}$ must be solved for, usually in some non-trivial way. An example is the \textit{Trapezoidal Method},
which takes the form
\begin{equation}
\label{trap}
\mathbf{y}^{k+1} = \mathbf{y}^k + \frac{h}{2}(\mathbf{f}(t_k,\mathbf{y}^k)+\mathbf{f}(t_{k+1},\mathbf{y}^{k+1})).
\end{equation}
The use of such a method is generally required if the equations are stiff and the solution is to be computed stably over an interval of
any significant length. Most software packages for scientific computing (e.g., Maple, Mathematica, Matlab) have built-in solvers specifically designed for stiff systems of ODE's.
See \cite{butcher} for much more information on numerical methods for ODE's, stiff and otherwise.

\section{EXERCISES}

\subsection{An RC Circuit}


Consider the RC circuit in Figure \ref{fig3}, components and current directions as labeled. Let $q_1(t)$ denote the charge on capacitor $C_1$ and $q_2(t)$ the charge on capacitor $C_2$.
\begin{figure}[ht]
\centerline{\includegraphics[width=0.3\textwidth]{circuit.jpeg}}
\caption{Two loop RC circuit.}
\label{fig3}
\end{figure}

Suppose the open switch is closed at time $t=0$ (and the capacitors $C_1$ and $C_2$ have some initial charge). From Kirchhoff's Voltage Law applied to the loop involving $R_1$ and $C_1$ we have
\begin{equation}
\label{rc1}
R_1 i_1(t) + q_1(t)/C_1 = 0.
\end{equation}
For the loop involving $C_1,C_2$, and $R_2$ we have
\begin{equation}
\label{rc2}
R_2i_3(t) + q_2(t)/C_2 - q_1(t)/C_1 = 0.
\end{equation}
At the top junction between $R_1,R_2,$ and $C_1$ we must have
\begin{equation}
\label{rc3}
i_1 - i_2 - i_3 = 0.
\end{equation}
Finally, we also have
\begin{align}
q_1' & = i_2\label{rc4}\\
q_2' & = i_3.\label{rc5}
\end{align}
Using (\ref{rc3}) to eliminate $i_1$ in (\ref{rc1})-(\ref{rc2}), then (\ref{rc4})-(\ref{rc5}) to replace $i_2,i_3$, and finally solving for $q'_1$ and $q'_2$ leads to
\begin{align}
q_1'(t) & = -\left (\frac{1}{R_2C_1}+\frac{1}{R_1C_2}\right )q_1(t) + \frac{1}{R_2C_2} q_2(t)\label{rcode1}\\
q_2'(t) & = \frac{1}{R_2C_1} q_1(t) - \frac{1}{R_2C_2} q_2(t)\label{rcode2}
\end{align}
a linear system of ODE's for functions $q_1(t),q_2(t)$.

Suppose $C_1=C_2=0.001$ farad, $R_1=1000$ ohms, and $R_2=0.1$ ohm. The ODE's (\ref{rcode1})-(\ref{rcode2}) become
\begin{align}
q_1'(t) & = -(1.0001\times 10^7)q_1(t) + (1.0\times 10^7)q_2(t)\label{rcode1b}\\
q_2'(t) & = (1.0\times 10^7)q_1(t) - (1.0\times 10^7)q_2(t).\label{rcode2b}
\end{align}
With initial condition $q_1(0)=0.000001$ and $q_2(0)=0.000001$ the true solution to (\ref{rcode1b})-(\ref{rcode2b}) is (to good approximation)
\begin{align}
q_1(t) & = (9.99975\times 10^{-7})e^{-499.99 t} + (2.5\times 10^{-11}) e^{(-2\times 10^7) t}\label{q1sol}\\
q_2(t) & = (1.0\times 10^{-6})e^{-499.99 t} - (2.5\times 10^{-11}) e^{(-2\times 10^7) t}.\label{q2sol}
\end{align}
Note that both of $q_1(t)$ and $q_2(t)$ contain a comparatively large amount of $e^{-499.99 t}$, a decaying exponential, and a very
small multiple of $e^{(-2\times 10^7) t}$, a much faster decaying exponential (this term decays $40,000$ times faster!).\\\ \\
{\bf RC Circuit Exercises}
\begin{enumerate}
\item Plot $q_1(t)$ and $q_2(t)$ on the range $0\leq t\leq 0.01$. (They should look pretty similar).

\item It seems that a step size of $h=0.0001$ would do a good job of tracking both of $q_1(t)$ and $q_2(t)$. Solve the system (\ref{rcode1b})-(\ref{rcode2b}) numerically
with Euler's method and $h=0.0001$ on the interval $0\leq t\leq 0.01$ (or as far as the numerics will go). What happens?

Decrease $h$ until the Euler iteration is stable over $0\leq t\leq 0.01$. How small does $h$ have to be?

\item Repeat the last problem but with the Implicit Euler method and step size $h=0.0001$. Compare to the true solution.

\item Try the Trapezoidal Method (\ref{trap}).  What step sizes yield stability?
\end{enumerate}

Note the disparity of time scales in the circuit's physical behavior. The $e^{(-2\times 10^7)t}$ stems from the capacitors discharging through the small
$R_2$ resistor, and they quickly assume an equal charge; this occurs on a time scale of $10^{-7}$ seconds. But both capacitors then discharge through the much larger $R_1$ resistor at a comparatively slow rate---this is the $e^{-499.99t}$ term in the solutions, and this discharges occurs on a scale of $10^{-3}$ seconds.
It may be this slower scale behavior that interests us, while the very rapid equalization of the capacitor charges that occurs when the switch is closed is of no interest.
Nonetheless, it is this short time-scale phenomena that dictates the very small step size Euler's method must take to remain stable. We are thus forced to take many very
small time steps in order to track the solution out to time $0.01$. This problem is not confined to Euler's Method, but shared by any explicit method.

\subsection{A Spring-Mass Problem}

Consider a simple spring-mass system consisting of two masses $m_1$ and $m_2$, as illustrated in Figure \ref{fig8}. Let's assume the masses are constrained to move horizontally and that the only forces acting on either mass are the springs themselves and friction (no gravity). The spring attaching $m_1$ to the wall has spring constant $k_1$, while the spring connecting $m_1$ to $m_2$ has constant $k_2$. Suppose both springs have natural length $1$.


\begin{figure}[ht]
\centerline{\includegraphics[width=0.4\textwidth]{springmass.jpg}}
\caption{Double spring-mass system.}
\label{fig8}
\end{figure}


Let $x_1(t)$ denote the displacement of the spring attaching $m_1$ to the wall from its natural length
(e.g., $x_1=0$ means the spring is at its natural length, neither stretched nor compressed) and $x_2(t)$ the same for the second spring.
Assume both springs obey Hooke's Law. The position of the first mass would then be $1+x_1(t)$ and the position of the second mass is $2+x_1(t)+x_2(t)$.
The velocity of the first mass is then $x'_1(t)$ and the velocity of the second mass is $x'_1(t)+x'_2(t)$.

The spring with constant $k_1$ exerts a force $-k_1x_1$ on $m_1$, while the spring with constant $k_2$ exerts a force $k_2(x_2-x_1)$. Suppose also that
$m_1$ experiences a frictional force proportional to its velocity, say $-c_1 x'_1$ for some constant $c_1>0$. The net force on $m_1$ is thus
$F_1 = -k_1x_1+k_2(x_2-x_1)-c_1 x'_1 = -(k_1+k_2)x_1 + k_2 x_2 -c_1 x'_1$. From Newton's Second Law, $F=ma$, we have
\begin{equation}
\label{spring1}
m_1 x_1'' = -(k_1+k_2)x_1 + k_2 x_2 -c_1 x'_1.
\end{equation}
The force exerted by the spring with constant $k_2$ on $m_2$ is $-k_2(x_2-x_1) = k_2 x_1 - k_2 x_2$. If the frictional force on $m_2$ is proportional
to the velocity of $m_2$ then this force is of the form $-c_2(x'_1 + x'_2)$ for some $c_2>0$, and so
\begin{equation}
\label{spring2}
m_2 x_2'' = k_2 x_1 - k_2 x_2 - c_2 (x'_1 + x'_2).
\end{equation}
With initial conditions $x_1(0)=a,x'_1(0)=v_1,x_2(0)=b,x'_2(0)=v_2$, equations (\ref{spring1})-(\ref{spring2}) form a couple pair of linear second order ODE's. There is a unique solution.

If we let $y_1=x_1,y_2=x'_1,y_3=x_2,$ and $y_4=x'_2$, the system (\ref{spring1})-(\ref{spring2}) can be formulated as a first order system of four differential equations in four
unknowns
\begin{align*}
y'_1 & = y_2\\
y'_2 & = -\frac{k_1+k_2}{m_1} y_1 - \frac{c_1}{m_1} y_2 + \frac{k_2}{m_1} y_3\\
y'_3 & = y_4\\
y'_4 & = \frac{k_2}{m_2} y_1 - \frac{c_2}{m_2}y_2 - \frac{k_2}{m_2} y_3 - \frac{c_2}{m_2} y_4
\end{align*}
or, more compactly, as $\mathbf{y}' = \mathbf{A}\mathbf{y}$ where
\begin{equation}
\label{Adef}
\mathbf{A} = \smfour 0 & 1 & 0 & 0\\
-\frac{k_1+k_2}{m_1} & -\frac{c_1}{m_1} & \frac{k_2}{m_1} & 0\\
0 & 0 & 0 & 1\\
\frac{k_2}{m_2} & -\frac{c_2}{m_2} & -\frac{k_2}{m_2} & -\frac{c_2}{m_2} \em.
\end{equation}
{\bf Spring-Mass Exercises}
\begin{enumerate}
\item Take $k_1=1,k_2=10^4,m_1=1,m_2=0.001,c_1=c_2=0.1$. Compute $\mathbf{A}$ and its eigenvalues (they are complex). What exponential functions govern the decay rates in this problem? What natural frequencies are present?

\item With the parameters as in (1) above and initial conditions $y_1(0)=0.1$ and $y_2(0)=y_3(0)=y_4(0)=0$.
it turns out that the true solution for $y_1(t)$ that governs the motion of the first mass is (to a few significant figures)
\begin{align*}
y_1(t) & \approx e^{-49.9 t}((9.98\times 10^{-5})\cos(3163.5 t) + (4.72\times 10^{-6})\sin(3163.5 t))\\
& + e^{-0.1499 t}(0.0999\cos(0.9882 t) + 0.0051\sin(0.9882 t)).
\end{align*}
Note the fast oscillating term is much smaller and decays much faster than the more slowly oscillating term.
Plot this for $0\leq t\leq 50$. Repeat on a smaller interval, e.g., $0\leq t\leq 5$. The graphs of $y_2(t),y_3(t),$ and $y_4(t)$ look fairly similar.
It seems a step size of $h=0.01$ would track the solution well.

\item Solve the system with Euler's method and $h=0.01$ out to time $t=5$, or however far it will go. Is it stable?

\item Repeat with the Implicit Euler method. Is the numerical solution stable? Accurate?

\item Repeat the last question with the Trapezoidal Method (\ref{trap}).
\end{enumerate}
\teacher{
\section*{NOTES FOR TEACHERS}

The notion of stiffness is rarely part of an introductory ODE course; it might be covered in an undergraduate numerical analysis course, but even this
isn't common. Nonetheless, stiffness is a common problem in many physical ODE models.
I first encountered a stiff system of ODE's when I started working in industry after college. Unfortunately, I'd never seen the concept at that time
and was left wondering why the fancy Runge-Kutta ODE solver that I'd written in a college numerical analysis class kept blowing up.

This material makes a great extra project in an introductory ODE course for students who like numerics and programming. It's nice way to give
students some extra practice with some of the basic numerical ODE methods in the setting of systems of ODE's (whose numeric treatment is often neglected
as we focus on scalar equations).

\section*{EXERCISE SOLUTIONS}

\noindent{\bf Section 2}
\begin{enumerate}
\item The DE (\ref{decay}) has $f(t,y) = \lambda y$, so Euler's method is easily seen to become $y_{k+1}=(1+\lambda h)y_k$, and $y_0=1$. Of course then
$y_1 = (1+\lambda h) y_0$, $y_2=(1+\lambda h)y_1 = (1+\lambda h)y_0^2,\ldots$. An obvious induction gives $y_k = (1+\lambda h)^k$.

\item The condition $|1+\lambda h|<1$ is the same as $-1<1+\lambda h < 1$, or $0<2+\lambda h < 2$. Since we assume $\lambda<0$ and $h>0$ the right inequality
is always true. Then $0<2+\lambda h$ can be rearranged to $h<-2/\lambda$ or $h<2/|\lambda|$. The argument is reversible too, so if $h>2/|\lambda|$ we find $|1+\lambda h|>1$
and the iterates will grow!

\item Improved Euler here
with $\lambda=-1$ gives $y(5)\approx 0.009095$ (true value $e^{-5}\approx 0.006738$.) The iterates remain positive and decay to zero (since $1+\lambda h=0.5$ here). When $\lambda=-3$ with $h=0.5$ we actually get the same iterates as for $\lambda=-1$ (sheer coincidence!) But when $\lambda=-5$ we obtain $y(5)\approx 128.4$; the iterates remain positive, but grow without bound (as $y_k=(13/8)^k$, it turns out).
\end{enumerate}

\noindent{\bf RC Circuit Solutions}
\begin{enumerate}
\item A plot of the exact solutions for $0\leq t\leq 0.01$ is shown in Figure \ref{fig4} on the left; $q_1(t)$ and $q_2(t)$ are almost identical.
\begin{figure}[ht]
\centerline{\includegraphics[width=0.4\textwidth]{rcresp1.jpg}\ \includegraphics[width=0.4\textwidth]{rcexer1.jpg}}
\caption{Charge on capacitors over time, $q_1$ and $q_2$ exact (left) and implicit Euler approximation with $h=0.0001$ (right).}
\label{fig4}
\end{figure}
Both functions decay smoothly to zero.

\item The numerical solution with $h=0.0001$ exceeds $10^{302}$ in magnitude by $t=0.01$! We need to decrease $h$ to about $h=10^{-7}$ before the iterates decay to zero.

\item Implicit Euler with $h=0.0001$ does a good job, as shown in the right panel of Figure \ref{fig4}.

\item Because the system is linear with $f(t,\mathbf{y}) = \mathbf{A}\mathbf{y}$
where
$$
\mathbf{A} = \smtwo -1.0001\times 10^{-7} & 1.0\times 10^7\\
1.0\times 10^7 & -1.0\times 10^7 \em
$$
we find that (\ref{trap}) can be written as
$$
(\mathbf{I}-\frac{h}{2}\mathbf{A})\mathbf{y}^{k+1} = (\mathbf{I}+\frac{h}{2}\mathbf{A})\mathbf{y}^{k}
$$
and this is a linear system, easy to solve for $\mathbf{y}^{k+1}$. The choice $h=0.0001$ does an excellent job. In fact, any choice for $h$ remains stable, though the
approximation is poor if $h$ is large.
\end{enumerate}

\noindent{\bf A Spring-Mass Problem Solutions}
\begin{enumerate}
\item The matrix is given by
$$
\mathbf{A} = \smfour 0 & 1 & 0 & 0\\
-10001 & -1/10 & 10000 & 0\\
0 & 0 & 0 & 1\\
10^7 & -100 & -10^7 & -100\em
$$
and has approximate eigenvalues $-49.9\pm 3163.5i$, $-0.14985 \pm 0.9882i$.
The two natural decay rates in this system are governed by $e^{-49.4t}$ and $e^{-0.149t}$, and natural frequencies $3163.5$ and $0.988$ radians per second ($503.5$ and $0.157$ Hz, respectively) are present.

The plots of $y_1(t)$ for $0\leq t\leq 50$ and $0\leq t\leq 5$ are shown in Figure \ref{fig13}.
\begin{figure}[ht]
\centerline{\includegraphics[width=0.4\textwidth]{mass1a.jpg}\ \includegraphics[width=0.4\textwidth]{mass1b.jpg}}
\caption{Displacement of mass $m_1$ for $0\leq t\leq 50$ (left) and $0\leq t\leq 5$ (right).}
\label{fig13}
\end{figure}

\item Euler's method with $h=0.01$ blows up rapidly, hitting $10^{144}$ for the magnitude of $y_1(t)$ by the time $t=1$.

\item Implicit Euler with $h=0.01$ does an excellent job, as shown in the left panel of Figure \ref{fig15}. The true solution $y_1(t)$ has value $8.369\times 10^{-3}$, while
implicit Euler estimates $7.829\times 10^{-3}$.
\begin{figure}[ht]
\centerline{\includegraphics[width=0.4\textwidth]{massimpeul.jpg}\ \includegraphics[width=0.4\textwidth]{masstrap.jpg}}
\caption{Implicit Euler solution (left) and Trapezoidal rule (right), both with $h=0.01$.}
\label{fig15}
\end{figure}
\item The Trapezoidal Method with $h=0.01$ does an even better job, as shown in the right panel of Figure \ref{fig15}. We obtain estimate $y_1(5)\approx 8.402\times 10^{-3}$.
%\section*{ACKNOWLEDGEMENTS}
}

\end{enumerate}
\begin{thebibliography}{99}

\bibitem{butcher} Butcher, J.C., 2016. \textit{Numerical Methods for Ordinary Differential Equations, Third Edition}. New York; John Wiley \& Sons, Inc.
\end{thebibliography}

\end{document}
