#!/usr/bin/env python
# coding: utf-8

# # Analysis of schemes for the diffusion equation
# <div id="diffu:pde1:analysis"></div>
# 
# 
# The numerical experiments in the sections [diffu:pde1:FE:experiments](#diffu:pde1:FE:experiments) and [diffu:pde1:theta:experiments](#diffu:pde1:theta:experiments)
# reveal that there are some
# numerical problems with the Forward Euler and Crank-Nicolson schemes:
# sawtooth-like noise is sometimes present in solutions that are,
# from a mathematical point of view, expected to be smooth.
# This section presents a mathematical analysis that explains the
# observed behavior and arrives at criteria for obtaining numerical
# solutions that reproduce the qualitative properties of the exact
# solutions. In short, we shall explain what is observed in
# Figures [diffu:pde1:FE:fig:F=0.5](#diffu:pde1:FE:fig:F=0.5)-[diffu:pde1:CN:fig:F=10](#diffu:pde1:CN:fig:F=10).
# <!-- [diffu:pde1:FE:fig:F=0.5](#diffu:pde1:FE:fig:F=0.5), -->
# <!-- [diffu:pde1:FE:fig:F=0.25](#diffu:pde1:FE:fig:F=0.25), -->
# <!-- [diffu:pde1:FE:fig:F=0.51](#diffu:pde1:FE:fig:F=0.51), -->
# <!-- [diffu:pde1:FE:fig:gauss:F=0.5](#diffu:pde1:FE:fig:gauss:F=0.5), -->
# <!-- [diffu:pde1:BE:fig:F=0.5](#diffu:pde1:BE:fig:F=0.5), -->
# <!-- [diffu:pde1:CN:fig:F=3](#diffu:pde1:CN:fig:F=3), -->
# <!-- and -->
# <!-- [diffu:pde1:CN:fig:F=10](#diffu:pde1:CN:fig:F=10). -->
# 
# ## Properties of the solution
# <div id="diffu:pde1:analysis:uex"></div>
# 
# 
# A particular characteristic of diffusive processes, governed
# by an equation like

# <!-- Equation labels as ordinary links -->
# <div id="diffu:pde1:eq"></div>
# 
# $$
# \begin{equation}
# u_t = \dfc u_{xx},
# \label{diffu:pde1:eq} \tag{1}
# \end{equation}
# $$

# is that the initial shape $u(x,0)=I(x)$ spreads out in space with
# time, along with a decaying amplitude.  Three different examples will
# illustrate the spreading of $u$ in space and the decay in time.
# 
# ### Similarity solution
# 
# The diffusion equation ([1](#diffu:pde1:eq)) admits solutions
# that depend on $\eta = (x-c)/\sqrt{4\dfc t}$ for a given value
# of $c$. One particular solution
# is

# <!-- Equation labels as ordinary links -->
# <div id="diffu:pdf1:erf:sol"></div>
# 
# $$
# \begin{equation}
# u(x,t) = a\,\mbox{erf}(\eta) + b,
# \label{diffu:pdf1:erf:sol} \tag{2}
# \end{equation}
# $$

# where

# <!-- Equation labels as ordinary links -->
# <div id="diffu:analysis:erf:def"></div>
# 
# $$
# \begin{equation}
# \mbox{erf}(\eta) = \frac{2}{\sqrt{\pi}}\int_0^\eta e^{-\zeta^2}d\zeta,
# \label{diffu:analysis:erf:def} \tag{3}
# \end{equation}
# $$

# is the *error function*, and $a$ and $b$ are arbitrary constants.
# The error function lies in $(-1,1)$, is odd around $\eta =0$, and
# goes relatively quickly to $\pm 1$:

# $$
# \begin{align*}
# \lim_{\eta\rightarrow -\infty}\mbox{erf}(\eta) &=-1,\\ 
# \lim_{\eta\rightarrow \infty}\mbox{erf}(\eta) &=1,\\ 
# \mbox{erf}(\eta) &= -\mbox{erf}(-\eta),\\ 
# \mbox{erf}(0) &=0,\\ 
# \mbox{erf}(2) &=0.99532227,\\ 
# \mbox{erf}(3) &=0.99997791
# \thinspace .
# \end{align*}
# $$

# As $t\rightarrow 0$, the error function approaches a step function centered
# at $x=c$. For a diffusion problem posed on the unit interval $[0,1]$,
# we may choose the step at $x=1/2$ (meaning $c=1/2$), $a=-1/2$, $b=1/2$.
# Then

# <!-- Equation labels as ordinary links -->
# <div id="diffu:analysis:pde1:step:erf:sol"></div>
# 
# $$
# \begin{equation}
# u(x,t) = \frac{1}{2}\left(1 -
# \mbox{erf}\left(\frac{x-\frac{1}{2}}{\sqrt{4\dfc t}}\right)\right) =
# \frac{1}{2}\mbox{erfc}\left(\frac{x-\frac{1}{2}}{\sqrt{4\dfc t}}\right),
# \label{diffu:analysis:pde1:step:erf:sol} \tag{4}
# \end{equation}
# $$

# where we have introduced the *complementary error function*
# $\mbox{erfc}(\eta) = 1-\mbox{erf}(\eta)$.
# The solution ([4](#diffu:analysis:pde1:step:erf:sol))
# implies the boundary conditions

# <!-- Equation labels as ordinary links -->
# <div id="diffu:analysis:pde1:p1:erf:uL"></div>
# 
# $$
# \begin{equation}
# u(0,t) = \frac{1}{2}\left(1 - \mbox{erf}\left(\frac{-1/2}{\sqrt{4\dfc t}}\right)\right),
# \label{diffu:analysis:pde1:p1:erf:uL} \tag{5} 
# \end{equation}
# $$

# <!-- Equation labels as ordinary links -->
# <div id="diffu:analysis:pde1:p1:erf:uR"></div>
# 
# $$
# \begin{equation}  
# u(1,t) = \frac{1}{2}\left(1 - \mbox{erf}\left(\frac{1/2}{\sqrt{4\dfc t}}\right)\right)
# \label{diffu:analysis:pde1:p1:erf:uR} \tag{6}
# \thinspace .
# \end{equation}
# $$

# For small enough $t$, $u(0,t)\approx 1$ and $u(1,t)\approx 0$, but as
# $t\rightarrow\infty$, $u(x,t)\rightarrow 1/2$ on $[0,1]$.
# 
# ### Solution for a Gaussian pulse
# 
# The standard diffusion equation $u_t = \dfc u_{xx}$ admits a
# Gaussian function as solution:

