% SIMIODESample.tex and simiode.cls for LaTeX2e
% Version 11 March 2013.
% 
% This class (simiode.cls) is a modification of that used by 
% the journal PRIMUS - Problems, Resources, and Issues in 
% Mathematics Undergraduate Studies. 
% 
% See http://www.tandfonline.com/loi/upri .
% 
% The class was developed by  Marc Frantz with the  suggestions of the 
% PRIMUS Editorial Board, especially John Thoo, and is much appreciated.
%
% This template "SIMIODEPaper.tex" is to be used with the 
% class file "simiode.cls".  You may rename the template.  
% The latest version of this class and template files 
% can be found at
% 
%                <http://www.SIMIODE.org> .
%
% For an explanation of using the class file "simiode.cls", 
% see the pdf ``Description and Illustration of Use: 
% A LaTeX Class File for SIMIODE'' found at  
% 
%                <http://www.SIMIODE.org> . 
%
% A "Quick Start Guide" appears below.
%
% For information on using TeX or LaTeX, go to %                 <http://www.tug.org>.

% This file must be processed with version LaTeX2e or higher.
%
\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}
%
% The \teacher{ } command wraps any text that will appear only in the teacher version of
% the document. See this file for examples. This has been done for you in the document below.   
% Rename the pdf after you compile each final version.
%

\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx}
\usepackage{url}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NOTE regarding double spacing %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Manuscripts should be submitted doubled spaced.  This is 
% accomplished via the \baselinestretch command below.
% However, you may obtain an idea of how your article will 
% look in final form, so that you can better plan your 
% presentation.  To see your manuscript formatted very 
% similarly to the on-line published format of  SIMIODE, 
% comment the following line out by inserting a percent sign % (%) at its beginning.

\renewcommand{\baselinestretch}{1.3}

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NOTE regarding macros %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Any macros defined for your paper should be contained in 
% the top matter. Likewise, any environment definitions such 
% as \newtheorem or \newenvironment should also go in the 
% top matter. 
% 
% Do not \input your macros as separate files within your 
% final version because more files make it harder for the 
% editors and users to keep track of and/or modify and 
% append your paper.


\newcommand{\AmSLaTeX}{$\cal A$\kern-.1667em\lower.5ex\hbox{$\cal
M$}\kern-.125em $\cal S$-\LaTeX}
\newcommand{\BibTeX}{{\sc Bib}\TeX}
\newcommand{\scls}{{\tt simiode.cls}}
\newcommand{\stemp}{{\tt SIMIODEPaper.tex}}
\newcommand{\SI}{SIMIODE}

\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{proposition}[theorem]{Proposition}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Top matter is above and 
% the beginning of the document below.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

\newpage

\title{\ifnum\TeacherEdition=1 TEACHER VERSION \else STUDENT VERSION \fi\\MACHINE REPLACEMENT - LAPLACE TRANSFORM}


\markboth{Modeling Machine Replacement}{Modeling Machine Replacement}






\author{Brian Winkel,
Director SIMIODE\\
Chardon OH 44024 USA\\
BrianWinkel@simiode.org}


\makeStitlePDFLaTex



\begin{abstract}
Students build an integro-differential equation model using a convolution for  machine replacement strategies for two different machine failure models: (1) exponentially distributed failure (student exercise) and (2) fixed time replacement. We discuss all the necessary probability notions which support the building of the convolution model and solve it for (2)  while asking students to solve a general situation and then apply their findings to  address (1).  

We apply the notion of convolution and Laplace transform to a machine maintenance problem in a large facility. We seek the replacement function $R(t)$ for the number of machines to be replaced in the time interval $[0,t]$ to meet a set demand function, and $N(t)$, for the number of machines in operation at time $t$. We are given probability information  on $F(t)$, the fraction  of the machines that were new at time $t_0$ which are still  in operation at time $t_0 + t$, for the two different  scenarios (1) and (2) above. 

We offer a separate Appendix with some introductory probability notions.
\end{abstract}



\section*{SCENARIO DESCRIPTION}


