% 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}


\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx, caption, subcaption}
\usepackage{url}
\usepackage{hyperref}
\renewcommand{\baselinestretch}{1.3}

\begin{document}

\newpage
\title{\ifnum\TeacherEdition=1 TEACHER VERSION 
\else STUDENT VERSION \fi\\Heat Diffusion}



\markboth{HEAT DIFFUSION}{Heat Diffusion}

\author{Kimberly Spayd and James Puckett\\
Department of Mathematics and Department of Physics\\ 
Gettysburg College\\
Gettysburg PA  17325 USA\\
kspayd@gettysburg.edu, jpuckett@gettysburg.edu
}

\makeStitlePDFLaTex


\begin{abstract}
This project guides students through experimental, analytical, and numerical techniques for understanding the heat (diffusion) equation with nonhomogeneous boundary conditions. In particular, students collect data and model a physical scenario in which heat energy diffuses through a long, thin metal rod that has one of its ends submerged in hot water and the other end in room-temperature water. Students then investigate the subtleties of analytical and numerical techniques used to solve the initial boundary value problem. Finally, the experimental data can be used to compare the theoretical model with the physical phenomenon.
\end{abstract}




\section*{SCENARIO DESCRIPTION}

\subsection*{Background}
The one-dimensional heat (diffusion) equation, 
\begin{equation}
u_t=ku_{xx}, \qquad x\in(0,L), \qquad t>0,
\label{heat}
\end{equation}
is a standard topic considered in introductory Partial Differential Equation courses. In \eqref{heat}, $u=u(x,t)$ is a scalar function that represents the temperature ($^\circ$C) at any time $t$ (s) and position $x$ (mm) measured from the left end of a long thin rod of length $L$. The parameter $k$ represents the thermal diffusivity of the rod material; for steel, $k=4$  mm$^2/$s, while aluminum has a thermal diffusivity of $k=97$ mm$^2/$s, corresponding to its propensity to conduct heat. The possible boundary conditions at $x=0$ and $x=L$ include prescribed temperatures (Dirichlet), insulation (Neumann), and Newton's Law of Cooling (Robin). Besides being an accessible entry point for modeling a physical phenomenon with differential equations and introducing separation of variables as a fundamental analytical technique for solving PDEs, \eqref{heat} provides the basic structure for studying the larger class of parabolic equations. 


Credit for the formulation of \eqref{heat} and its solution belong to Jean-Baptiste Joseph Fourier. At the beginning of the nineteenth century, Fourier proposed analytical solutions which were sums of sine and cosine functions, what we now call Fourier series \cite{APS2010, Freeman1955}. At the time, the veracity of Fourier's heat diffusion model and solution technique were heavily doubted but have since been bolstered by further analysis \cite{APS2010}. History has made him a household name, at least for some mathematicians, because of the pervasiveness and power of his work. Other applications that can be modeled by \eqref{heat} include the spreading of a pollutant in water \cite{Haberman2013}, financial option pricing \cite{Highman2015} and digital image restoration \cite{Schonlieb}. Solution methods based on Fourier series are powerful tools in analyzing particular PDEs on bounded domains. 



\subsection*{Experiment}
The physical phenomenon of heat diffusion in a long thin metal rod is measurable with the experimental design illustrated in Figure \ref{fig:diagram}. The ends of the rod are submerged in water baths at different temperatures and the heat from the hot water bath travels through the metal to the cooler end. The temperature of the rod is measured at four locations; those measurements are sent to a Raspberry Pi, which processes the raw data and sends the collated data to be displayed on the computer screen.
\begin{figure}[h!]
\centering
\includegraphics[width=6in]{heat_flow_diagram.eps}
\caption{ \label{fig:diagram} The experimental apparatus to observe heat diffusion}
\label{fig:diagram}
\end{figure}


