\magnification=\magstep1
\vsize=9truein
\hsize=6.5truein
%\parindent=0.8truein
\def\C{{\bf C}}
\def\R{{\bf R}}
\def\N{{\bf N}}
\def\span{{\rm span}}
\font\sc=cmcsc10
\font\smalltt=cmtt10 scaled840
\font\smseven=cmr7 scaled840

\newcount\qnum\newcount\qprt\qnum=0\qprt=0
\newcount\ano
\def\qn{\ano=0\global\qprt=0\global\advance\qnum by1\medskip
    \par\noindent\item{\bf\the\qnum.}}
\def\cqnn{\ano=0\global\advance\qprt by1\par
    \par\noindent\itemitem{\bf(\romannumeral\qprt)}}
\def\qnn{\ano=0\global\advance\qprt by1\medskip
    \par\noindent\itemitem{\bf(\romannumeral\qprt)}}

%verbatim environment
\def\uncatcodespecials{\def\do##1{\catcode`##1=12 }\dospecials}
\newcount\lineno
\def\listing#1{\par\begingroup\setupverbatim #1\endgroup}
\def\setupverbatimc{\parindent=1truecm\tt
%   \smalltt\baselineskip=10pt plus 1pt minus 1pt
    \def\par{\leavevmode\endgraf}\catcode`\`=\active\obeylines
    \uncatcodespecials\obeyspaces
    \everypar{\global\advance\lineno by 1\llap{\smseven\the\lineno\ \ }}}
\def\setupverbatimb{\parindent=1truecm
    \global\lineno=0\smalltt\baselineskip=10pt plus 1pt minus 1pt
    \def\par{\leavevmode\endgraf}\catcode`\`=\active\obeylines
    \uncatcodespecials\obeyspaces
    \everypar{\advance\lineno by 1{}}}
\def\setupverbatim{\global\lineno=0\setupverbatimc}

{\obeyspaces\global\let =\ }

\def\beginverbatim{\par\begingroup\setupverbatim\doverbatim}
{\catcode`\|=0\catcode`\\=12\catcode`|`=|active|gdef`{|relax|lq}%
|obeylines|gdef|doverbatim^^M#1\endverbatim{#1|endgroup}}

\def\beginverbatimb{\par\begingroup\setupverbatimb\doverbatim}
\def\listing#1{\medskip
    \begingroup\setupverbatim\input#1\endgroup\medskip
    \noindent\ignorespaces}
\def\listout#1{\medskip
    \begingroup\setupverbatimb\input#1\endgroup\medskip
    \noindent\ignorespaces}
\def\macinput#1{{\catcode`\$=10\catcode`@=3 \input #1 }}