# <!-- Equation labels as ordinary links -->
# <div id="diffu:pde1:sol:Gaussian"></div>
# 
# $$
# \begin{equation}
# u(x,t) = \frac{1}{\sqrt{4\pi\dfc t}} \exp{\left({-\frac{(x-c)^2}{4\dfc t}}\right)}
# \label{diffu:pde1:sol:Gaussian} \tag{7}
# \thinspace .
# \end{equation}
# $$

# At $t=0$ this is a Dirac delta function, so for computational
# purposes one must start to view the solution at some time $t=t_\epsilon>0$.
# Replacing $t$ by $t_\epsilon +t$ in ([7](#diffu:pde1:sol:Gaussian))
# makes it easy to operate with a (new) $t$ that starts at $t=0$
# with an initial condition with a finite width.
# The important feature of ([7](#diffu:pde1:sol:Gaussian)) is that
# the standard deviation $\sigma$ of a sharp initial Gaussian pulse
# increases in time according to $\sigma = \sqrt{2\dfc t}$, making
# the pulse diffuse and flatten out.
# 
# <!-- Mention combinations of such kernels to build up a general analytical sol? -->
# <!-- Or maybe an exercise for verification. -->
# 
# ### Solution for a sine component
# 
# Also, ([1](#diffu:pde1:eq)) admits a solution of the form

# <!-- Equation labels as ordinary links -->
# <div id="diffu:pde1:sol1"></div>
# 
# $$
# \begin{equation}
# u(x,t) = Qe^{-at}\sin\left( kx\right)
# \label{diffu:pde1:sol1} \tag{8}
# \thinspace .
# \end{equation}
# $$

# The parameters $Q$ and $k$ can be freely chosen, while
# inserting ([8](#diffu:pde1:sol1)) in ([1](#diffu:pde1:eq)) gives the constraint

# $$
# a = -\dfc k^2
# \thinspace .
# $$

# A very important feature is that the initial shape $I(x)=Q\sin\left( kx\right)$
# undergoes a damping $\exp{(-\dfc k^2t)}$, meaning that
# rapid oscillations in space, corresponding to large $k$, are very much
# faster dampened than slow oscillations in space, corresponding to small
# $k$. This feature leads to a smoothing of the initial condition with time.
# (In fact, one can use a few steps of the diffusion equation as
# a method for removing noise in signal processing.)
# To judge how good a numerical method is, we may look at its ability to
# smoothen or dampen the solution in the same way as the PDE does.
# 
# The following example illustrates the damping properties of
# ([8](#diffu:pde1:sol1)). We consider the specific problem

# $$
# \begin{align*}
# u_t &= u_{xx},\quad x\in (0,1),\ t\in (0,T],\\ 
# u(0,t) &= u(1,t) = 0,\quad t\in (0,T],\\ 
# u(x,0) & = \sin (\pi x) + 0.1\sin(100\pi x)
# \thinspace .
# \end{align*}
# $$

# The initial condition has been chosen such that adding
# two solutions like ([8](#diffu:pde1:sol1)) constructs
# an analytical solution to the problem:

# <!-- Equation labels as ordinary links -->
# <div id="diffu:pde1:sol2"></div>
# 
# $$
# \begin{equation}
# u(x,t) = e^{-\pi^2 t}\sin (\pi x) + 0.1e^{-\pi^2 10^4 t}\sin (100\pi x)
# \label{diffu:pde1:sol2} \tag{9}
# \thinspace .
# \end{equation}
# $$

# [Figure](#diffu:pde1:fig:damping) illustrates the rapid damping of
# rapid oscillations $\sin (100\pi x)$ and the very much slower damping of the
# slowly varying $\sin (\pi x)$ term. After about $t=0.5\cdot10^{-4}$ the rapid
# oscillations do not have a visible amplitude, while we have to wait
# until $t\sim 0.5$ before the amplitude of the long wave $\sin (\pi x)$
# becomes very small.
# 
# <!-- dom:FIGURE: [fig-diffu/diffusion_damping.png, width=800] Evolution of the solution of a diffusion problem: initial condition (upper left), 1/100 reduction of the small waves (upper right), 1/10 reduction of the long wave (lower left), and 1/100 reduction of the long wave (lower right). <div id="diffu:pde1:fig:damping"></div> -->
# <!-- begin figure -->
# <div id="diffu:pde1:fig:damping"></div>
# 
# <p>Evolution of the solution of a diffusion problem: initial condition (upper left), 1/100 reduction of the small waves (upper right), 1/10 reduction of the long wave (lower left), and 1/100 reduction of the long wave (lower right).</p>
# <img src="fig-diffu/diffusion_damping.png" width=800>
# 
# <!-- end figure -->
# 
# 
# <!-- x/sqrt(t) solution, kernel with integral -->
# 
# ## Analysis of discrete equations
# 
# A counterpart to ([8](#diffu:pde1:sol1)) is the complex representation
# of the same function:

# $$
# u(x,t) = Qe^{-at}e^{ikx},
# $$

# where $i=\sqrt{-1}$ is the imaginary unit.
# We can add such functions, often referred to as wave components,
# to make a Fourier representation
# of a general solution of the diffusion equation:

# <!-- Equation labels as ordinary links -->
# <div id="diffu:pde1:u:Fourier"></div>
# 
# $$
# \begin{equation}
# u(x,t) \approx \sum_{k\in K} b_k e^{-\dfc k^2t}e^{ikx},
# \label{diffu:pde1:u:Fourier} \tag{10}
# \end{equation}
# $$

# where $K$ is a set of an infinite number of $k$ values needed to construct
# the solution. In practice, however, the series is truncated and
# $K$ is a finite set of $k$ values
# needed to build a good approximate solution.
# Note that ([9](#diffu:pde1:sol2)) is a special case of
# ([10](#diffu:pde1:u:Fourier)) where $K=\{\pi, 100\pi\}$, $b_{\pi}=1$,
# and $b_{100\pi}=0.1$.
# 
# The amplitudes $b_k$ of the individual Fourier waves must be determined
# from the initial condition. At $t=0$ we have $u\approx\sum_kb_k\exp{(ikx)}$
# and find $K$ and $b_k$ such that

# <!-- Equation labels as ordinary links -->
# <div id="_auto1"></div>
# 
# $$
# \begin{equation}
# I(x) \approx \sum_{k\in K} b_k e^{ikx}\thinspace .
# \label{_auto1} \tag{11}
# \end{equation}
# $$

