% maple lab 2 S2006
\documentclass[leqno]{article}
\input{linalg.sty}
\def\rf#1{(\ref{#1})}
\def\BF#1{{\bf #1}}
\def\Prob#1#2{\par\bigskip\par\noindent
{\Large\sf Problem #1}. {\bf #2}\par}
\parskip=3pt
\def\MAPLE{{\tt maple}}
% \begin{topics}{template} ... \end{topics}
% used like enumerate, itemize, description environments
% The item[] is never optional!
%
\iffalse
\newcommand{\topicslabel}[1]{\mbox{#1}\hfil}
\newenvironment{topics}[1]{
    \begin{list}{}{%
      \let\makelabel\topicslabel
      \settowidth{\labelwidth}{\topicslabel{#1}}
      \setlength{\leftmargin}{\labelwidth}
      \advance\leftmargin by 1.5em}
}{\end{list}}
\fi
\pagestyle{plain}
\setlength{\textheight}{9.8in}
\begin{document}
\thispagestyle{empty}

 {\bf Name~}\underline{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}\hfill
 {Class Time~}\underline{~~~~~~~~}\newline


\begin{center}
\bf Math 2280
\bf Maple Project 2: Newton Cooling \\
{\small\bf January 2006}
\end{center}


{\bf Due date}: See the internet due dates.

{\bf References}: Edwards-Penney, pages 55--57. Code in \MAPLE{}
appears in \verb|2280mapleL2-S2006.txt| at URL\newline
\verb|http://www.math.utah.edu/~gustafso/|. This document:
\verb|2280mapleL2-S2006.pdf|.


{\bf The Problem and Model}. The project is to determine an
approximation for the indoor temperature $u(t)$ in an unheated
building. The model uses Newton's cooling law, insulation data $k$
and a formula for the ambient outside temperature $A(t)$ (see the
background section {\em infra}):
\begin{equation}\label{newtoncoolingone}
\begin{array}{l} u'(t)+ku(t)=kA(t),\\
  u(0)=u_0. \end{array}
\end{equation}

{\bf Assumptions and Notation}. Let the daily temperature in Salt
Lake City vary from a low of $m=22$ degrees to a high of $M=48$
degrees with the low and high at $3$am and $3$pm, respectively.
The building is assumed insulated, the average insulation factor
being represented by a Newton cooling constant $k$, usually $0.2
\le k \le 0.5$. The ambient outside temperature is called $A(t)$
and the inside temperature is called $u(t)$, where $t$ is in hours
measured from midnight $t=0$. The number $u_0$ is the temperature
at midnight. The model assumes that no heat sources or heat sinks
are present inside the building. Symbol $u_h(t)$ refers to the
general solution of the homogeneous equation $u'+ku=0$. Symbol
$u_p(t)$ stands for a particular solution of $u'+ku=kA(t)$.


\Prob{L2.1}{(Solution Formulae for $u_p$ and $u$)} Derive by hand,
using \MAPLE{} integral tables, an explicit symbolic formula for
$u_p(t)$. Display a final formula for $u=u_h+u_p$ which depends
only on $t$, $u_0$, $k$ and $\omega$. Check your hand answer for
$u$ in \MAPLE{}. The only \MAPLE{} assist in this problem is the
answer check and the \MAPLE{} integral table.

\Prob{L2.2}{(Steady-State Periodic Solution)} Derive by hand a
formula for the steady-state periodic solution $u_{\mbox{ss}}$ of
$u'+ku=kA(t)$. There is no \MAPLE{} assist in this problem, it is
all done by hand from the answer in L2.1. The results:
$$
\begin{array}{rcl}
u_{\mbox{ss}}(t) &=& 35 -\frac{13k}{k^2+\omega_0^2}(k\cos\omega_0
(t-3)+\omega_0\sin\omega_0 (t-3)), \\
 u(t) &=& Ce^{-kt} + u_{\mbox{ss}}(t),\\
C &=&u_0-u_{\mbox{ss}}(0). \end{array}
$$

\Prob{L2.3}{(Indoor-Outdoor Variation)} Compare in a \MAPLE{}
graphic the indoor ($u(t)$) and outdoor ($A(t)$) temperature
oscillations over a $48$-hour period assuming $k=0.35$, $u_0=74$,
$\omega=\pi/12$. Compute the indoor and outdoor temperature
variations from this $48$-hour plot. Find the phase delay from
this same plot.

\Prob{L2.4}{(Freezing Pipes)} Assume $\omega=\pi/12$ and the
insulation constant $k$ ranges from $0.2$ to $0.4$. Suppose the
inside temperature is $74$ degrees at midnight when the furnace is
turned off. Report approximate ranges of hours and insulation
constants, during the first 72 hours and $0.2 < k < 0.4$, for
which the indoor temperature is at or below $29$ degrees. Justify
your logic used to find the ranges, in a short paragraph.
Illustrate with a computer graphic. \vfill

{{\sf Staple this page on top of all hand and maple work.}
\newline
\medskip
Sample maple code, problem notes and background are on the next
page \ldots}
\newpage

\section*{Background}
{\bf Newton Cooling Model}.
Newton's law of cooling is:
\begin{quote}\em
The rate of change
of the indoor temperature is proportional to the difference between the
ambient and indoor temperatures.
\end{quote}
This implies that $du/dt$ is proportional to $A-u$. Let $A(t)$ be
the ambient outside temperature and let $k$ denote the
proportionality constant in Newton's cooling law. A model for the
indoor temperature $u(t)$ is given by $du/dt=k(A-u)$, $u(0)=u_0$,
which can be rearranged as
\begin{equation}\label{newtoncooling}
u'(t)+ku(t)=kA(t),\quad u(0)=u_0.
\end{equation}
The number $k$ is called the {\bf insulation
constant}.

{\bf Ambient Temperature Model}. Let $M=48$, $m=22$,
$\omega_0=\pi/12$. The ambient temperature formula
$$A(t)=\frac12(M+m)-\frac12(M-m)\cos \omega_0 (t-3)$$
satisfies $M=\max A(t)=A(15)$, $m=\min A(t)=A(3)$
and $A(t)$ is $t$-periodic of period $24$ hours.

For use in \MAPLE{}, a function of two variables
$$AA(t,\omega)=\frac12(M+m)-\frac12(M-m)\cos \omega (t-3)$$
is used to maintain the variable name $\omega$ in displays.

{\bf Indoor Temperature $u(t)$}.
The integrating factor method for linear equations
applies to find the general solution by these steps:

\def\AA{XX$u=u_0e^{-kt}+e^{-kt}\int_0^t kA(r)e^{kr}\,dr$}
\begin{topics}{\AA}
\item[$u'+ku=kA(t)$]  Copy the differential equation.
The integrating factor is $e^{kt}$.

\item[$\dd\frac{(e^{kt}u)'}{e^{kt}}=kA(t)$]  Replace the left
side.

\item[$(e^{kt}u)'=kA(t)e^{kt}$]  Clear fractions.

\item[$e^{kt}u=u_0+\int_0^t kA(x)e^{kx}\,dx$]
Integrate both sides with respect to $t$. Apply the Fundamental Theorem of Calculus. Use
$u(0)=u_0$.

\item[$u=u_0e^{-kt}+e^{-kt}\int_0^t kA(x)e^{kx}\,dx$]
Divide to isolate the {\em explicit solution} $u$.

\end{topics}

Let $u_h(t)=u_0e^{-kt}$, a solution of the homogeneous
differential equation $u'+ku=0$. Let $u_p(t)=ke^{-kt}\int_0^t
e^{kx}A(x)dx$, a particular solution of the nonhomogeneous
differential equation $u'+ku=kA(t)$. Then the indoor temperature
$u = u_h + u_p$ depends on the time $t$, the initial temperature
$u_0$, the insulation constant $k$ and the frequency $\omega$
(which is fixed at $\pi/12$). Write $u=u(t,u_0,k,\omega)$ to
emphasize the dependence. In \MAPLE{}, advantages exist for adding
the variable name $\omega$, which is later set to value
$\omega_0=\pi/12$. Write $u$ as \fbox{\tt U(t,u0,k,omega)} for use
in \MAPLE{}.

{\bf Steady-state solution}. The steady-state solution
$u_{\mbox{ss}}$ is an expression obtained from the general
solution formula $u=u_h+u_p$ by dropping all terms containing a
negative exponential. It depends on $t$, $k$ and $\omega$. This
expression is known to be a solution of the differential equation,
independent of any arbitrary constants, which represents the {\bf
observed temperature} after a long time.  In short,
$u_{\mbox{ss}}$ is independent of $u_0$.


%\noindent
%\begin{tabular}{p{3.35in}} \hline\\
%\end{tabular}
%\smallskip

\newpage

\section*{Problem Notes}



\subsection*{Notes on L2.1}

The integration problem can be solved by hand using the book's
integral table:
$$u_p(t)=ke^{-kt}\int_0^t e^{kx}\left(35-13\cos \left(\omega(x-3)\right)\right)dx.$$
A change of variable $u=x-3$ allows use of integral table entry $\int
e^{au}\cos bu\,du$. The symbol $\omega$ will be set to $\pi/12$, but for
simplicity, use symbol $\omega$ throughout.

A better integral table method is {\tt maple}. It uses a table
lookup method superior to the hand method above, in that no change
of variable is required to obtain the answer. Most people will use
this method because it is accurate and avoids primary pitfalls in
using a book table. A coding example appears below, which moves
the factor $ke^{-kt}$ inside the integrand, in order to
automatically perform simplifications.
\begin{verbatim}
 # maple integral table lookup
 unassign('k','omega','t','x');
 integrand:=(35-13*cos(omega*(x-3)))*k*exp(k*x-k*t);
 F:=unapply(integrand,(x,t,k,omega));
 int(F(x,t,k,omega),x=0..t);
\end{verbatim}
The formulae for $u_h$ and $u_p$ are used again in L2.2 in order
to derive the steady-state solution. Your answer {\em must
contain} symbols \fbox{\tt t, u0, k, omega}.

The answer check in \MAPLE{} is organized as follows. The complications
of setting $\omega=\pi/12$ are avoided here by leaving $\omega$ as a
symbol, since it does not affect the answer check.
\begin{verbatim}
 # Test LHS=RHS for u'+ku=kA.
 unassign('t','u0','omega','k'):
 myANS:=your hand-derived formula for u=uh+up:
 LHS:=diff(myANS,t)+k*myANS:
 RHS:=k*(35-13*cos(omega*(t-3))):
 simplify(expand(LHS-RHS));
\end{verbatim}
A successful test of $LHS=RHS$ produces answer {\em zero}, or an
expression that reduces to zero.
%
% unassign('t','u0','omega','k'):
%  myANS:=int((35-13*cos(omega*(x-3)))*k*exp(k*x-k*t),x=0..t);
%  LHS:=diff(myANS,t)+k*myANS:
%  RHS:=k*(35-13*cos(omega*(t-3))):
%  simplify(expand(LHS-RHS));
%
% myANS := -(-35*k^2-35*omega^2+13*k^2*cos(omega*t-3*omega)+13*omega*sin(omega*t-3*omega)*k
% +35*exp(-k*t)*k^2+35*exp(-k*t)*omega^2-13*exp(-k*t)*k^2*cos(3*omega)
% +13*exp(-k*t)*omega*sin(3*omega)*k)/(k^2+omega^2)



\subsection*{Notes on L2.2}

The steady-state solution is derived from the hand-generated
symbolic solution $u=u_h+u_p$ in L2.1 by dropping all terms that
contain $e^{-kt}$. The answer, where $\omega_0=\pi/12$:
$$u_{\mbox{ss}}=35
-\frac{13k}{k^2+\omega_0^2}(k\cos\omega_0
(t-3)+\omega_0\sin\omega_0 (t-3)).$$ To {\em check your answer},
compare it against the formula displayed in L2.2.

To get \MAPLE{} to report the above formula, it is essential to
evaluate everything with $\omega$ as a {\em symbol}, to wit, use
in \MAPLE{} the statement \fbox{\tt unassign('omega'):}. {\bf
Beware} of writing \fbox{\tt omega:=Pi/12}, which defines $\omega$
to be a constant, unless you undo the effect immediately
afterward. Maple's \fbox{\tt unassign} can be used to erase
assignments.

Maple is unable to evaluate \fbox{limit(exp(-k*t),t=infinity)} as
zero, because the expression contains $k$, an unevaluated symbol.
However, it may be possible to simulate the result as in the
following example, which extracts the steady state part
$2\sin(\omega(t-3))$ manually.
\begin{verbatim}
 expr:=exp(-k*t)*cos(omega*(t-3))+exp(-2*k*t)*sin(omega*(t-3))+2*sin(omega*(t-3));
 SS:=subs(exp(-k*t)=0,exp(-2*k*t)=0,expr);
\end{verbatim}

%t:='t':omega:='omega':k:='k':
%AA:=(t,omega)->35-13*cos(omega*(t-3));
%uss:=35-(13*k/(k^2+omega^2))*(k*cos(omega*(t-3))+omega*sin(omega*(t-3))):
%USS:=unapply(uss,(t,omega,k));
%simplify(diff(USS(t,omega,k),t)+k*USS(t,omega,k)-k*AA(t,omega));

\newpage

\subsection*{Notes on L2.3}

The outside temperature  $A(t)$ ({\tt AA(t,omega)} in \MAPLE{})
and the solution $u=U(t,u_0,k,\omega)$ obtained in L2.2 can be
coded in \MAPLE{} as follows:
\begin{verbatim}
 unassign('t','u0','k','omega'):
 AA:=unapply(35-13*cos(omega*(t-3)),(t,omega));
 uss:=35-(13*k/(k^2+omega^2))*(k*cos(omega*(t-3))+omega*sin(omega*(t-3))):
 uss0:=subs(t=0,uss);
 U:=unapply((u0-uss0)*exp(-k*t)+uss,(t,u0,k,omega)):
\end{verbatim}
The two curves are placed onto the \BF{plot} by this
\MAPLE{} command:
\begin{verbatim}
 plot({U(t,74,0.35,Pi/12),AA(t,Pi/12)},t=0..48);
\end{verbatim}
Missing curves? Probably, the missing curves are defined to
contain an unevaluated variable name, like \fbox{\tt pi} instead
of \fbox{\tt Pi}. While {\tt Pi} is the constant 3.14159,
lowercase symbol {\tt pi} is a variable name: {\em case} is
significant in \MAPLE{}.

Curves vanish on the printer? Add \fbox{\tt color=black} to the
plot command, e.g., \fbox{\tt plot(sin(x),x=0..Pi,color=black);}.

Click a mouse button on the high and low spots in the graphic.
Somewhere on the \MAPLE{} worksheet the coordinates of the click
are displayed. This method gives decimal approximations to the max
and min values.

The indoor temperature variation is just the maximum minus the minimum,
as computed from the plot of $u(t)$. Beware: $u$ has a maximum at $t=0$.
Similarly, the outdoor temperature variation is computed from the graph
of $A(t)$.

The phase shift is computed as $|T_2-T_1|$, where $A(T_1)=\max
A(t)$ and $u(T_2)=\max u(t)$. Look at the graphic to find sane
answers for $T_1$ and $T_2$. See the textbook for a more complete
discussion of the ideas. Beware: $T_1$ and $T_2$ are values for
time in the first $72$ hours (abscissa values, not ordinate
values). The smallest phase shift obtainable from the graphic is
to be reported, because the transient solution $Ce^{-ht}$ may
affect the graph of $u(t)=Ce^{-ht}+u_{\mbox{ss}}(t)$ for a few
hours.



\subsection*{Notes on L2.4}

A computer algebra assist for this problem can be found in
\MAPLE{}'s function {\tt implicitplot}. This function can plot the
equation $u(t,74,k,\pi/12)=29$ over the domain $0\le t \le 72$,
$0.2 \le k \le 0.4$. From this plot, and the 3D-plot
$z=u(x,74,y,\pi/12)$, the question is easily answered.
\begin{verbatim}
 with(plots): unassign('t','u0','k','omega'):
 U:=(t,u0,k,omega)->your answer of uh+up from L2.1:
 implicitplot(U(t,74,k,Pi/12)=29,t=0..72,k=0.2..0.4);
 plot3d({U(t,74,k,Pi/12),29},t=0..72,k=0.2..0.4);
\end{verbatim}
Zoom in on the implicit plot by using a smaller time domain,
suggested by the larger plot. The relation between the implicit
plot and the 3D plot is seen by slicing the 3D plot at height
$z=29$ to obtain a {\em bread slice} parallel to the $xy$-plane.
The implicit plot depicts the bread slice by projecting it onto
the $xy$-plane.

Physically, inside temperature $29$ degrees is reached a few hours
{\em after} the outside temperature drops below $29$ degrees.
During $72$ hours, there are three such inside temperature drops,
verified from the 3D-plot, where $z$=temperature. For $k<0.38$,
inside temperature $29$F is not reached in the first 23 hours. For
$0.38 < k < 0.4$, inside temperature $29$ degrees is reached
between $5.5$ and $7$ hours. Please report answers in this
format:
\begin{quote}\sf
No freezing $0$-$23$ hours, $0.2<k < 0.38$, \\
 Freezing $5.5$-$7$ hours for $0.38 < k < 0.4$, \\
 Freezing $28$-$32.7$ hours for $0.2<k<0.4$.
\end{quote}
Use as many ranges as needed to describe the situation; between 5
and 7 will be required.
\end{document}