The use of the Laplace transform in solving differential equations permits one to change the nature of the solution strategy from solving a differential or derivative equation in time ($t$) to solving an algebraic equation in the so-called ``frequency'' domain ($s$). Then return to the time domain for a solution of the original differential equation.  This process has three steps.  
\begin{enumerate}
  \item The starter step is transforming the differential equation in some time dependent function, $y(t)$, using  the Laplace transforms of the  parts of the differential equation to obtain an algebraic equation in, $Y(s)$, the Laplace transform of $y(t)$. This step has a prescribed approach using a set of developed rules, either found in tables in texts or implemented in a computer algebra system. 
  \item The middle step is isolating or solving algebraically for the transform, $Y(s)$, of the desired solution function, $y(t)$, from the original differential equation. 
  \item The final step is going ``back" from the frequency domain, $s$, function $Y(s)$, to the time domain, $t$, function $y(t)$.  This is often the hard part. Indeed, it is possible that we cannot complete this inverse process and might have to resort to a numerical inverse Laplace transform. This latter possibility is beyond most differential equations courses and we will not consider that here.
\end{enumerate}

Various functions such as  the step or Heaviside function, the Dirac function, or the convolution of functions lend themselves to Laplace transform analysis and are often found in many engienering applications. However, in the field of industrial engineering, where  we are often given probabilistic information, e.g.,  how long a part will last, the  notion of the convolution occurs quite naturally. Here we  see it arise in a machine replacement application in which we have some probability distribution of machine failure and we seek to determine a replacement  procedure so that we keep a specified number of machines in operation in our facility.



\section*{Machine failure and replacement}


We consider the situation where we have a number of identical machines in a plant, e.g., stamping machines in a metal shop, printing machines at a publisher, or weaving machines at a textile factory, etc. 
The following modeling opportunity, with guidance added, originates from \cite[pp. 261-262]{FinneyOstberg1968} and we quote and paraphrase some of the original problem here. This is an excellent  application of differential equations and Laplace transforms to industrial  engineering.

\vspace*{.2 true cm} 
A  factory requires $N(t)$ machines to be in operation at time $t$. It is known that in any time interval $[t_0,t_0 + t]$ only a fraction, $F(t)$,  of the machines that were operating at time $t_0$ will still be in operation at time $t_0 + t$. Moreover, $F(t)$ depends only upon $t$ and not  $t_0$, i.e. survival of a given machine depends upon the length of time the machine is in service, not when the machine actually starts service, $t_0$. We wish to find a replacement function $R(t)$, where $R(t)$ is the number of machine replacements needed in the time interval $[0,t]$ to assure that $N(t)$ machines are in operation at time $t$. We will assume that the number of machines is so large that $N$, $F$, and $R$ may as well be continuous.

Here are activities to help you realize the power of convolution and the action of Laplace transforms in industrial engineering applications.

\begin{itemize}
\item [1)] Under the assumption that $R$ is differentiable show that the number of
replacements required in the time interval $[v_k,v_k + \Delta v_k]$ is approximately
$R'(v_k)\Delta v_k$,
where $\Delta v_k = v_{k + 1} - v_k$, and $0 < v_1 < \ldots < v_n = t$. Hint:  What does each term in the  expression, $\frac{R(v_k+ \Delta v_k) - R(v_k)}{\Delta v_k} \cdot \Delta v_k$, represent?
 
\item [2)] Use the result in activity (1) to show that


\begin{equation}\label{eq:sum}
N(t) \approx N(0)F(t) +\sum_{k=1}^n F(t - v_k)R'(v_k)\Delta v_k\, ,
\end{equation}

\noindent and hence, in the limit


\begin{equation}\label{eq:int}
N(t) = N(0)F(t) + \int_0^t F(t - v)R'(v)\, dv\, .
\end{equation}

\noindent Hint: What does $F(t - v_k)(R(v_k+ \Delta v_k) - R(v_k))$ represent and what does the sum of all such terms over each respective time interval $[0,v_1]$, $[v_1,v_2]$, \ldots $[v_{n-1}, v_n]$ mean?  Remember our discussion about  independent and mutually exclusive  events.
\end{itemize}






\subsection*{Laplace Transform and Convolution}

We first define  the {\it Laplace transform\/} of a function $f(t)$:

\begin{equation}\label{eq:LT}
 {\cal L}(f)(s) = \int_0^{\infty} f(t) e^{-s\, t} \, dt\, .
\end{equation}