# (The relevant formulas for $b_k$ come from Fourier analysis, or
# equivalently, a least-squares method for approximating $I(x)$
# in a function space with basis $\exp{(ikx)}$.)
# 
# Much insight about the behavior of numerical methods can be obtained
# by investigating how a wave component $\exp{(-\dfc k^2
# t)}\exp{(ikx)}$ is treated by the numerical scheme. mathcal{I}_t appears that
# such wave components are also solutions of the schemes, but the
# damping factor $\exp{(-\dfc k^2 t)}$ varies among the schemes.  To
# ease the forthcoming algebra, we write the damping factor as
# $A^n$. The exact amplification factor corresponding to $A$ is $\Aex =
# \exp{(-\dfc k^2\Delta t)}$.
# 
# 
# ## Analysis of the finite difference schemes
# <div id="diffu:pde1:analysis:details"></div>
# 
# We have seen that a general solution of the diffusion equation
# can be built as a linear combination of basic components

# $$
# e^{-\dfc k^2t}e^{ikx} \thinspace .
# $$

# A fundamental question is whether such components are also solutions of
# the finite difference schemes. This is indeed the case, but the
# amplitude $\exp{(-\dfc k^2t)}$ might be modified (which also happens when
# solving the ODE counterpart $u'=-\dfc u$).
# We therefore look for numerical solutions of the form

# <!-- Equation labels as ordinary links -->
# <div id="diffu:pde1:analysis:uni"></div>
# 
# $$
# \begin{equation}
# u^n_q = A^n e^{ikq\Delta x} = A^ne^{ikx},
# \label{diffu:pde1:analysis:uni} \tag{12}
# \end{equation}
# $$

# where the amplification factor $A$
# must be determined by inserting the component into an actual scheme.
# Note that $A^n$ means $A$ raised to the power of $n$, $n$ being the
# index in the time mesh, while the superscript $n$ in $u^n_q$ just
# denotes $u$ at time $t_n$.
# 
# ### Stability
# 
# The exact amplification factor is $\Aex=\exp{(-\dfc^2 k^2\Delta t)}$.
# We should therefore require $|A| < 1$ to have a decaying numerical
# solution as well. If
# $-1\leq A<0$, $A^n$ will change sign from time level to
# time level, and we get stable, non-physical oscillations in the numerical
# solutions that are not present in the exact solution.
# 
# 
# ### Accuracy
# 
# To determine how accurately a finite difference scheme treats one
# wave component ([12](#diffu:pde1:analysis:uni)), we see that the basic
# deviation from the exact solution is reflected in how well
# $A^n$ approximates $\Aex^n$,
# or how well $A$ approximates $\Aex$.
# We can plot $\Aex$ and the various expressions for $A$, and we can
# make Taylor expansions of $A/\Aex$ to see the error more analytically.
# 
# <!-- We shall in particular investigate the error $\Aex - A$ in the -->
# <!-- amplification factor. -->
# 
# ### Truncation error
# 
# As an alternative to examining the accuracy of the damping of a wave
# component, we can perform a general truncation error analysis as
# explained in "Truncation error analysis": ""
# [[Langtangen_deqbook_trunc]](#Langtangen_deqbook_trunc). Such results are more general, but
# less detailed than what we get from the wave component analysis.  The
# truncation error can almost always be computed and represents the
# error in the numerical model when the exact solution is substituted
# into the equations. In particular, the truncation error analysis tells
# the order of the scheme, which is of fundamental importance when
# verifying codes based on empirical estimation of convergence rates.
# 
# ## Analysis of the Forward Euler scheme
# <div id="diffu:pde1:analysis:FE"></div>
# 
# 
# <!-- 2DO: refer to vib and wave -->
# 
# 
# The Forward Euler finite difference scheme for $u_t = \dfc u_{xx}$ can
# be written as

# $$
# [D_t^+ u = \dfc D_xD_x u]^n_q\thinspace .
# $$

# Inserting a wave component ([12](#diffu:pde1:analysis:uni))
# in the scheme demands calculating the terms

# $$
# e^{ikq\Delta x}[D_t^+ A]^n = e^{ikq\Delta x}A^n\frac{A-1}{\Delta t},
# $$

# and

# $$
# A^nD_xD_x [e^{ikx}]_q = A^n\left( - e^{ikq\Delta x}\frac{4}{\Delta x^2}
# \sin^2\left(\frac{k\Delta x}{2}\right)\right)
# \thinspace .
# $$

# Inserting these terms in the discrete equation and
# dividing by $A^n e^{ikq\Delta x}$ leads to

# $$
# \frac{A-1}{\Delta t} = -\dfc \frac{4}{\Delta x^2}\sin^2\left(
# \frac{k\Delta x}{2}\right),
# $$

# and consequently

# <!-- Equation labels as ordinary links -->
# <div id="_auto2"></div>
# 
# $$
# \begin{equation}
# A = 1 -4F\sin^2 p
# \label{_auto2} \tag{13}
# \end{equation}
# $$

# where

# <!-- Equation labels as ordinary links -->
# <div id="_auto3"></div>
# 
# $$
# \begin{equation}
# F = \frac{\dfc\Delta t}{\Delta x^2}
# \label{_auto3} \tag{14}
# \end{equation}
# $$

# is the *numerical Fourier number*, and $p=k\Delta x/2$.
# The complete numerical solution is then

# <!-- Equation labels as ordinary links -->
# <div id="_auto4"></div>
# 
# $$
# \begin{equation}
# u^n_q = \left(1 -4F\sin^2 p\right)^ne^{ikq\Delta x}
# \thinspace .
# \label{_auto4} \tag{15}
# \end{equation}
# $$

# ### Stability
# 
# We easily see that $A\leq 1$. However, the $A$ can be less than $-1$,
# which will lead
# to growth of a numerical wave component. The criterion $A\geq -1$ implies

# $$
# 4F\sin^2 (p/2)\leq 2
# \thinspace .
# $$

# The worst case is when $\sin^2 (p/2)=1$, so a sufficient criterion for
# stability is

# <!-- Equation labels as ordinary links -->
# <div id="_auto5"></div>
# 
# $$
# \begin{equation}
# F\leq {\frac{1}{2}},
# \label{_auto5} \tag{16}
# \end{equation}
# $$

# or expressed as a condition on $\Delta t$:

# <!-- Equation labels as ordinary links -->
# <div id="_auto6"></div>
# 
# $$
# \begin{equation}
# \Delta t\leq \frac{\Delta x^2}{2\dfc}\thinspace .
# \label{_auto6} \tag{17}
# \end{equation}
# $$