If you are collecting your own data, make initial measurements of the length and width of the rod, as well as where the temperature sensors are. Sketch your experimental setup so that you can refer to it later. Before collecting any data, consider the following questions:
\begin{itemize}
\item
What type of metal is your rod made of? What is its expected thermal diffusivity?
\item
How well do you think the metal rod in front of you corresponds to one spatial dimension?
\item
What are the boundary conditions in this set-up? That is, what are the temperatures of the water baths at each end of the rod? Do you think they will remain fixed throughout the demonstration or do you think they will fluctuate? Why?
\item
For short and long periods of time, what do you expect to happen when you submerge the ends of the rod in each water bath?
\end{itemize}
Use the code provided in {\tt PythonCodeRecord.pdf} to start recording data from the temperature sensors. Then submerge the ends of the rod in the water baths. The code provided in \newline
{\tt PythonCodeGraph.pdf} creates a graph of temperature against time for each of the temperature sensors. Collect data until the graph exhibits plateaus at each sensor location for a few hundred seconds. Stop the data recording and save the collected temperatures to use in answering the questions below.


If you are not collecting your own data, see the provided Excel file {\tt HeatEquationData.xls}
in which you can find the time and temperature data that are shown in Figure \ref{fig:data}. To produce this data set, we used an aluminum rod of length $L=300$ mm and square cross-sectional width $3.2$ mm. Foam tubing, with a thickness of $25$ mm, was wrapped around the rod and sensors to provide some insulation. The ambient temperature in the room was $24^\circ$C and the hot water bath measured $53^\circ$C.
\begin{figure}[h!]
\centering
\includegraphics[width=5.5in]{heat_raw_1.eps}
\caption{ \label{fig:data}Temperature readings from four sensors along the rod over the course of the heat diffusion demonstration \cite{Spayd2016}}
\label{fig:data}
\end{figure}


\subsection*{Mathematical Model}
The following assumptions are made in the derivation of \eqref{heat}:
\begin{itemize}
\item
the rod is perfectly insulated so no heat is lost through the surface;
\item
the rod material is uniform;
\item
there is no internal heat generation;
\item
there is no heat applied to the rod on the lateral surface, only on the ends.
\end{itemize}
Reflect on whether these four assumptions are met, approximated, or neglected in the experiment carried out above. How do you think the degree to which these assumptions match the experimental setup will affect the model's accuracy? See \cite{Haberman2013} for more details regarding assumptions and governing equations for the model. 



\subsubsection*{Analysis}
\begin{enumerate}
\item
Write the initial boundary value problem for the heat equation corresponding to the experimental setup observed.
\item
Determine the steady-state solution, $u_s(x)$, of the initial boundary value problem you wrote above.
\item
Explain how and why the analytical technique of separation of variables, i.e. assuming $u(x,t)=X(x)T(t)$, fails with the initial boundary value problem above.
\item
Instead of solving for $u(x,t)$, consider the function $v(x,t)=u(x,t)-u_s(x)$, which measures the displacement of the time-dependent solution from the steady-state solution. Write the associated initial boundary value problem for $v(x,t)$.
\item
Use separation of variables to solve for $v(x,t)$. Explain why the technique is successful with the initial boundary value problem for $v$, in contrast to the initial boundary value problem for $u$.
\item
Now determine the solution $u(x,t)$ of the first initial boundary value problem.
\end{enumerate}


\subsubsection*{Numerics}
\begin{enumerate}
\item
Use Matlab to graph the Fourier series solution for $u(x,t)$ that you obtained in Step (6) of the previous section on Analysis. Discuss the choices necessary to do this.
\item
Read Matlab's help documentation and then use its PDE solver {\tt pdepe} to simulate the solution of the full one-dimensional heat equation problem corresponding to the experimental setup above. Describe how this numerical simulation compares with the plot of your analytical solution, particularly near the initial condition. Discuss a possible cause for the difference between the two plots.
\item
Overlay the data points from the four temperature sensors in the experiment on top of the simulations from the previous question. Discuss how well (or not) the experimental data points match up with the numerical solution. Compare the data and model over time.
\end{enumerate}


\subsubsection*{Errors and Fine-Tuning}
\begin{enumerate}
\item
Describe some possible sources of error between the physical experiment and mathematical model.
\item
Identify any parameters that could be changed to fine-tune the model so that the model better matches the data points. Explain your choices. 
\item
In Matlab, experiment with values for the parameters you identified above to see if there is an improvement in how well the model fits the experimental data. Summarize your findings.
\end{enumerate}