Notice the product $s\, t$ in the exponential term. Now  due to the difficulties in handling units in exponents we regard the unit for $s$ as being $1/t$. Hence the Laplace transform is thought of as acting on the frequency domain, $s$, rather than the time domain, $t$.

If we have two functions $f$ and $g$ we define the {\it  convolution\/} of $f$ and $g$,  written as $f*g$, using an asterisk or star, to be the integral of the product of the two functions with one  domain reversed and shifted, i.e.

\begin{equation}\label{eq:convo}
(f*g)(t) = \int_{-\infty}^{\infty} f(\tau) g(t - \tau)\, d\tau \, . 
\end{equation}

There are lots of applications of the convolution, in particular in electrical engineering, but in our context we are using it in an industrial engineering setting to solve a problem of maintaining the correct number of machines in a plant.

There are many elegant properties of the convolution, e.g., it is a commutative operation: 
 $$(f*g)(t) = \int_{-\infty}^{\infty} f(\tau) g(t - \tau)\, d\tau \, = \int_{-\infty}^{\infty} g(\tau) f(t - \tau)\, d\tau = (g*f)(t)\, . $$ 

A very useful property of the convolution with respect to the Laplace transform of a function is the following:

\begin{equation}\label{eq:Conv}
{\cal L}(f*g)(s) = {\cal L}(f)(s)\cdot {\cal L}(g)(s)\, ,
\end{equation}

We do not offer a proof of (\ref{eq:Conv}) here but suggest to the reader that this can be found in a number of places (text and web) under the formal name of Convolution Theorem and involves changing of variables and switching the order of integration. 

Now in the context of our problem these notions of convolution and Laplace transform  of the integral in (\ref{eq:int}) yield the following

\begin{equation}\label{eq:transum}   
{\cal L}\left(\int_0^t F(t-v)R'(v)\, dv\right)(s) = {\cal L}(F)(s)\cdot {\cal L}(R')(s) = {\cal F}(s)\cdot (s {\cal R}(s) - R(0))\, .
\end{equation}
\noindent In light of   (\ref{eq:Conv})  it looks like we are really ready for our next issue if we can do inverse Laplace transforms. Let's go for it!

\begin{itemize}

\item [3)] Use (\ref{eq:int}) to find a formal expression for the Laplace transform of $R(t)$ given that $R(0) = 0$. Why is $R(0) = 0$? We will use this result extensively with a variety of machine requirements, $N(t)$, and machine failure rates or distributions, $F(t)$.

\item [4)] Find $R(t)$ when $N(t) = N_0$, a constant, and $F(t) = e^{-a t}$, $a > 0$. The former statement simply means that machine  needs are constant. This is quite often the case, as there is usually a steady demand and a steady supply of workers to run the machines. 

$F(t) = e^{-a t}$, $a > 0$,  basically says the operating life of a machine starting up at time $t = 0$ and running over the entire period $[0,t]$ exponentially decays, i.e. the fraction of the machines that were new at time $t=0$   and are still running at time $t$ decays and is given by $F(t)$ or the probability of a specific machine NOT failing in the time interval $[0,t]$ is  $F(t)$. In this case, $G(t) = 1 - F(t) = 1 -  e^{-a t}$ is called the {\it exponential distribution\/}, where $G(t)$ is the probability that a machine which has started at time $t=0$ WILL FAIL  in the time interval $[0,t]$, i.e. is NOT running at time $t$ or put another way, the machine fails in the time interval $[0,t]$. We say that $G(t)$ is the {\it distribution function of the random variable, $X$\/}, where $X$ is the first time of failure for a given machine, so ${\rm Prob}[X < t] = G(t) = 1 - F(t) = 1 - e^{-a t}$.  

Now $G'(t)= a e^{-a t}$, loosely speaking, is the probability that the time of failure is $t$. The integral of $G'(t)$ over an interval is the probability that the time of failure occurs somewhere in that interval.  $G'(t)$ is called the {\it probability density function\/}.  If we add up the probability of the mutually exclusive  outcomes $t$ for machine failure times the time or value, $t$, of each outcome we get the {\it expected value\/} or {\it average\/} value of the random variable, $X$, i.e., the average time of failure. 



In the continuous distribution case this is an  integral, using the {\it probability density function}, $G'(t) = a e^{-a t}$:
\begin{equation}\label{eq:EV}
\int_0^{\infty} t G'(t) \, dt = \int_0^{\infty} t (a e^{-a t}) \, dt = \frac{1}{a}\, .
\end{equation}
 