# Note that halving the spatial mesh size, $\Delta x \rightarrow {\frac{1}{2}}
# \Delta x$, requires $\Delta t$ to be reduced by a factor of $1/4$.
# The method hence becomes very expensive for fine spatial meshes.
# 
# <!-- 2DO: verification based on exact solutions -->
# 
# ### Accuracy
# 
# Since $A$ is expressed in terms of $F$ and the parameter we now call
# $p=k\Delta x/2$, we should also express $\Aex$ by $F$ and $p$. The exponent
# in $\Aex$ is $-\dfc k^2\Delta t$, which equals $-F k^2\Delta x^2=-F4p^2$.
# Consequently,

# $$
# \Aex = \exp{(-\dfc k^2\Delta t)} = \exp{(-4Fp^2)}
# \thinspace .
# $$

# All our $A$ expressions as well as $\Aex$ are now functions of the two
# dimensionless parameters $F$ and $p$.
# 
# Computing
# the Taylor series expansion of $A/\Aex$ in terms of $F$
# can easily be done with aid of `sympy`:

# In[1]:


def A_exact(F, p):
    return exp(-4*F*p**2)

def A_FE(F, p):
    return 1 - 4*F*sin(p)**2

from sympy import *
F, p = symbols('F p')
A_err_FE = A_FE(F, p)/A_exact(F, p)
print(A_err_FE.series(F, 0, 6))


# The result is

# $$
# \frac{A}{\Aex} = 1 - 4 F \sin^{2}p + 2F p^{2} - 16F^{2} p^{2} \sin^{2}p + 8 F^{2} p^{4} + \cdots
# $$

# Recalling that $F=\dfc\Delta t/\Delta x^2$, $p=k\Delta x/2$, and that
# $\sin^2p\leq 1$, we
# realize that the dominating terms in $A/\Aex$ are at most

# $$
# 1 - 4\dfc \frac{\Delta t}{\Delta x^2} +
# \dfc\Delta t - 4\dfc^2\Delta t^2
# + \dfc^2 \Delta t^2\Delta x^2 + \cdots
# \thinspace .
# $$

# ### Truncation error
# 
# We follow the theory explained in
# "Truncation error analysis": ""
# [[Langtangen_deqbook_trunc]](#Langtangen_deqbook_trunc). The recipe is to set up the
# scheme in operator notation and use formulas from
# "Overview of leading-order error terms in finite difference formulas": ""
# [[Langtangen_deqbook_trunc]](#Langtangen_deqbook_trunc) to derive an expression for
# the residual. The details are documented in
# "Linear diffusion equation in 1D": ""
# [[Langtangen_deqbook_trunc]](#Langtangen_deqbook_trunc). We end up with a truncation error

# $$
# R^n_i = \Oof{\Delta t} + \Oof{\Delta x^2}\thinspace .
# $$

# Although this is not the true error $\uex(x_i,t_n) - u^n_i$, it indicates
# that the true error is of the form

# $$
# E = C_t\Delta t + C_x\Delta x^2
# $$

# for two unknown constants $C_t$ and $C_x$.
# 
# 
# ## Analysis of the Backward Euler scheme
# <div id="diffu:pde1:analysis:BE"></div>
# 
# Discretizing $u_t = \dfc u_{xx}$ by a Backward Euler scheme,

# $$
# [D_t^- u = \dfc D_xD_x u]^n_q,
# $$

# and inserting a wave component ([12](#diffu:pde1:analysis:uni)),
# leads to calculations similar to those arising from the Forward Euler scheme,
# but since

# $$
# e^{ikq\Delta x}[D_t^- A]^n = A^ne^{ikq\Delta x}\frac{1 - A^{-1}}{\Delta t},
# $$

# we get

# $$
# \frac{1-A^{-1}}{\Delta t} = -\dfc \frac{4}{\Delta x^2}\sin^2\left(
# \frac{k\Delta x}{2}\right),
# $$

# and then

# <!-- Equation labels as ordinary links -->
# <div id="diffu:pde1:analysis:BE:A"></div>
# 
# $$
# \begin{equation}
# A = \left(1  + 4F\sin^2p\right)^{-1}
# \label{diffu:pde1:analysis:BE:A} \tag{18}
# \thinspace .
# \end{equation}
# $$

# The complete numerical solution can be written

# <!-- Equation labels as ordinary links -->
# <div id="_auto7"></div>
# 
# $$
# \begin{equation}
# u^n_q = \left(1  + 4F\sin^2 p\right)^{-n}
# e^{ikq\Delta x} \thinspace .
# \label{_auto7} \tag{19}
# \end{equation}
# $$

# ### Stability
# 
# We see from ([18](#diffu:pde1:analysis:BE:A)) that $0<A<1$, which means
# that all numerical wave components are stable and non-oscillatory
# for any $\Delta t >0$.
# 
# ### Truncation error
# 
# The derivation of the truncation error for the Backward Euler scheme is almost
# identical to that for the Forward Euler scheme. We end up with

# $$
# R^n_i = \Oof{\Delta t} + \Oof{\Delta x^2}\thinspace .
# $$

# ## Analysis of the Crank-Nicolson scheme
# <div id="diffu:pde1:analysis:CN"></div>
# 
# The Crank-Nicolson scheme can be written as

# $$
# [D_t u = \dfc D_xD_x \overline{u}^x]^{n+\frac{1}{2}}_q,
# $$

# or

# $$
# [D_t u]^{n+\frac{1}{2}}_q = \frac{1}{2}\dfc\left( [D_xD_x u]^{n}_q +
# [D_xD_x u]^{n+1}_q\right)
# \thinspace .
# $$

# Inserting ([12](#diffu:pde1:analysis:uni)) in the time derivative approximation
# leads to

# $$
# [D_t A^n e^{ikq\Delta x}]^{n+\frac{1}{2}} = A^{n+\frac{1}{2}} e^{ikq\Delta x}\frac{A^{\frac{1}{2}}-A^{-\frac{1}{2}}}{\Delta t} = A^ne^{ikq\Delta x}\frac{A-1}{\Delta t}
# \thinspace .
# $$

# Inserting ([12](#diffu:pde1:analysis:uni)) in the other terms
# and dividing by
# $A^ne^{ikq\Delta x}$ gives the relation

# $$
# \frac{A-1}{\Delta t} = -\frac{1}{2}\dfc\frac{4}{\Delta x^2}
# \sin^2\left(\frac{k\Delta x}{2}\right)
# (1 + A),
# $$

# and after some more algebra,

# <!-- Equation labels as ordinary links -->
# <div id="_auto8"></div>
# 
# $$
# \begin{equation}
# A = \frac{ 1 - 2F\sin^2p}{1 + 2F\sin^2p}
# \thinspace .
# \label{_auto8} \tag{20}
# \end{equation}
# $$

# The exact numerical solution is hence