\teacher{
\section*{NOTES FOR TEACHERS}
This is a multi-part group project (approximately four students per group in our experience) that is mainly done outside of class over a couple of weeks. Students are guided through the above questions in a way that reinforces connections between the physical phenomenon of heat diffusion and the mathematical model while considering strengths and weaknesses of the modeling process. Through the project, our students gain experience with modeling using real data as well as using analytical and numerical approaches for understanding solution behavior. In the end, they write a report which summarizes their answers to the questions given in the sections above, providing an opportunity to practice technical writing skills.


Our usual approach is to start this project at the very beginning of the semester, before students have acquired any knowledge of solution techniques for PDEs. We have our students perform the experiment described above during class to become familiar with the physical phenomenon of heat diffusion and develop some intuition for what to expect from the corresponding mathematical model. (We note that one of the authors is a physicist and is the point person for the experimental design and data collection while the other author, a mathematician, leads the modeling portion of this project.) Equipment to build the apparatus seen in Figure \ref{fig:diagram} is listed below. Running the experiment takes approximately twenty to thirty minutes, depending on the type of metal rod used. We recommend that the students run the experiment during the class meeting immediately preceding an introduction to the heat equation and its derivation. Then any physical intuition from the experiment can help orient and guide students through the development of the mathematical model. 


During the following couple of weeks, the in-class focus can be on constructing, analyzing and solving homogeneous initial boundary value problems (IBVPs) with separation of variables, while leaving the bulk of this project for students to complete outside of class. With this group project in hand, students are able to naturally and independently extend the separation of variables technique to inhomogeneous IBVPs. If we find that students face considerable challenges making this extension, we also point them towards the appropriate chapter in the textbook \cite{Haberman2013} for additional support.


In general, students do quite well with answering the questions in the  Analysis and Errors and Fine Tuning sections of this project. The Numerics section can require more attention because of our students' limited experience with programming. We provide the students with initial Matlab scripts to help them get started with the first component of the Numerics section. To further support them with this section, we set aside one class period before the project summary is due to specifically answer Matlab related questions. We meet the students in a computer lab to guide and troubleshoot as necessary, with special attention to the third Numerics question. By the end of the project, students successfully complete all parts and indicate a positive impact on their learning \cite{Spayd2016}.

\subsection*{Experiment}
The experiment uses the following supplies in a single setup:
\begin{itemize}
\item
one Raspberry Pi;
\item
four DS18B20 digital temperature sensors;
\item
one USB power supply;
\item
one USB A to micro B cable;
\item
one 4.7 k ohm resistor;
\item
one metal rod, 1/8 x 1/8 x 24 inches.
\end{itemize}
A supply list with links is available at \cite{Puckett2019}. 
In class, we set up as many as four stations, each with a different type of metal rod, so that each student has an opportunity to interact with the experiment in some way; groups of approximately four students work well to this end. As the students are arriving at the lab and getting settled, we recommend that they sketch the physical setup in front of them and take detailed notes to refer back to later. In particular, they need to measure where the temperature sensors are located on the rod so they can accurately plot their data along with the numerical simulations below. The questions included in the Experiment section are meant to help students initially engage with the project and create some interaction in each group.




The provided Python code in the files mentioned in the Experiment section above process data from the temperature sensors and produce a graph as seen in Figure \ref{fig:data}. This graph updates in real time and students can see the heat diffusion in transient and steady-state regimes. We find this visualization to be extremely helpful for the students when linking the lab experience with analysis of the steady-state solution.



\subsection*{Mathematical Model}
The listed assumptions offered in the Mathematical Model section above can act as discussion points in the lab or classroom after the experiment. Students find it relatively easy to contribute in this way and tend to be more invested in the process of modeling, including recognizing a need to make simplifying assumptions, as a result. The length and depth of any such conversation is entirely left to the discretion of the instructor.


\subsection*{Analysis}
It is assumed that the students are learning about or separation of variables with homogeneous boundary conditions in the classroom. The analytical questions included here are meant to stretch and deepen student understanding of the topic by involving inhomogeneous boundary conditions. In our experience, students are able to navigate through these questions without significant problems. We include solutions here and below for the instructor's convenience.
\begin{enumerate}
\item
With an ambient temperature in the lab ($u=T_L$ $^\circ$C) and temperatures of the water baths $u=T_0$ $^\circ$C, $T_L$ $^\circ$C, we have the initial boundary value problem
\begin{subequations}
\label{NHheat}
\begin{align}
u_t&=ku_{xx}, \quad x\in(0,L), \quad t>0,\\
u(0,t)&=T_0, \ \  \qquad t>0,\\
u(L,t)&=T_L, \ \  \qquad t>0,\\
u(x,0)&=T_A, \qquad x\in(0,L)
\end{align}
\end{subequations}
where time $t$ is measured in seconds and space $x$ is measured in millimeters from the left end of the rod. The value of the thermal diffusivity parameter $k$ (mm$^2/$s) will depend on the specific metal being used. For the instructor's ease, we include thermal diffusivity values for common metals below.
\begin{itemize}
\item
Copper: $k=111$ mm$^2$/s
\item
Aluminum: $k=97$ mm$^2$/s
\item
Steel: $k=4$ mm$^2$/s
\end{itemize}

\item
Letting $u_t=0$ and integrating twice with the boundary values, we have
$u_s(x)=\frac{T_L-T_0}{L}x+T_0$.
\item
Homogeneous boundary conditions for \eqref{heat} allow us to assume that the corresponding spatial ordinary differential equation (ODE) also has homogeneous boundary conditions and that the corresponding temporal ODE is not identically zero. With inhomogeneous boundary conditions, we cannot make any useful assumptions about how the spatial and temporal factors contribute to those values. As a result, we cannot efficiently solve for eigenvalues of the spatial boundary value problem.
\item
The initial boundary value problem for $v$ is (again where, $v(x,t)=u(x,t)-u_s(x)$):
\begin{subequations}
\label{Hheat}
\begin{align}
v_t&=kv_{xx}, \quad x\in(0,L), \quad t>0,\\
v(0,t)&=0, \ \  \qquad t>0,\\
v(L,t)&=0, \ \  \qquad t>0,\\
v(x,0)&=T_A-\left(\frac{T_L-T_0}{L}x+T_0\right), \quad x\in(0,L)
\end{align}
\end{subequations}
again where time $t$ is measured in seconds, space $x$ is measured in millimeters from the left end of the rod, and  $k$ (mm$^2/$s) is the value of the thermal diffusivity parameter for the specific metal used.
\item
Now, with the homogeneous boundary conditions (3b-3c) for $v$, separation of variables can be successfully applied to obtain $\displaystyle v(x,t)=\sum_{n=1}^{\infty}a_n \sin\left(\frac{n\pi x}{L}\right) e^{-k(n\pi/L)^2t}$ where \\
$\displaystyle a_n=\frac{2}{L}\int_0^L \left[T_A- \left(\frac{T_L-T_0}{L}x+T_0\right)\right]\sin\left(\frac{n\pi x}{L}\right)dx$ \cite{Haberman2013}.
\item
Finally, $u(x,t)=v(x,t)+u_s(x)$ with the coefficients given above.
\end{enumerate}


\subsection*{Numerics}

\begin{enumerate}
\item
See Figure \ref{fig:gibbs}. We provide the Matlab file {\tt HeatEqnInClass.m} to our students and discuss it
during class in order to visualize analytical solutions of homogeneous IBVPs; they easily adapt it for this project. 
Students must decide how many terms of the Fourier series to include in the graph. They often stick with $n=500$, which is used in class discussions. While working on this project, students have not yet learned about the Gibbs phenomenon, which is introduced several weeks after the project is due. This particular question in the project is intended to foreshadow the Gibbs phenomenon and we purposefully return to these graphs when the time comes in class.
\begin{figure}[h!]
\centering
\begin{subfigure}{\textwidth}
\centering
\includegraphics[width=3.5in]{gibbs1.eps}
\caption{Truncated analytical solution}
\end{subfigure}
\begin{subfigure}{\textwidth}
\centering
\includegraphics[width=3.5in]{Gibbs_zoom.eps}
\caption{Closer look at oscillation in (a)}
\end{subfigure}
\caption{Graphs used to visualize the truncated analytic solution ($n=500$) of  \eqref{NHheat} \cite{Spayd2016}}
\label{fig:gibbs}
\end{figure}
\item
See Figure \ref{fig:sim}. The oscillations in Figure \ref{fig:gibbs} are no longer present in Figure \ref{fig:sim}. In drawing this contrast, students are able to determine that the truncation of the analytical solution causes the oscillations.
\begin{figure}[h!]
\centering
\includegraphics[width=3.5in]{sim.eps}
\caption{Numerical simulation of  \eqref{NHheat} using Matlab's {\tt pdepe} solver \cite{Spayd2016}}
\label{fig:sim}
\end{figure}

\item
See Figure \ref{fig:snapshots}. The Matlab code provided in the file {\tt HeatEqnSimulation.m} is for the instructor's convenience. In it, the experimental data is imported and then plotted on top of the numerical simulation.
\begin{figure}[h!]
\centering
\begin{subfigure}{4.9 cm}
\centering
\includegraphics[width=4.9 cm]{Heat_eqn_plot_ax.png}
\caption{$t_1=60.3379$ s}
\end{subfigure}
\begin{subfigure}{4.9 cm}
\centering
\includegraphics[width=4.9 cm]{Heat_eqn_plot_1_ax.png}
\caption{$t_2=302.1575$ s}
\end{subfigure}
\begin{subfigure}{4.9 cm}
\centering
\includegraphics[width=4.9 cm]{Heat_eqn_plot_2_ax.png}
\caption{$t_3=799.3076$ s}
\end{subfigure}
\caption{Collected temperature data (dots) and numerical simulation (curve) at times $t_1, t_2, t_3$; these are snapshots of the animation produced for the third question in the Numerics section \cite{Spayd2016}.}
\label{fig:snapshots}
\end{figure}

\end{enumerate}


\subsection*{Errors and Fine-Tuning}

\begin{enumerate}
\item
There are many possible answers; some may include measurement error from the temperature sensors, lack of insulation (or imperfect insulation with the provided data in {\tt HeatEquationData.xls}), and changing temperatures for the water baths. Our students with lab experience tend to identify these easily.
\item
The tunable parameters include the thermal diffusivity, $k$, the boundary conditions $T_0, \ T_L$ and the initial condition $T_A$. While the thermal diffusivity is a given constant, using a different value can be explained by thermal loss due to a lack of insulation in the experiment (or imperfect insulation with the provided data in {\tt HeatEquationData.xls}). The boundary and initial temperatures naturally include error from reading inaccuracies. Additionally, some students may comment that thermal diffusivity may be a space-dependent function while the boundary conditions may be time-dependent functions. This adds complexity to the model that may be able to account for some of the discrepancy between the data and model.
\item
Students can rerun the Matlab code with different values of the parameters to see whether they improve the fit. In our opinion, it is beyond the scope of this project for students to analyze or numerically simulate the heat equation with non-constant parameters. Our students usually give adequate qualitative summaries based on visual interpretations; so far, no group has provided a statistically rigorous analysis to validate the model.


\end{enumerate}


\begin{thebibliography}{99}

\bibitem{APS2010} Chodos, A. Editor. 2010. This Month in Physics History March 21, 1768: Birth of Jean-Baptiste Joseph Fourier. {\em APS News}, 19(3): 2.

\bibitem{Freeman1955} Freeman, A. Translator. 1955. {\em The Analytical Theory of Heat} By Joseph Fourier. New York: Dover.

\bibitem{Haberman2013} Haberman, R. 2013. {\em Applied Partial Differential Equations with Fourier Series and Boundary Value Problems, Fifth Edition}. Boston  MA: Pearson.

\bibitem{Highman2015} Highman, N. J. Editor, 2015 {\em The Princeton Companion to Applied Mathematics}. Princeton  NJ: Princeton.

\bibitem{Puckett2019} 
Puckett, J. 2019. \url{http://public.gettysburg.edu/~jpuckett/heatEqn.html}. Accessed 14 October 2019.

\bibitem{Schonlieb} Sch{\"o}nlieb, C.-B. 2019.  \href{https://www.mathworks.com/company/newsletters/articles/applying-modern-pde-techniques-to-digital-image-restoration.html}{https://www.mathworks.com/company/newsletters/articles/applying-modern-pde-techniques-to-digital-image-restoration.html}. Accessed 14 October 2019.

\bibitem{Spayd2016} 
Spayd, K. and J.~Puckett. 2016. A Three-Fold Approach to the Heat Equation: Data, Modeling, Numerics. {\em PRIMUS}. 26(10): 938--951.

\end{thebibliography}

}





\end{document}