\item [5)] Offer up some plots of $R(t)$, as found in activity (4),  for the various parameters in this situation, e.g., $N_0 = 100$, $a = 0.1$ Then alter $a$ and explain the significance of these parameters, again, and of your plots. Interpret  what is happening in this situation using your plots.

\item [6)] Suppose we vary the demand for machines, i.e., $N(t)$ is something other than constant. Say,  we keep increasing the demand over time, $N(t) = N_0 + b t$.  Using exponential failure distribution, $F(t) = e^{-a t}$, which we have used before, solve for $R(t)$, the machine replacement rate and offer some plots using values $N_0 = 100$ with $a = 0.01$ and $b = 0.1$ and then $b = 0.2$. Compare the plots and explain the difference. Tell why these differences are reasonable.
\item [7)] What if  we  add some trend activity ($c\,t$) and also add seasonality ($\sin\left( \frac{2 \pi}{12} t\right)$) to  the demand for machines, i.e., $N(t)$ is something other than constant, e.g.  we keep increasing the demand over time, $N(t) = N_0\left(1 + b \sin\left( \frac{2 \pi}{12} t\right) \right)
 + c t$?  Use exponential failure distribution, $F(t) = e^{-a t}$, which we have used before, and solve for $R(t)$, the machine replacement rate.  Then   offer some plots using values $N_0 = 100$ with $a = 0.01$ and $b = 0.05$ (seasonal variation of 5\% in machine demand) and $c = 0.2$ (steady increase of 0.2 machines per unit time). Compare the plots here with those of activity (6) and explain the difference. Tell why these differences are reasonable.
\end{itemize}

\subsection*{Another possible replacement scheme}

\begin{itemize}
\item [8)] Now instead of waiting for machines to fail naturally or on their own, i.e., using exponential distribution for modeling machine failure ($F(t) = e^{-a t}$, $a >0$), suppose we just plan to replace the machines - all of them at once - at a fixed time.   Find a formal expression for $R(t)$  when $N(t)$ is arbitrary and every machine has the same finite life span $L$.  This is the case of  sweeping through and replacing all machines at once, as is often done with light bulbs in offices or plants.

This can be modeled using a Heaviside or unit step function (\ref{eq:UnitStep}) to model the machine failure rate, i.e. 