# <!-- Equation labels as ordinary links -->
# <div id="_auto9"></div>
# 
# $$
# \begin{equation}
# u^n_q = \left(\frac{ 1 - 2F\sin^2p}{1 + 2F\sin^2p}\right)^ne^{ikq\Delta x}
# \thinspace .
# \label{_auto9} \tag{21}
# \end{equation}
# $$

# ### Stability
# 
# The criteria $A>-1$ and $A<1$ are fulfilled for any $\Delta t >0$.
# Therefore, the solution cannot grow, but it will oscillate if
# $1-2F\sin^p < 0$. To avoid such non-physical oscillations, we must demand
# $F\leq\frac{1}{2}$.
# 
# ### Truncation error
# 
# The truncation error is derived in
# "Linear diffusion equation in 1D": ""
# [[Langtangen_deqbook_trunc]](#Langtangen_deqbook_trunc):

# $$
# R^{n+\frac{1}{2}}_i = \Oof{\Delta x^2} + \Oof{\Delta t^2}\thinspace .
# $$

# ## Analysis of the Leapfrog scheme
# <div id="diffu:pde1:analysis:leapfrog"></div>
# 
# 
# An attractive feature of the Forward Euler scheme is the explicit
# time stepping and no need for solving linear systems. However, the
# accuracy in time is only $\Oof{\Delta t}$. We can get an explicit
# *second-order* scheme in time by using the Leapfrog method:

# $$
# [D_{2t} u = \dfc D_xDx u + f]^n_q\thinspace .
# $$

# Written out,

# $$
# u_q^{n+1} = u_q^{n-1} + \frac{2\dfc\Delta t}{\Delta x^2}
# (u^{n}_{q+1} - 2u^n_q + u^n_{q-1}) + f(x_q,t_n)\thinspace .
# $$

# We need some formula for the first step, $u^1_q$, but for that we can use
# a Forward Euler step.
# 
# Unfortunately, the Leapfrog scheme is always unstable for the
# diffusion equation. To see this, we insert a wave component $A^ne^{ikx}$
# and get

# $$
# \frac{A - A^{-1}}{\Delta t} = -\dfc \frac{4}{\Delta x^2}\sin^2 p,
# $$

# or

# $$
# A^2 + 4F \sin^2 p\, A - 1 = 0,
# $$

# which has roots

# $$
# A = -2F\sin^2 p \pm \sqrt{4F^2\sin^4 p + 1}\thinspace .
# $$

# Both roots have $|A|>1$ so the amplitude always grows, which is not in
# accordance with the physics of the problem.
# However, for a PDE with a first-order derivative in space, instead of
# a second-order one, the Leapfrog scheme performs very well.
# 
# ## Summary of accuracy of amplification factors
# 
# We can plot the various amplification factors against $p=k\Delta x/2$
# for different choices of the $F$ parameter. Figures
# [diffu:pde1:fig:A:err:C20](#diffu:pde1:fig:A:err:C20), [diffu:pde1:fig:A:err:C0.5](#diffu:pde1:fig:A:err:C0.5), and
# [diffu:pde1:fig:A:err:C0.1](#diffu:pde1:fig:A:err:C0.1) show how long and small waves are
# damped by the various schemes compared to the exact damping. As long
# as all schemes are stable, the amplification factor is positive,
# except for Crank-Nicolson when $F>0.5$.
# 
# <!-- dom:FIGURE: [fig-diffu/diffusion_A_F20_F2.png, width=800] Amplification factors for large time steps. <div id="diffu:pde1:fig:A:err:C20"></div> -->
# <!-- begin figure -->
# <div id="diffu:pde1:fig:A:err:C20"></div>
# 
# <p>Amplification factors for large time steps.</p>
# <img src="fig-diffu/diffusion_A_F20_F2.png" width=800>
# 
# <!-- end figure -->
# 
# 
# <!-- dom:FIGURE: [fig-diffu/diffusion_A_F05_F025.png, width=800] Amplification factors for time steps around the Forward Euler stability limit. <div id="diffu:pde1:fig:A:err:C0.5"></div> -->
# <!-- begin figure -->
# <div id="diffu:pde1:fig:A:err:C0.5"></div>
# 
# <p>Amplification factors for time steps around the Forward Euler stability limit.</p>
# <img src="fig-diffu/diffusion_A_F05_F025.png" width=800>
# 
# <!-- end figure -->
# 
# 
# <!-- dom:FIGURE: [fig-diffu/diffusion_A_F01_F001.png, width=800] Amplification factors for small time steps. <div id="diffu:pde1:fig:A:err:C0.1"></div> -->
# <!-- begin figure -->
# <div id="diffu:pde1:fig:A:err:C0.1"></div>
# 
# <p>Amplification factors for small time steps.</p>
# <img src="fig-diffu/diffusion_A_F01_F001.png" width=800>
# 
# <!-- end figure -->
# 
# 
# The effect of negative amplification factors is that $A^n$ changes
# sign from one time level to the next, thereby giving rise to
# oscillations in time in an animation of the solution. We see from
# [Figure](#diffu:pde1:fig:A:err:C20) that for $F=20$, waves with
# $p\geq \pi/4$ undergo a damping close to $-1$, which means that the
# amplitude does not decay and that the wave component jumps up and down
# (flips amplitude) in time. For $F=2$ we have a damping of a factor of
# 0.5 from one time level to the next, which is very much smaller than
# the exact damping. Short waves will therefore fail to be effectively
# dampened.  These waves will manifest themselves as high frequency
# oscillatory noise in the solution.
# 
# A value $p=\pi/4$ corresponds to four mesh points per wave length of
# $e^{ikx}$, while $p=\pi/2$ implies only two points per wave length,
# which is the smallest number of points we can have to represent the
# wave on the mesh.
# 
# To demonstrate the oscillatory behavior of the Crank-Nicolson scheme,
# we choose an initial condition that leads to short waves with
# significant amplitude. A discontinuous $I(x)$ will in particular serve
# this purpose: Figures [diffu:pde1:CN:fig:F=3](#diffu:pde1:CN:fig:F=3) and
# [diffu:pde1:CN:fig:F=10](#diffu:pde1:CN:fig:F=10) correspond to $F=3$ and $F=10$,
# respectively, and we see how short waves pollute the overall solution.
# 
# ## Analysis of the 2D diffusion equation
# <div id="diffu:2D:analysis"></div>
# 
# Diffusion in several dimensions is treated later, but it is appropriate to
# include the analysis here. We first consider the 2D diffusion equation

# $$
# u_{t} = \dfc(u_{xx} + u_{yy}),
# $$