\input epsf

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\qn [Thomas Problem 1.2.1] \ Write a code to approximately
solve the boundary value problem
$$\left\{\eqalign{
	v_t&=\nu v_{xx},\quad x\in (0,1),\quad t>0\cr
	v(x,0)&=f(x),\quad x\in [0,1]\cr
	v(0,t)&=a(t),\quad t\ge 0\cr
	v(1,t)&=b(t),\quad t\ge 0\cr
}\right.$$
Use $f(x)=\sin 2\pi x$, $a(t)=b(t)=0$, $K=10$ and $\nu=1/6$.  Find
solutions at $t=0.06$, $t=0.1$, $t=0.9$ and $t=50.0$.  For the first
three values of $t$, use $\Delta t=0.02$.  To speed the solution of
the last value of $t$, you might use a larger value for $\Delta t$.
Determine how large you can choose $\Delta t$ and still get results
that might still be correct.  Compare and contrast your solutions
to the exact solutions.
\bigskip
\noindent
We compute the exact solution using separation of variables.
First look for eigenfunctions $X(x)$ of the 
operator $\partial^2/\partial x^2$ such that $X(0)=X(1)=0$.
Solving
$$X''(x)=\lambda X(x)$$ 
yields
$$X(x)=
A\exp(\sqrt{\lambda} x)+B\exp(-\sqrt{\lambda} x).$$
Now $X(0)=0$ implies $B=-A$ so that
$$X(t)=
A\exp(\sqrt{\lambda} x)-A\exp(-\sqrt{\lambda} x),$$
and on the other hand $X(1)=0$ 
implies $\exp(\sqrt{\lambda})=\exp(-\sqrt{\lambda})$.
Therefore $\exp(2\sqrt{\lambda})=1$ implies
$2\sqrt{\lambda}=2n \pi i$ and
we obtain eigenfunctions $X_n(x)$ with eigenvalues $\lambda_n$
given by
$$
	X_n(x)=\sin(n\pi x)\qquad\hbox{and}\qquad
\lambda_n=-n^2\pi^2\qquad\hbox{for}\qquad
n\in\N.$$

Consider a solution of the form
$$v(x,t)=\sum_{n\in\N} T_n(t)X_n(x).$$
Plugging this into the differential equation obtains
$$
	\sum_{n\in\N} T_n'(t)X_n(x)= \sum_{n\in\N} \nu\lambda_n T_n(t)X_n(x)
$$
and since the eigenfunctions $X_n(x)$ are linearly independent, it follows
that
$$
	T_n'(t)=\nu\lambda_n T_n(t)\qquad\hbox{for all}\qquad n\in\N.$$
Solving this system yields that
$ T_n(t)=T_n(0)\exp(-\nu n^2\pi^2 t) $
and consequently
$$v(x,t)=\sum_{n\in\N} T_n(0)\exp(-\nu n^2\pi^2 t)\sin(n\pi x).$$

The initial condition
$$
    v(x,0)=\sum_{n\in\N} T_n(0)\sin(n\pi x)=f(x)=\sin(2\pi x)
$$
implies
$$
	T_n(0)=\cases {1& if $n=2$,\cr
			0&otherwise.}
$$
Therefore, the exact solution is
$$
	v(x,t)=
	\exp(-\nu 4\pi^2t)\sin(2\pi x).
$$

We use the finite difference scheme
$$
\cases{
	u_k^{n+1}=u_k^n+\rho\delta^2 u_k^n&for 
		$n=0,\ldots,N-1$ and $k=1,\dots,K-1$\cr
	u_k^0=f(x_k)&for $k=0,\ldots,K$\cr
	u_0^n=a(t_n)&for $n=1,\ldots,N$\cr
	u_K^n=b(t_n)&for $n=1,\ldots,N$\cr
}$$
where $x_k=k\Delta x$, $t_n=n\Delta t$, $\Delta x=1/K$, $\Delta t=T/N$
and $\rho=\nu\Delta t/\Delta x^2$.  We chose the values of $K$, $N$ and 
$T$ given in the table
$$
\vbox{\halign{\quad \hfil#\quad&&\quad\hfil #\quad\cr
	$K$&$N$&$T$&$\Delta x$&$\Delta t$&$\rho$\cr
\noalign{\smallskip\hrule\smallskip}
	10&3&0.06&0.1&0.02&0.333333\cr
	10&5&0.1&0.1&0.02&0.333333\cr
	10&45&0.9&0.1&0.02&0.333333\cr
	10&2000&50&0.1&0.025&0.416667\cr
	10&1750&50&0.1&0.0285714&0.47619\cr
	10&1500&50&0.1&0.0333333&0.555556\cr
}}
$$
In the following program, and all other programs included in this
homework assignment, the output will be displayed in four columns.
This output describes the approximation obtained by the numerical
scheme, the exact solution and the error between the approximation
and the exact solution.  Specifically we have
\medskip
\item{1.} The column{\tt\ xk }gives the value of $x_k$.
\smallskip
\item{2.} The column{\tt\ uk }gives the approximation $u_k^N$ for $v(x_k,T)$.
\smallskip
\item{3.} The column{\tt\ vk }gives the exact value of $v(x_k,T)$.
\smallskip
\item{4.} The column{\tt\ vk-uk }gives the error $u_k^N-v(x_k,T)$.
\medskip\noindent
The program is
\medskip
\listing{prog1.c}
\noindent
The output is
\listout{prog1.out}
\noindent
Note that the last computation is clearly wrong because rather
than decaying, the solution increases.  In general, the scheme
is unstable and therefore unreliable when $\rho>1/2$.
\vfill
\eject

\qn [Thomas Problem 1.3.1] \  Solve Problem 1.2.1 using the leapfrog scheme
$$
	u_k^{n+1}=u_k^{n-1}+{2\nu\Delta t\over \Delta x^2}
	\big(u_{k+1}^n-2u_k^n+u_{k-1}^n\big)
$$
with $\Delta t=0.02$.  For convenience, use values from the exact
solution at $\Delta t$ to get the leapfrog scheme started.
\medskip
\noindent
The program is
\listing{prog2.c}
\noindent
The output is
\listout{prog2.out}
\noindent
Note that this method oscillates between $\pm \infty$ at $T=50$ which is
a noticably poor result as the solution is supposed to converge to 0.
\vfill
\eject

\qn [Thomas Problem 1.3.2] \  Use $\Delta t=\Delta x^2/(6\nu)$
for the solution to Problem 1.2.1.
\medskip
\noindent
The program is
\listing{prog3.c}
\noindent
The output is
\listout{prog3.out}
\noindent
Note that this choice of $\Delta t$ results in significantly more 
accurate accurate 
solutions than the choices used in Problem 1.2.1.
\vfill
\eject

\qn [Thomas Problem 1.5.1] \ Solve the problem
$$\left\{\eqalign{
	v_t+av_x&=\nu v_{xx},\quad x\in (0,1),\quad t>0\cr
	v(x,0)&=f(x),\quad x\in [0,1]\cr
	v(0,t)&=v(1,t)=0,\quad t\ge 0\cr
}\right.$$
using the difference scheme
$$\cases{
	u_k^{n+1}=u_k^n-R\delta_0 u_k^n+\rho\delta^2 u_k^n&
		for $n=0,\ldots,N$ and $k=1,\ldots,K-1$\cr
	u_k^0=f(x_k)& for $k=0,\ldots,K$\cr
	u_0^n=u_K^n=0& for $n=1,\ldots,N$\cr
}$$
where $R=a\Delta t/(2\Delta x)$ and $\rho=\nu\Delta t/\Delta x^2$.
Use $\nu=1.0$, $a=2$, $K=20$ and $f(x)=\sin 4\pi x$.  Find
solutions at $t=0.06$, $t=0.1$, $t=0.9$ using $\Delta t=0.001.$
\medskip
\noindent
We first find the exact solution.  We begin by
finding eigenfunctions $X(x)$ so that
$$\nu X''(x) -a X'(x)=\lambda X(x)
\qquad\hbox{and}\qquad X(0)=X(1)=0.$$
The characteristic equation is
$$
	\nu r^2-ar-\lambda=0
$$
from which we determine that
$$
	r_1={a+\sqrt{a^2+4\nu\lambda}\over 2\nu}
\qquad\hbox{and}\qquad
	r_2={a-\sqrt{a^2+4\nu\lambda}\over 2\nu}.
$$
It follows that
$$
	X(x)=A\exp(r_1 x)+B\exp(r_2 x).
$$
Since $X(0)=0$ we obtain $B=-A$.  Consequently
$$
	X(x)=A\exp(r_1 x)-A\exp(r_2 x).
$$
Now $X(1)=0$ implies
$ \exp(r_1)=\exp(r_2) $ or that $\exp(r_1-r_2)=1$.
It follows that
$$
	r_1-r_2={\sqrt{a^2+4\nu\lambda}\over \nu}=2n\pi i
$$
and so we obtain eigenfunction $X_n(x)$ with eigenvalues $\lambda_n$ 
given by 
$$
	X_n(x)=\exp(ax/2)\sin(n\pi x)\quad\hbox{and}\quad
	\lambda_n=-\nu n^2\pi^2-{a^2\over 4\nu}\quad\hbox{for}\quad n\in \N.
$$

Consider a solution of the form
$$v(x,t)=\sum_{n\in\N} T_n(t)X_n(x).$$
Plugging this into the differential equation obtains
$$
	\sum_{n\in\N} T_n'(t)X_n(x)=
\sum_{n\in\N} \lambda_n T_n(t)X_n(x)
$$
and since the eigenfunctions $X_n(x)$ are linearly independent, it follows
that
$$
	T_n'(t)=\lambda_n T_n(t)\qquad\hbox{for all}\qquad n\in\N.$$
Solving this system yields that
$$ T_n(t)=T_n(0)
\exp\Big(-\nu n^2\pi^2t-{a^2 t\over 4\nu} \Big)
 $$
and consequently
$$v(x,t)=\sum_{n\in\N} T_n(0)
\exp\Big(-\nu n^2\pi^2t-{a^2 t\over 4\nu} \Big)
\exp(ax/2)
	\sin(n\pi x).$$

The initial condition
$$
    v(x,0)=\sum_{n\in\N} T_n(0)\exp(ax/2)\sin(n\pi x)=f(x)=\sin(4\pi x)
$$
implies
$$
    \sum_{n\in\N} T_n(0)\sin(n\pi x)=\exp(-ax/2)\sin(4\pi x)
$$
Multiplying by $\sin(m\pi x)$ and integrating yields that
$$\eqalign{
	T_m(0)&=2\int_0^1 \exp(-ax/2)\sin(4\pi x)\sin(m\pi x)\,dx\cr
			&= {128 a m (-1+(-1)^m\exp(-a/2))\over
		16m^4\pi^4+(8a^2\pi^2-512\pi^4)m^2+a^4+128a^2\pi^2+4096\pi^4}.
}$$
Therefore, the exact solution is a complicated series.  We turn to
Maple to help with this series.  The maple script
\listing{maple4.mpl}
produces the output
\listout{maple4.out}
We replace{\tt\ pow(-1.0,1.0*m) }with{\tt\ (1-2*(m\%2)) }to obtain
the exact solution which, for convenience, we keep in a
separate file of utility programs
\listing{util4.c}
\noindent
with headers
\listing{util4.h}
\noindent
The main routine is
\listing{prog4.c}
\noindent
The output is
\listout{prog4.out}
\vfill
\eject

\qn [Thomas Problem 1.5.2] \ Repeat Problem 1.5.1 with $\nu=0.01$.
\medskip\noindent
The program is
\listing{prog5.c}
\noindent
The output is
\listout{prog5.out}
\vfill
\eject

\qn [Thomas Problem 1.5.3] \ Repeat Problem 1.5.1 and 1.5.2 with K=40.
\medskip\noindent
The program is
\listing{prog6.c}
\noindent
The output is
\listout{prog6.out}
\noindent
Note that in order to verify the answers with $K=40$ one needs also
to decrease the $\Delta t$ time steps.
\vfill
\eject

\qn [Thomas Problem 1.5.7] \ Solve the following initial boundary
value problem
$$\cases{
	v_t=\nu v_{xx}+F(x,t),&for $x\in(0,1)$ and $t>0$\cr
	v(x,0)=f(x),&for $x\in[0,1]$\cr
	v(0,t)=a(t),&for $t\ge 0$\cr
	v(1,t)=b(t),&for $t\ge 0$\cr
}$$
where $\nu=0.1$, $f(x)=0$, $a(t)=b(t)=0$ 
and $F(x,t)=\sin(2\pi x)\sin(4\pi t)$.
Use $K=10$ and $\Delta t=0.05$ and find solutions
at $t=0.1$, $t=0.9$ and $t=2.0$.  Resolve the problem
in the above choice of parameters using $K=20$ and $\Delta t=0.01$.
\medskip\noindent
First we derive the exact solution.  Following the techniques used
in Problem 1.2.1 let $X_n(x)=\sin(n\pi x)$ and $\lambda_n=-n^2\pi^2$
and write
$$
	v(x,t)=\sum_{n\in\N} T_n(t)X_n(x).
$$
Plugging this into the differential equations yields
$$
	\sum_{n\in\N} T'_n(t)X_n(x)=\sum_{n\in\N}\nu\lambda_n T_n(t)X_n(x)+
		\sin(4\pi t)X_2(x).
$$
Therefore
$$\cases{
	T'_n(t)=\nu\lambda_n T_n(t)&for $n\ne 2$\cr
\noalign{\smallskip}
	T'_2(t)=\nu\lambda_2 T_2(t)+\sin(4\pi t)&for $n=2$.
}$$
Solution for $T_n(t)$ when $n\ne 2$ is the same as in Problem 1.2.1.
When $n=2$ we procede as follows:
$$
	{d\over dt}\big( T_2(t) \exp(-\nu\lambda_2t)\big)
		=\sin(4\pi t)\exp(-\nu\lambda_2t).
$$
Therefore
$$\eqalign{
	T_2(t)\exp(-\nu\lambda_2 t)-T_2(0)&=
		\int_0^t \sin(4\pi s)\exp(-\nu\lambda_2 s)\,ds\cr
%	&={4\pi-(4\pi\cos(4\pi t)+\nu\lambda_2\sin(4\pi t))
%		\exp(-\nu\lambda_2 t)\over
%		\nu^2\lambda_2^2+16\pi^2}\cr
	&={1-(\cos(4\pi t)-\nu\pi\sin(4\pi t))\exp(\nu 4\pi^2 t)
		\over 4\pi(\nu^2\pi^2+1)}.\cr
}$$
It follows that
$$\cases{
	T_n(t)=
T_n(0)\exp(-\nu n^2\pi^2t)&for $n\ne 2$\cr
\noalign{\medskip}
	T_2(t)=
T_2(0)\exp(-\nu 4\pi^2t)
+\displaystyle {\exp(-\nu 4\pi^2 t)-\cos(4\pi t)+\nu\pi\sin(4\pi t)
		\over 4\pi (\nu^2\pi^2+1)}& for $n=2$.\cr
}$$
Solving for the initial condition 
$$
	\sum_{n\in\N} T_n(0)\sin(n\pi x)=0.
$$
yields than $T_n(0)=0$ for all $n$.  Therefore the solution is
$$
	v(x,t)={\exp(-\nu 4\pi^2 t)-\cos(4\pi t)+\nu\pi\sin(4\pi t)
        \over 4\pi (\nu^2\pi^2+1)} \sin(2\pi x).
$$

Integrating both sides of the equation from $t_n$ to $t_{n+1}$ yields
$$
    v(x,t_{n+1})
	=v(x,t_n)+\int_{t_n}^{t_{n+1}} \big(\nu v_{xx}(x,s)+F(x,s)\big)\,ds
$$
Approximating integral using the left endpoint yields
$$
	v(x,t_{n+1})\approx v(x,t_n)+\Delta t \big(\nu v_{xx}(x,t_n)+F(x,t_n)\big)
$$
and the explicit method
$$
	u_k^{n+1}=u_k^n+R \delta^2 u_k^n + F(x_k,t_n).
$$
Alternatively, for better accuracy, we can approximate the second term 
in the intergral using the midpoint method.  This yields an
explicit scheme of the form
$$
	u_k^{n+1}=u_k^n+R \delta^2 u_k^n + F(x_k,t_{n+1/2}).
$$
The program is
\listing{prog7.c}
\noindent
The output is
\listout{prog7.out}
\noindent
Note that our numeric output verifies that the midpoint 
is generally more accurate than the left endpoint for
approximating the contribution of the force $F(x,t)$.
\vfill
\eject
\end