\begin{equation}\label{eq:UnitStep}
u_L(t) =  \left\{  
\begin{array}{ll}
0 & \mbox{if $t< L$}\\
1 &  \mbox{if $t\geq L$}
\end{array}
\right. \, .
\end{equation}

We use a machine failure distribution  of:

\begin{equation}\label{eq:StepFailure}
F(t) = 1 - u_L(t) =   \left\{  
\begin{array}{ll}
1 & \mbox{if $t< L$}\\
0 &  \mbox{if $t\geq L$}
\end{array}
\right. \, ,
\end{equation}
which says that the probability of a given machine lasting until time $t = L$ is 1 or certain, while the probability of a given machine lasting longer then time $L$ is 0. 
\end{itemize}  





One problem with this approach is that there may be a machine which fails before we do our total replacement. However, usually this replacement is way before any given machine is likely to falter. Usually this replacement is done  so as not to inconvenience workers in a shift at a time outside the working shift,  say at night or over the weekend.

From activity (3) we recall solving our general integro-differential equation (\ref{eq:int}) and finding (you did, did you not?) the transform, ${\cal R}(s)$,  for $R(t)$ in terms of the other transforms:  


\begin{equation}\label{eq:TransRSol}
{\cal R}(s) = \frac{{\cal N}(s)- N(0) {\cal F}(s)}{s {\cal F}(s)} \, .
\end{equation}

\noindent If in (\ref{eq:TransRSol}) we use the machine failure distribution, (\ref{eq:StepFailure}), then we can find ${\cal R}(s)$ in this case:  

\begin{equation}\label{eq:FinRSol}
{\cal R}(s) = \frac{N(0)}{s \left(-1+e^{L s}\right)} \, .
\end{equation}

This is a sticky frequency domain transform to invert, so we need to make a change. We first notice that the expression $\frac{1}{1-e^{L s}} $ can be expanded into an infinite geometric series with ratio $e^{L s}$, but this will not converge in general as $L>0$ so we multiple top and bottom of (\ref{eq:FinRSol}) by $e^{-L s}$ to obtain (\ref{eq:FinalRSol}):

\begin{equation}\label{eq:FinalRSol}
{\cal R}(s) = \frac{e^{- L s}N(0)}{s \left(1-e^{-L s}\right)} \, .
\end{equation}


\begin{equation}\label{eq:expand}
\frac{1}{1-e^{-L s}} = \left(  1+e^{-L s}+e^{-2 L s}+e^{3 L
   s}+e^{-4 L s}+e^{-5 L s}+e^{-6 L
   s} \ldots \right) = \sum_{i=0}^{\infty} e^{-i L s} \, .
\end{equation}

This means that we have a usable form (in a series, albeit) for the transform, ${\cal R}(s)$ of our replacement function, $R(t)$.

\begin{equation}\label{eq:FinRRSol}
{\cal R}(s) =  \frac{N(0) e^{- L s}}{s \left(1-e^{-L s}\right)} = N(0)\left( \sum_{i=1}^{\infty} \frac{e^{-i L s}}{s}\right) \, .
\end{equation}

Now, in (\ref{eq:FinRRSol}) each term, $\displaystyle{\frac{e^{-i L s}}{s}}$, $i = 1$, 2, 3, \ldots , has an inverse Laplace Transform, namely 
\begin{equation}\label{eq:piece}
{\cal L}^{-1}\left(\frac{e^{-i L s}}{s}\right)(t) = u_{i L}(t) = {\rm UnitStep}(t-i L) \, .
\end{equation}

\begin{itemize}
\item [9)] Verify (\ref{eq:TransRSol}) and (\ref{eq:FinalRSol}) for ${\cal R}(s)$ and then for activity (8) show that
 \begin{equation}\label{eq:Nsol}
R(t) = N(0) \sum_{i=1}^\infty u_{iL}(t) = N(0) \sum_{i=1}^\infty {\rm UnitStep} (t - iL) \, .
\end{equation}


\noindent Another way of writing (\ref{eq:Nsol}) is

\begin{equation}\label{eq:NsolAlt}
R(t) = \left\{   
\begin{array}{ll}
N_0 - N_0 = 0  & \mbox{if $0 \leq t < L$}\\
2N_0 - N_0 = N_0 & \mbox{if $L \leq t < 2L$}\\
3N_0 - N_0  = 2 N_0& \mbox{if $2L \leq t < 3L$}\\
 \quad \vdots\\
\end{array}
\right.
\end{equation}

\noindent Now consider the cases in activities (9) and (10) where $L = 2$ and $N_0 = 100$. For each activity, offer plots of $R(t)$ for $0 \leq t < 40$. Explain what these plots signify.
  \end{itemize}


\section*{APPENDIX -- Essential Probability Details}




\subsubsection*{Independent events}

Consider two events: A in which we toss a fair two-sided coin with one side Head and one side Tail and B in which we toss a fair six-sided die, with sides numbered 1, 2, 3, 4, 5, and 6, and note the value on the side facing up. These are called {\it independent events\/}, meaning the outcome of one does not effect the outcome of the other, so the probability of a ``Head'' on a coin toss followed by a ``3''  on a toss of the die is just the product of the two probabilities, namely probability of a ``Head'' on a coin toss times the  probability of  a ``3'' on a die toss.  Think about this when considering each term in the sum found in the right hand side of (\ref{eq:sum}) below and use this terminology in  addressing the questions below.


\subsubsection*{Mutually exclusive events}



The two outcomes of ``3 or 4'' and ``5 or 6'' in event B, i.e. die toss, are said to be {\it mutually exclusive events} as there is no overlap between them. Thus, to determine the probability we get a ``3 or 4'' OR  ``5 or 6'' we add the two probabilities of the two mutually exclusive events, $A$ and $B$, i.e., (probability we get a ``3 or 4'')  $+$  (probability we get ``5 or 6''). Think about this when considering the addition of terms  in the sum found in the right hand side of (\ref{eq:sum}) below and use this terminology in  addressing the questions that follow.


Consider an experiment  in which we toss  two fair six-sided die, each with sides numbered 1, 2, 3, 4, 5, and 6, and add the values on the sides facing up. Consider two events $A$ and $B$ where $A$ is the event where the sum of the two faces is 6 and $B$ is the event where the sum of the two faces is 8.   These events are examples of two events which are {\it mutually exclusive events\/}, i.e. their intersection is empty or they have nothing in common. Thus to determine the  probability that either $A$ or $B$ occur we simply add the probability that $A$ occurs to the probability that $B$ occurs. We note the difference between this situation with one involving two events:  $A$ again and $C$ - the sum is a non-zero multiple of 3. These are not mutually exclusive events,  so in order to determine the probability that either $A$ or $B$ occur we have to subtract off the probability that both occur, i.e., ${\rm Prob}(A \cup B) = {\rm Prob}(A) +  {\rm Prob}(B) - {\rm Prob}(A \cap B)$.  We will be dealing with mutually exclusive events in our situation, i.e., $A$ and $B$ such that $A \cap B = \phi$, the empty set.  


\section*{NOTES FOR TEACHERS}
   
We stumbled on this problem years ago in this classic text \cite[pp. 261-262]{FinneyOstberg1968}.  We believe this is a nice application which combines elementary probability concepts, Laplace transforms, and convolution to solve a real industrial engineering problem. We have not seen this application elsewhere in the differential equation text literature, but we have found that students, while challenged by the introduction of probabilistic notions in their differential equations course(!),  enjoy seeing this application.

In the STATEMENT for students we have tried to weave new notions and introduce students to the essential probability ideas so that the situation has a vocabulary and  a set of tools. We attempt to get students through the analysis and also have them reflect on the applicability of the mathematics and the meaning in each situation.

Let us look at what a student would have to do to consider her/himself successful in this scenario. We take it activity by activity.

\begin{itemize}
\item [1)] Since $R(t)$ is the number of machine replacements  needed in the time interval, $[0,t]$, then $R(t+a) - R(t)$ would be the number of machine replacements needed in the time interval $[t,t+a]$. Thus $R(\tau_k+ \Delta \tau_k) - R(\tau_k)$ would be the number of machine replacements needed in the time interval $[\tau_k,\tau_k+\Delta \tau_k]$. Replacing $R(\tau_k+ \Delta \tau_k) - R(\tau_k)$ by the fraction term (essentially multiplying top and bottom by $\Delta \tau_k$) yields  $\frac{R(\tau_k+ \Delta \tau_k) - R(\tau_k)}{\Delta \tau_k} \cdot \Delta \tau_k$,      which as $\Delta \tau_k \to 0$, as it must to give us a continuous model, gives us $R'(\tau_k) \Delta \tau_k$.
\item [2)] One should ask, ``How would we get $N(t)$ machines in operation at time $t$?'' We could break up the time interval $[0, t]$  to produce mutually exclusive events over the subintervals $[\tau_k, \tau_k + \Delta \tau_k]$ and just add them up. What do each of the terms look like?  The big concept with regard to mutually exclusive events here is  the fact that  machines are replaced in the time interval $[\tau_k, \tau_k + \Delta \tau_k]$ and  continue to operate until the end of the entire time interval, i.e. until time $t$. Now we just have to count the number of machines replaced in each of the small intervals, $[\tau_k, \tau_k + \Delta \tau_k]$ which continue to work until time $t$.   Each machine in operation at time $t$ got to that position in one and only one of these ways. In addition there are the machines that were running at the start of the time interval $[0,t]$, there were $N(0)$ such machines, and did not need any repair, i.e. the fraction of theses machines which kept on running over the entire time interval, $F(t)$.