# which has Fourier component solutions of the form

# $$
# u(x,y,t) = Ae^{-\dfc k^2t}e^{i(k_x x + k_yy)},
# $$

# and the schemes have discrete versions of this Fourier component:

# $$
# u^{n}_{q,r} = A\xi^{n}e^{i(k_x q\Delta x + k_y r\Delta y)}\thinspace .
# $$

# ### The Forward Euler scheme
# 
# For the Forward Euler discretization,

# $$
# [D_t^+u = \dfc(D_xD_x u + D_yD_y u)]_{q,r}^n,
# $$

# we get

# $$
# \frac{\xi - 1}{\Delta t}
# =
# -\dfc\frac{4}{\Delta x^2}\sin^2\left(\frac{k_x\Delta x}{2}\right) -
# \dfc\frac{4}{\Delta y^2}\sin^2\left(\frac{k_y\Delta y}{2}\right)\thinspace .
# $$

# Introducing

# $$
# p_x = \frac{k_x\Delta x}{2},\quad p_y = \frac{k_y\Delta y}{2},
# $$

# we can write the equation for $\xi$ more compactly as

# $$
# \frac{\xi - 1}{\Delta t}
# =
# -\dfc\frac{4}{\Delta x^2}\sin^2 p_x -
# \dfc\frac{4}{\Delta y^2}\sin^2 p_y,
# $$

# and solve for $\xi$:

# <!-- Equation labels as ordinary links -->
# <div id="diffu:2D:analysis:xi"></div>
# 
# $$
# \begin{equation}
# \xi = 1 - 4F_x\sin^2 p_x - 4F_y\sin^2 p_y\thinspace .
# \label{diffu:2D:analysis:xi} \tag{22}
# \end{equation}
# $$

# The complete numerical solution for a wave component is

# <!-- Equation labels as ordinary links -->
# <div id="diffu:2D:analysis:FE:numexact"></div>
# 
# $$
# \begin{equation}
# u^{n}_{q,r} = A(1 - 4F_x\sin^2 p_x - 4F_y\sin^2 p_y)^n
# e^{i(k_xq\Delta x + k_yr\Delta y)}\thinspace .
# \label{diffu:2D:analysis:FE:numexact} \tag{23}
# \end{equation}
# $$

# For stability we demand $-1\leq\xi\leq 1$, and $-1\leq\xi$ is the
# critical limit, since clearly $\xi \leq 1$, and the worst case
# happens when the sines are at their maximum. The stability criterion
# becomes

# <!-- Equation labels as ordinary links -->
# <div id="diffu:2D:analysis:FE:stab"></div>
# 
# $$
# \begin{equation}
# F_x + F_y \leq \frac{1}{2}\thinspace .
# \label{diffu:2D:analysis:FE:stab} \tag{24}
# \end{equation}
# $$

# For the special, yet common, case $\Delta x=\Delta y=h$, the
# stability criterion can be written as

# $$
# \Delta t \leq \frac{h^2}{2d\dfc},
# $$

# where $d$ is the number of space dimensions: $d=1,2,3$.
# 
# ### The Backward Euler scheme
# 
# The Backward Euler method,

# $$
# [D_t^-u = \dfc(D_xD_x u + D_yD_y u)]_{q,r}^n,
# $$

# results in

# $$
# 1 - \xi^{-1} = - 4F_x \sin^2 p_x - 4F_y \sin^2 p_y,
# $$

# and

# $$
# \xi = (1 + 4F_x \sin^2 p_x + 4F_y \sin^2 p_y)^{-1},
# $$

# which is always in $(0,1]$. The solution for a wave component becomes

# <!-- Equation labels as ordinary links -->
# <div id="diffu:2D:analysis:BN:numexact"></div>
# 
# $$
# \begin{equation}
# u^{n}_{q,r} = A(1 + 4F_x\sin^2 p_x + 4F_y\sin^2 p_y)^{-n}
# e^{i(k_xq\Delta x + k_yr\Delta y)}\thinspace .
# \label{diffu:2D:analysis:BN:numexact} \tag{25}
# \end{equation}
# $$

# ### The Crank-Nicolson scheme
# 
# With a Crank-Nicolson discretization,

# $$
# [D_tu]^{n+\frac{1}{2}}_{q,r} =
# \frac{1}{2} [\dfc(D_xD_x u + D_yD_y u)]_{q,r}^{n+1} +
# \frac{1}{2} [\dfc(D_xD_x u + D_yD_y u)]_{q,r}^n,
# $$

# we have, after some algebra,

# $$
# \xi = \frac{1 - 2(F_x\sin^2 p_x + F_x\sin^2p_y)}{1 + 2(F_x\sin^2 p_x + F_x\sin^2p_y)}\thinspace .
# $$

# The fraction on the right-hand side is always less than 1, so stability
# in the sense of non-growing wave components is guaranteed for all
# physical and numerical parameters. However,
# the fraction can become negative and result in non-physical
# oscillations. This phenomenon happens when

# $$
# F_x\sin^2 p_x + F_x\sin^2p_y > \frac{1}{2}\thinspace .
# $$

# A criterion against non-physical oscillations is therefore

# $$
# F_x + F_y \leq \frac{1}{2},
# $$

# which is the same limit as the stability criterion for the Forward Euler
# scheme.
# 
# The exact discrete solution is

# <!-- Equation labels as ordinary links -->
# <div id="diffu:2D:analysis:CN:numexact"></div>
# 
# $$
# \begin{equation}
# u^{n}_{q,r} = A
# \left(
# \frac{1 - 2(F_x\sin^2 p_x + F_x\sin^2p_y)}{1 + 2(F_x\sin^2 p_x + F_x\sin^2p_y)}
# \right)^n
# e^{i(k_xq\Delta x + k_yr\Delta y)}\thinspace .
# \label{diffu:2D:analysis:CN:numexact} \tag{26}
# \end{equation}
# $$

