%%
%% Automatically generated file from DocOnce source
%% (https://github.com/hplgit/doconce/)
%%

% #define PREAMBLE

% #ifdef PREAMBLE


\documentclass[%
oneside,                 % oneside: electronic viewing, twoside: printing
final,                   % draft: marks overfull hboxes, figures with paths
10pt]{article}

\listfiles               %  print all files needed to compile this document

\usepackage{relsize,makeidx,color,setspace,amsmath,amsfonts,amssymb}
\usepackage[table]{xcolor}
\usepackage{bm,ltablex,microtype}

\usepackage[pdftex]{graphicx}

% Packages for typesetting blocks of computer code
\usepackage{fancyvrb,framed,moreverb}

% Define colors
\definecolor{orange}{cmyk}{0,0.4,0.8,0.2}
\definecolor{tucorange}{rgb}{1.0,0.64,0}
\definecolor{darkorange}{rgb}{.71,0.21,0.01}
\definecolor{darkgreen}{rgb}{.12,.54,.11}
\definecolor{myteal}{rgb}{.26, .44, .56}
\definecolor{gray}{gray}{0.45}
\definecolor{mediumgray}{gray}{.8}
\definecolor{lightgray}{gray}{.95}
\definecolor{brown}{rgb}{0.54,0.27,0.07}
\definecolor{purple}{rgb}{0.5,0.0,0.5}
\definecolor{darkgray}{gray}{0.25}
\definecolor{darkblue}{rgb}{0,0.08,0.45}
\definecolor{darkblue2}{rgb}{0,0,0.8}
\definecolor{lightred}{rgb}{1.0,0.39,0.28}
\definecolor{lightgreen}{rgb}{0.48,0.99,0.0}
\definecolor{lightblue}{rgb}{0.53,0.81,0.92}
\definecolor{lightblue2}{rgb}{0.3,0.3,1.0}
\definecolor{lightpurple}{rgb}{0.87,0.63,0.87}
\definecolor{lightcyan}{rgb}{0.5,1.0,0.83}

\colorlet{comment_green}{green!50!black}
\colorlet{string_red}{red!60!black}
\colorlet{keyword_pink}{magenta!70!black}
\colorlet{indendifier_green}{green!70!white}

% Backgrounds for code
\definecolor{cbg_gray}{rgb}{.95, .95, .95}
\definecolor{bar_gray}{rgb}{.92, .92, .92}

\definecolor{cbg_yellowgray}{rgb}{.95, .95, .85}
\definecolor{bar_yellowgray}{rgb}{.95, .95, .65}

\colorlet{cbg_yellow2}{yellow!10}
\colorlet{bar_yellow2}{yellow!20}

\definecolor{cbg_yellow1}{rgb}{.98, .98, 0.8}
\definecolor{bar_yellow1}{rgb}{.98, .98, 0.4}

\definecolor{cbg_red1}{rgb}{1, 0.85, 0.85}
\definecolor{bar_red1}{rgb}{1, 0.75, 0.85}

\definecolor{cbg_blue1}{rgb}{0.87843, 0.95686, 1.0}
\definecolor{bar_blue1}{rgb}{0.7,     0.95686, 1}

\usepackage{minted}
\usemintedstyle{default}

\usepackage[T1]{fontenc}
%\usepackage[latin1]{inputenc}
\usepackage{ucs}
\usepackage[utf8x]{inputenc}

% Set helvetica as the default font family:
\RequirePackage{helvet}
\renewcommand\familydefault{phv}

\usepackage{lmodern}         % Latin Modern fonts derived from Computer Modern

% Hyperlinks in PDF:

\usepackage{hyperref}
\hypersetup{
    breaklinks=true,
    colorlinks=true,
    linkcolor=seccolor,
    urlcolor=seccolor,
    citecolor=black,
    filecolor=black,
    %filecolor=blue,
    pdfmenubar=true,
    pdftoolbar=true,
    bookmarksdepth=3   % Uncomment (and tweak) for PDF bookmarks with more levels than the TOC
    }
%\hyperbaseurl{}   % hyperlinks are relative to this root

\setcounter{tocdepth}{2}  % levels in table of contents

% Tricks for having figures close to where they are defined:
% 1. define less restrictive rules for where to put figures
\setcounter{topnumber}{2}
\setcounter{bottomnumber}{2}
\setcounter{totalnumber}{4}
\renewcommand{\topfraction}{0.95}
\renewcommand{\bottomfraction}{0.95}
\renewcommand{\textfraction}{0}
\renewcommand{\floatpagefraction}{0.75}
% floatpagefraction must always be less than topfraction!
% 2. ensure all figures are flushed before next section
\usepackage[section]{placeins}
% 3. enable begin{figure}[H] (often leads to ugly pagebreaks)
%\usepackage{float}\restylefloat{figure}

% --- fancyhdr package for fancy headers ---
\usepackage{fancyhdr}
\fancyhf{} % sets both header and footer to nothing
\renewcommand{\headrulewidth}{1pt}
% section name to the left (L) and page number to the right (R)
% on even (E) pages, the other way around on odd pages
% (switch twoside to onside in documentclass to just have odd pages)
\fancyhead[LE,RO]{\rightmark} % section
\fancyhead[RE,LO]{\thepage}
% Ensure copyright on titlepage (article style) and chapter pages (book style)
\fancypagestyle{plain}{
  \fancyhf{}
  \fancyfoot[C]{{\footnotesize \copyright\ 2016, Hans Petter Langtangen. Released under CC Attribution 4.0 license}}
%  \renewcommand{\footrulewidth}{0mm}
  \renewcommand{\headrulewidth}{0mm}
}
% Ensure copyright on titlepages with \thispagestyle{empty}
\fancypagestyle{empty}{
  \fancyhf{}
  \fancyfoot[C]{{\footnotesize \copyright\ 2016, Hans Petter Langtangen. Released under CC Attribution 4.0 license}}
  \renewcommand{\footrulewidth}{0mm}
  \renewcommand{\headrulewidth}{0mm}
}

\pagestyle{fancy}


\usepackage{framed,wrapfig}

% --- begin definitions of admonition environments ---

% Admonition style "yellowicon" has colored background, yellow icons, and no farme
% Admon "notice"
\colorlet{yellowicon_notice_background}{yellow!5}
% \fboxsep sets the space between the text and the box
\newenvironment{noticeshaded}
{\def\FrameCommand{\fboxsep=3mm\colorbox{yellowicon_notice_background}}
 \MakeFramed {\advance\hsize-\width \FrameRestore}}{\endMakeFramed}

\newenvironment{notice_yellowiconadmon}[1][Notice]{
\begin{noticeshaded}

\begin{wrapfigure}{l}{0.07\textwidth}
\vspace{-13pt}
\includegraphics[width=0.07\textwidth]{latex_figs/small_yellow_notice}
\end{wrapfigure} \textbf{#1}\par
\nobreak\ignorespaces
}
{
\end{noticeshaded}
}

% Admonition style "yellowicon" has colored background, yellow icons, and no farme
% Admon "summary"
\colorlet{yellowicon_summary_background}{yellow!5}
% \fboxsep sets the space between the text and the box
\newenvironment{summaryshaded}
{\def\FrameCommand{\fboxsep=3mm\colorbox{yellowicon_summary_background}}
 \MakeFramed {\advance\hsize-\width \FrameRestore}}{\endMakeFramed}

\newenvironment{summary_yellowiconadmon}[1][Summary]{
\begin{summaryshaded}

\begin{wrapfigure}{l}{0.07\textwidth}
\vspace{-13pt}
\includegraphics[width=0.07\textwidth]{latex_figs/small_yellow_summary}
\end{wrapfigure} \textbf{#1}\par
\nobreak\ignorespaces
}
{
\end{summaryshaded}
}

% Admonition style "yellowicon" has colored background, yellow icons, and no farme
% Admon "warning"
\colorlet{yellowicon_warning_background}{yellow!5}
% \fboxsep sets the space between the text and the box
\newenvironment{warningshaded}
{\def\FrameCommand{\fboxsep=3mm\colorbox{yellowicon_warning_background}}
 \MakeFramed {\advance\hsize-\width \FrameRestore}}{\endMakeFramed}

\newenvironment{warning_yellowiconadmon}[1][Warning]{
\begin{warningshaded}

\begin{wrapfigure}{l}{0.07\textwidth}
\vspace{-13pt}
\includegraphics[width=0.07\textwidth]{latex_figs/small_yellow_warning}
\end{wrapfigure} \textbf{#1}\par
\nobreak\ignorespaces
}
{
\end{warningshaded}
}

% Admonition style "yellowicon" has colored background, yellow icons, and no farme
% Admon "question"
\colorlet{yellowicon_question_background}{yellow!5}
% \fboxsep sets the space between the text and the box
\newenvironment{questionshaded}
{\def\FrameCommand{\fboxsep=3mm\colorbox{yellowicon_question_background}}
 \MakeFramed {\advance\hsize-\width \FrameRestore}}{\endMakeFramed}

\newenvironment{question_yellowiconadmon}[1][Question]{
\begin{questionshaded}

\begin{wrapfigure}{l}{0.07\textwidth}
\vspace{-13pt}
\includegraphics[width=0.07\textwidth]{latex_figs/small_yellow_question}
\end{wrapfigure} \textbf{#1}\par
\nobreak\ignorespaces
}
{
\end{questionshaded}
}

% Admonition style "yellowicon" has colored background, yellow icons, and no farme
% Admon "block"
\colorlet{yellowicon_block_background}{yellow!5}
% \fboxsep sets the space between the text and the box
\newenvironment{blockshaded}
{\def\FrameCommand{\fboxsep=3mm\colorbox{yellowicon_block_background}}
 \MakeFramed {\advance\hsize-\width \FrameRestore}}{\endMakeFramed}

\newenvironment{block_yellowiconadmon}[1][Block]{
\begin{blockshaded}

 \textbf{#1}\par
\nobreak\ignorespaces
}
{
\end{blockshaded}
}

% --- end of definitions of admonition environments ---

% prevent orhpans and widows
\clubpenalty = 10000
\widowpenalty = 10000

% http://www.ctex.org/documents/packages/layout/titlesec.pdf
\usepackage{titlesec}  % needed for colored section headings
%\usepackage[]{titlesec}  % reduce the spacing around section headings

% --- section/subsection headings with blue color ---
\definecolor{seccolor}{cmyk}{.9,.5,0,.35}  % siamltexmm.sty section color
\titleformat{name=\section}
{\color{seccolor}\normalfont\Large\bfseries}
{\color{seccolor}\thesection}{1em}{}
\titleformat{name=\subsection}
{\color{seccolor}\normalfont\large\bfseries}
{\color{seccolor}\thesubsection}{1em}{}
\titleformat{name=\paragraph}[runin]
{\color{seccolor}\normalfont\normalsize\bfseries}
{}{}{\indent}

% let the header have a thick gray hrule with section and page in blue above
\renewcommand{\headrulewidth}{1pt}
\renewcommand{\headrule}{{\color{gray!50}%
\hrule width\headwidth height\headrulewidth \vskip-\headrulewidth}}
\fancyhead[LE,RO]{{\color{seccolor}\nouppercase{\rightmark}}} %section
\fancyhead[RE,LO]{{\color{seccolor}\thepage}}


% --- color every two table rows ---
\let\oldtabular\tabular
\let\endoldtabular\endtabular
\definecolor{appleblue}{rgb}{0.93,0.95,1.0}  % Apple blue
\renewenvironment{tabular}{\rowcolors{2}{white}{appleblue}%
\oldtabular}{\endoldtabular}


% --- end of standard preamble for documents ---


% insert custom LaTeX commands...

\raggedbottom
\makeindex
\usepackage[totoc]{idxlayout}   % for index in the toc
\usepackage[nottoc]{tocbibind}  % for references/bibliography in the toc



\begin{document}

% matching end for #ifdef PREAMBLE
% #endif

\newcommand{\exercisesection}[1]{\subsection*{#1}}







% ----------------- title -------------------------

\title{{\color{seccolor} Experiments with Schemes for Exponential Decay}}

% ----------------- author(s) -------------------------

\author{Hans Petter Langtangen\footnote{Email: \texttt{hpl@simula.no}. Center for Biomedical Computing, Simula Research Laboratory and Department of Informatics, University of Oslo.}}

% ----------------- end author(s) -------------------------

\date{Jul 2, 2016}
\maketitle

\begin{quote}\small
This report investigates the accuracy of three finite difference
schemes for the ordinary differential equation $u'=-au$ with the
aid of numerical experiments. Numerical artifacts are in particular
demonstrated.
\end{quote}

\tableofcontents


\vspace{1cm} % after toc





% !split
\section{Mathematical problem}
\label{math:problem}


\index{model problem} \index{exponential decay}

We address the initial-value problem

\begin{align}
u'(t) &= -au(t), \quad t \in (0,T], \label{ode}\\
u(0)  &= I,                         \label{initial:value}
\end{align}
where $a$, $I$, and $T$ are prescribed parameters, and $u(t)$ is
the unknown function to be estimated. This mathematical model
is relevant for physical phenomena featuring exponential decay
in time, e.g., vertical pressure variation in the atmosphere,
cooling of an object, and radioactive decay.

\section{Numerical solution method}
\label{numerical:problem}


\index{mesh in time} \index{$\theta$-rule} \index{numerical scheme}
\index{finite difference scheme}

We introduce a mesh in time with points $0 = t_0 < t_1 \cdots < t_{N_t}=T$.
For simplicity, we assume constant spacing $\Delta t$ between the
mesh points: $\Delta t = t_{n}-t_{n-1}$, $n=1,\ldots,N_t$. Let
$u^n$ be the numerical approximation to the exact solution at $t_n$.

The $\theta$-rule \cite{Iserles_2009}
is used to solve (\ref{ode}) numerically:

\[
u^{n+1} = \frac{1 - (1-\theta) a\Delta t}{1 + \theta a\Delta t}u^n,
\]
for $n=0,1,\ldots,N_t-1$. This scheme corresponds to

\begin{itemize}
  \item The \href{{http://en.wikipedia.org/wiki/Forward_Euler_method}}{Forward Euler}
    scheme when $\theta=0$

  \item The \href{{http://en.wikipedia.org/wiki/Backward_Euler_method}}{Backward Euler}
    scheme when $\theta=1$

  \item The \href{{http://en.wikipedia.org/wiki/Crank-Nicolson}}{Crank-Nicolson}
    scheme when $\theta=1/2$
\end{itemize}


\section{Implementation}




The numerical method is implemented in a Python function
\cite{Langtangen_2014} \texttt{solver} (found in the \href{{http://bit.ly/29ayDx3}}{\nolinkurl{model.py}} Python module file):


\begin{minted}[fontsize=\fontsize{9pt}{9pt},linenos=false,mathescape,baselinestretch=1.0,fontfamily=tt,xleftmargin=2mm]{python}
def solver(I, a, T, dt, theta):
    """Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt."""
    dt = float(dt)            # avoid integer division
    Nt = int(round(T/dt))     # no of time intervals
    T = Nt*dt                 # adjust T to fit time step dt
    u = zeros(Nt+1)           # array of u[n] values
    t = linspace(0, T, Nt+1)  # time mesh

    u[0] = I                  # assign initial condition
    for n in range(0, Nt):    # n=0,1,...,Nt-1
        u[n+1] = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)*u[n]
    return u, t
\end{minted}


\section{Numerical experiments}


\index{numerical experiments}

A set of numerical experiments has been carried out,
where $I$, $a$, and $T$ are fixed, while $\Delta t$ and
$\theta$ are varied. In particular, $I=1$, $a=2$,
$\Delta t = 1.25, 0.75, 0.5, 0.1$.
Figure~\ref{fig:BE} contains four plots, corresponding to
four decreasing $\Delta t$ values. The red dashed line
represent the numerical solution computed by the Backward
Euler scheme, while the blue line is the exact solution.
The corresponding results for the Crank-Nicolson and
Forward Euler methods appear in Figures~\ref{fig:CN}
and~\ref{fig:FE}.



\index{Backward Euler method}


\begin{figure}[!ht]  % fig:BE
  \centerline{\includegraphics[width=0.9\linewidth]{BE.pdf}}
  \caption{
  The Backward Euler method for decreasing time step values. \label{fig:BE}
  }
\end{figure}
%\clearpage % flush figures fig:BE




\index{Crank-Nicolson method}


\begin{figure}[!ht]  % fig:CN
  \centerline{\includegraphics[width=0.9\linewidth]{CN.pdf}}
  \caption{
  The Crank-Nicolson method for decreasing time step values. \label{fig:CN}
  }
\end{figure}
%\clearpage % flush figures fig:CN




\index{Forward Euler method}


\begin{figure}[!ht]  % fig:FE
  \centerline{\includegraphics[width=0.9\linewidth]{FE.pdf}}
  \caption{
  The Forward Euler method for decreasing time step values. \label{fig:FE}
  }
\end{figure}
%\clearpage % flush figures fig:FE




\section{Error vs $\Delta t$}

\index{error vs time step}

How the error

\[ E^n = \left(\int_0^T (Ie^{-at} - u^n)^2dt\right)^{\frac{1}{2}}\]
varies with $\Delta t$ for the three numerical methods
is shown in Figure~\ref{fig:error}.



\begin{warning_yellowiconadmon}[Observe:]
The data points for the three largest $\Delta t$ values in the
Forward Euler method are not relevant as the solution behaves
non-physically.
\end{warning_yellowiconadmon} % title: Observe:




\begin{figure}[!ht]  % fig:error
  \centerline{\includegraphics[width=0.9\linewidth]{error.pdf}}
  \caption{
  Variation of the error with the time step. \label{fig:error}
  }
\end{figure}
%\clearpage % flush figures fig:error



The $E$ numbers corresponding to Figure~\ref{fig:error}
are given in the table below.


\begin{quote}
\begin{tabular}{rrrr}
\hline
\multicolumn{1}{c}{ $\Delta t$ } & \multicolumn{1}{c}{ $\theta=0$ } & \multicolumn{1}{c}{ $\theta=0.5$ } & \multicolumn{1}{c}{ $\theta=1$ } \\
\hline
1.25       & 7.4630     & 0.2161       & 0.2440     \\
0.75       & 0.6632     & 0.0744       & 0.1875     \\
0.50       & 0.2797     & 0.0315       & 0.1397     \\
0.10       & 0.0377     & 0.0012       & 0.0335     \\
\hline
\end{tabular}
\end{quote}




\begin{summary_yellowiconadmon}[Summary]
\begin{enumerate}
\item $\theta =1$: $E\sim \Delta t$ (first-order convergence).

\item $\theta =0.5$: $E\sim \Delta t^2$ (second-order convergence).

\item $\theta =1$ is always stable and gives qualitatively corrects results.

\item $\theta =0.5$ never blows up, but may give oscillating solutions
   if $\Delta t$ is not sufficiently small.

\item $\theta =0$ suffers from fast-growing solution if $\Delta t$ is
   not small enough, but even below this limit one can have oscillating
   solutions (unless $\Delta t$ is sufficiently small).
\end{enumerate}


\end{summary_yellowiconadmon} % title: Summary





\bibliographystyle{plain}
\bibliography{.publish_references}




% #ifdef PREAMBLE
\cleardoublepage\phantomsection  % trick to get correct link to Index
\printindex

\end{document}
% #endif