Now in each interval, say $[\tau_k, \tau_k + \Delta \tau_k]$, if a machine were replaced  we are interested in its running from time $\tau_k$ to time $t$. Here we use $\tau_k$ as the representative time in the time interval $[\tau_k, \tau_k + \Delta \tau_k]$ for replacement and since we are going to be taking the limit as $\Delta \tau_k \to 0$ (as we do in all integration or accumulation processes) this will be adequate. But how do we represent the fact that the machine must then run from $\tau_k$ to $t$? We use, $F(t - \tau_k)$, the fraction of the machines  that were new $t-\tau_k$ time units ago and  are still running after these $t-\tau_k$ units of time.  In this context one sees the notion of convolution arising naturally.

Finally, we can add up all the ways a machine operating at time $t$ ($N(t)$ of them) got there:

\begin{equation}\label{eq:sumagain}
N(t) \approx N(0)F(t) +\sum_{k=1}^n \underbrace{F(t - \tau_k)}_{{\rm runs \,\,\,from\,\,}\tau_k {\rm \,\, to\,\,} t}\cdot \overbrace{R'(\tau_k)\Delta \tau_k}^{{\rm replaced \,\,\, in\,\,} [\tau_k,\tau_k+\Delta \tau_k]}\, ,
\end{equation}

\noindent and in the limit as $\Delta \tau_k \to 0$ we obtain the integral equation found in (\ref{eq:int}).
\item [3)] This is where the Laplace transform  becomes very helpful as we  transform both sides of (\ref{eq:int}) and solve for ${\cal R}(s)$ the Laplace transform of $R(t)$, from whence we can use inverse Laplace methods to finally obtain $R(t)$.