# ## Explanation of numerical artifacts
# 
# 
# The behavior of the solution generated by Forward Euler discretization in time (and centered
# differences in space) is summarized at the end of
# the section [diffu:pde1:FE:experiments](#diffu:pde1:FE:experiments). Can we, from the analysis
# above, explain the behavior?
# 
# We may start by looking at [Figure](#diffu:pde1:FE:fig:F=0.51)
# where $F=0.51$. The figure shows that the solution is unstable and
# grows in time. The stability limit for such growth is $F=0.5$ and
# since the $F$ in this simulation is slightly larger, growth is
# unavoidable.
# 
# [Figure](#diffu:pde1:FE:fig:F=0.5) has unexpected features:
# we would expect the solution of the diffusion equation to be
# smooth, but the graphs in [Figure](#diffu:pde1:FE:fig:F=0.5)
# contain non-smooth noise. Turning to [Figure](#diffu:pde1:FE:fig:gauss:F=0.5),  which has a quite similar
# initial condition, we see that the curves are indeed smooth.
# The problem with the results in [Figure](#diffu:pde1:FE:fig:F=0.5)
# is that the initial condition is discontinuous. To represent it, we
# need a significant amplitude on the shortest waves in the mesh.
# However, for $F=0.5$, the shortest wave ($p=\pi/2$) gives
# the amplitude in the numerical solution as $(1-4F)^n$, which oscillates
# between negative and positive values at subsequent time levels
# for $F>\frac{1}{4}$. Since the shortest waves have visible amplitudes in
# the solution profile, the oscillations becomes visible. The
# smooth initial condition in [Figure](#diffu:pde1:FE:fig:gauss:F=0.5),
# on the other hand, leads to very small amplitudes of the shortest waves.
# That these waves then oscillate in a non-physical way for
# $F=0.5$ is not a visible effect. The oscillations
# in time in the amplitude $(1-4F)^n$ disappear for $F\leq\frac{1}{4}$,
# and that is why also the discontinuous initial condition always leads to
# smooth solutions in [Figure](#diffu:pde1:FE:fig:F=0.25), where
# $F=\frac{1}{4}$.
# 
# Turning the attention to the Backward Euler scheme and the experiments
# in [Figure](#diffu:pde1:BE:fig:F=0.5), we see that even the discontinuous
# initial condition gives smooth solutions for $F=0.5$ (and in fact all other
# $F$ values). From the exact expression of the numerical amplitude,
# $(1  + 4F\sin^2p)^{-1}$, we realize that this factor can never flip between
# positive and negative values, and no instabilities can occur. The conclusion
# is that the Backward Euler scheme always produces smooth solutions.
# Also, the Backward Euler scheme guarantees that the solution cannot grow
# in time (unless we add a source term to the PDE, but that is meant to
# represent a physically relevant growth).
# 
# Finally, we have some small, strange artifacts when simulating the
# development of the initial plug profile with the Crank-Nicolson scheme,
# see [Figure](#diffu:pde1:CN:fig:F=10), where $F=3$.
# The Crank-Nicolson scheme cannot give growing amplitudes, but it may
# give oscillating amplitudes in time. The critical factor is
# $1 - 2F\sin^2p$, which for the shortest waves ($p=\pi/2$) indicates
# a stability limit $F=0.5$. With the discontinuous initial condition, we have
# enough amplitude on the shortest waves so their wrong behavior is visible,
# and this is what we see as small instabilities in
# [Figure](#diffu:pde1:CN:fig:F=10). The only remedy is to lower the $F$ value.
# 
# 
# 
# 
# # Exercises
# 
# 
# 
# <!-- --- begin exercise --- -->
# 
# ## Exercise 1: Explore symmetry in a 1D problem
# <div id="diffu:exer:1D:gaussian:symmetric"></div>
# 
# This exercise simulates the exact solution ([7](#diffu:pde1:sol:Gaussian)).
# Suppose for simplicity that $c=0$.
# 
# 
# **a)**
# Formulate an initial-boundary value problem that has
# ([7](#diffu:pde1:sol:Gaussian)) as solution in the domain $[-L,L]$.
# Use the exact solution ([7](#diffu:pde1:sol:Gaussian)) as Dirichlet
# condition at the boundaries.
# Simulate the diffusion of the Gaussian peak. Observe that the
# solution is symmetric around $x=0$.
# 
# **b)**
# Show from ([7](#diffu:pde1:sol:Gaussian)) that $u_x(c,t)=0$.
# Since the solution is symmetric around $x=c=0$, we can solve the
# numerical problem in frac{1}{2} of the domain, using a *symmetry boundary condition*
# $u_x=0$ at $x=0$. Set up the
# initial-boundary value problem in this case. Simulate the
# diffusion problem in $[0,L]$ and compare with the solution in a).
# 
# 
# <!-- --- begin solution of exercise --- -->
# **Solution.**

# $$
# \begin{align*}
# u_t &= \dfc u_xx,\\ 
# u_x(0,t) &= 0,\\ 
# u(L,t)& =\frac{1}{\sqrt{4\pi\dfc t}} \exp{\left({-\frac{x^2}{4\dfc t}}\right)}\thinspace .
# \end{align*}
# $$

# <!-- --- end solution of exercise --- -->
# 
# Filename: `diffu_symmetric_gaussian`.
# 
# <!-- --- end exercise --- -->
# 
# 
# 
# 
# <!-- --- begin exercise --- -->
# 
# ## Exercise 2: Investigate approximation errors from a $u_x=0$ boundary condition
# <div id="diffu:exer:1D:ux:onesided"></div>
# 
# We consider the problem solved in [Exercise 1: Explore symmetry in a 1D problem](#diffu:exer:1D:gaussian:symmetric)
# part b). The boundary condition $u_x(0,t)=0$ can be implemented in
# two ways: 1) by a standard symmetric finite difference $[D_{2x}u]_i^n=0$,
# or 2) by a one-sided difference $[D^+u=0]^n_i=0$.
# Investigate the effect of these two conditions on the
# convergence rate in space.
# 
# <!-- --- begin hint in exercise --- -->
# 
# **Hint.**
# If you use a Forward Euler scheme, choose a discretization parameter
# $h=\Delta t = \Delta x^2$ and assume the error goes like $E\sim h^r$.
# The error in the scheme is $\Oof{\Delta t,\Delta x^2}$ so one should
# expect that the estimated $r$ approaches 1. The question is if
# a one-sided difference approximation to $u_x(0,t)=0$ destroys this
# convergence rate.
# 
# <!-- --- end hint in exercise --- -->
# Filename: `diffu_onesided_fd`.
# 
# <!-- --- end exercise --- -->
# 
# 
# 
# 
# <!-- --- begin exercise --- -->
# 
# ## Exercise 3: Experiment with open boundary conditions in 1D
# <div id="diffu:exer:1D:openBC"></div>
# 
# We address diffusion of a Gaussian function
# as in [Exercise 1: Explore symmetry in a 1D problem](#diffu:exer:1D:gaussian:symmetric),
# in the domain $[0,L]$,
# but now we shall explore different types of boundary
# conditions on $x=L$. In real-life problems we do not know
# the exact solution on $x=L$ and must use something simpler.
# 
# 
# **a)**
# Imagine that we want to solve the problem numerically on
# $[0,L]$, with a symmetry boundary condition $u_x=0$ at $x=0$,
# but we do not know the exact solution and cannot of that
# reason assign a correct Dirichlet condition at $x=L$.
# One idea is to simply set $u(L,t)=0$ since this will be an
# accurate approximation before the diffused pulse reaches $x=L$
# and even thereafter it might be a satisfactory condition if the exact $u$ has
# a small value.
# Let $\uex$ be the exact solution and let $u$ be the solution
# of $u_t=\dfc u_{xx}$ with an initial Gaussian pulse and
# the boundary conditions $u_x(0,t)=u(L,t)=0$. Derive a diffusion
# problem for the error $e=\uex - u$. Solve this problem
# numerically using an exact Dirichlet condition at $x=L$.
# Animate the evolution of the error and make a curve plot of
# the error measure