We take the Laplace transform of both sides of (\ref{eq:int}), while noting the presence of the convolution in the integral $\int_0^t F(t - v)R'(v)\, dv$ and the fact that $R(0) = 0$. This latter fact holds true because at the start ($t=0$) we have all machines running and will need no new replacement at time $t=0$. Thus we have the transform equation (\ref{eq:Laplace}) in the frequency domain, $s$:

\begin{equation}\label{eq:Laplace}
{\cal N}(s) = N(0) {\cal F}(s) + {\cal F}(s)\cdot (s {\cal R}(s) - R(0))\, , \qquad R(0) = 0
\end{equation}

\noindent Solving for ${\cal R}(s)$ we obtain
\begin{equation}\label{eq:Trans}
{\cal R}(s) = \frac{{\cal N}(s)- N(0) {\cal F}(s)}{s {\cal F}(s)} \, .
\end{equation}
\item [4)]When $N(t) = N_0$, i.e. a most realistic situation occurs in which there is a need for a fixed number of machines at all times, we   merely have to apply (\ref{eq:Trans}) and perform the inverse Laplace transform on ${\cal R}(s)$ to obtain $R(t)$. In this case we use $F(t) = e^{-a t}$, $a > 0$. Thus, we can go back to our original formulation in (\ref{eq:int}), substituting the appropriate pieces to include $F(t) = e^{-a t}$, and obtain for ${\cal R}(s)$:

\begin{equation}\label{eq:p4}
{\cal R}(s) = \frac{a N_0}{s^2} \, .
\end{equation}
Now using inverse Laplace transforms we find that $R(t) = a N_0 t$. What does this really mean? It means that as time goes on (since the machines are failing more frequently due to the distribution $F(t) = e^{-a t}$ there will be a need for more replacements. A bigger value for $a$ means that as machines have an average survival time of $\frac{1}{a}$ and, as such, need to be replaced more often as time goes on, we see that here  $R(t)$ is proportional to $a$. Recall $R(t)$ is the number of machines needed in the time interval $[0,t]$ so we expect it to increase as $t$ increases also and $R(t) = a N_0 t$ IS increasing, linearly with a coefficient of $a$ and the required constant number of machines in operation, $N_0$.

\item [5-9)] We have placed extensive analysis into a {\it Mathematica\/} notebook, 7-8-T-Mma-MachineReplace- mentLaplaceTransform-TeacherVersion.nb 
\end{itemize}


\begin{thebibliography}{1}



\bibitem{FinneyOstberg1968} Finney, R. L. and D. E. Ostberg. 1968. {\it Elementary Differential Equations With Linear Algebra\/}. Reading MA:   Addison-Wesley.

\end{thebibliography}



\end{document}

 
\item In this case we need to find  a formal expression for $R(t)$  when $N(t)$ is arbitrary and every machine has the same finite life span $L$. To represent the latter we offer the following function for $F(t)$:

\begin{equation}\label{eq:fStep}
F(t) = 1 - u_L(t), {\rm \,\,where\,\,} u_L(t) = u(t-L) = \left\{
\begin{array}{c} 
0, {\rm \,\,if}\,\,t < L\\
1, {\rm \,\,if\,\,} t \geq L
\end{array}\right .
\end{equation}
Here $u_L(t) = u(t-L)$ is a step function.
\end{enumerate}
\newpage 