# $$
# E(t)=\sqrt{\frac{\int_0^L e^2dx}{\int_0^L udx}}\thinspace .
# $$

# Is this a suitable error measure for the present problem?
# 
# **b)**
# Instead of using $u(L,t)=0$ as approximate boundary condition for
# letting the diffused Gaussian pulse move out of our finite domain,
# one may try $u_x(L,t)=0$ since the solution for large $t$ is
# quite flat. Argue that this condition gives a completely wrong
# asymptotic solution as $t\rightarrow 0$. To do this,
# integrate the diffusion equation from $0$ to $L$, integrate
# $u_{xx}$ by parts (or use Gauss' divergence theorem in 1D) to
# arrive at the important property

# $$
# \frac{d}{dt}\int_{0}^L u(x,t)dx = 0,
# $$

# implying that $\int_0^Ludx$ must be constant in time, and therefore

# $$
# \int_{0}^L u(x,t)dx = \int_{0}^LI(x)dx\thinspace .
# $$

# The integral of the initial pulse is 1.
# 
# **c)**
# Another idea for an artificial boundary condition at $x=L$
# is to use a cooling law

# <!-- Equation labels as ordinary links -->
# <div id="diffu:pde1:Gaussian:xL:cooling"></div>
# 
# $$
# \begin{equation}
# -\dfc u_x = q(u - u_S),
# \label{diffu:pde1:Gaussian:xL:cooling} \tag{27}
# \end{equation}
# $$

# where $q$ is an unknown heat transfer coefficient and $u_S$ is
# the surrounding temperature in the medium outside of $[0,L]$.
# (Note that arguing that $u_S$ is approximately $u(L,t)$ gives
# the $u_x=0$ condition from the previous subexercise that is
# qualitatively wrong for large $t$.)
# Develop a diffusion problem for the error in the solution using
# ([27](#diffu:pde1:Gaussian:xL:cooling)) as boundary condition.
# Assume one can take $u_S=0$ "outside the domain" since
# $\uex\rightarrow 0$ as $x\rightarrow\infty$.
# Find a function $q=q(t)$ such that the exact solution
# obeys the condition ([27](#diffu:pde1:Gaussian:xL:cooling)).
# Test some constant values of $q$ and animate how the corresponding
# error function behaves. Also compute $E(t)$ curves as defined above.
# 
# 
# Filename: `diffu_open_BC`.
# 
# <!-- --- end exercise --- -->
# 
# 
# 
# 
# <!-- --- begin exercise --- -->
# 
# ## Exercise 4: Simulate a diffused Gaussian peak in 2D/3D
# 
# 
# **a)**
# Generalize ([7](#diffu:pde1:sol:Gaussian)) to multi dimensions by
# assuming that one-dimensional solutions can be multiplied to solve
# $u_t = \dfc\nabla^2 u$. Set $c=0$ such that the peak of
# the Gaussian is at the origin.
# 
# **b)**
# One can from the exact solution show
# that $u_x=0$ on $x=0$, $u_y=0$ on $y=0$, and $u_z=0$ on $z=0$.
# The approximately correct condition $u=0$ can be set
# on the remaining boundaries (say $x=L$, $y=L$, $z=L$), cf. [Exercise 3: Experiment with open boundary conditions in 1D](#diffu:exer:1D:openBC).
# Simulate a 2D case and make an animation of the diffused Gaussian peak.
# 
# **c)**
# The formulation in b) makes use of symmetry of the solution such that we
# can solve the problem in the first quadrant (2D) or octant (3D) only.
# To check that the symmetry assumption is correct, formulate the problem
# without symmetry in a domain $[-L,L]\times [L,L]$ in 2D. Use $u=0$ as
# approximately correct boundary condition. Simulate the same case as
# in b), but in a four times as large domain. Make an animation and compare
# it with the one in b).
# 
# Filename: `diffu_symmetric_gaussian_2D`.
# 
# <!-- --- end exercise --- -->
# 
# 
# 
# 
# <!-- --- begin exercise --- -->
# 
# ## Exercise 5: Examine stability of a diffusion model with a source term
# <div id="diffu:exer:uterm"></div>
# 
# Consider a diffusion equation with a linear $u$ term:

# $$
# u_t = \dfc u_{xx} + \beta u\thinspace .
# $$

# **a)**
# Derive in detail the Forward Euler, Backward Euler,
# and Crank-Nicolson schemes for this type of diffusion model.
# Thereafter, formulate a $\theta$-rule to summarize the three schemes.
# 
# **b)**
# Assume a solution like ([8](#diffu:pde1:sol1)) and find the relation
# between $a$, $k$, $\dfc$, and $\beta$.
# 
# <!-- --- begin hint in exercise --- -->
# 
# **Hint.**
# Insert ([8](#diffu:pde1:sol1)) in the PDE problem.
# 
# <!-- --- end hint in exercise --- -->
# 
# **c)**
# Calculate the stability of the Forward Euler scheme. Design
# numerical experiments to confirm the results.
# 
# <!-- --- begin hint in exercise --- -->
# 
# **Hint.**
# Insert the discrete counterpart to ([8](#diffu:pde1:sol1)) in the
# numerical scheme. Run experiments at the stability limit and slightly above.
# 
# <!-- --- end hint in exercise --- -->
# 
# **d)**
# Repeat c) for the Backward Euler scheme.
# 
# **e)**
# Repeat c) for the Crank-Nicolson scheme.
# 
# **f)**
# How does the extra term $bu$ impact the accuracy of the three schemes?
# 
# <!-- --- begin hint in exercise --- -->
# 
# **Hint.**
# For analysis of the accuracy,
# compare the numerical and exact amplification factors, in
# graphs and/or by Taylor series expansion.
# 
# <!-- --- end hint in exercise --- -->
# 
# Filename: `diffu_stability_uterm`.
# 
# <!-- --- end exercise --- -->