Warming, Cooling, and Urban Ozone Pollution

Part 1. Background: Newton's Law of Cooling

Explain why the solution of the initial value problem

	(dT/dt) = -k (T - a), T(0) = (Subscript[T, 0])
is
	T = a + ((Subscript[T, 0]) - a) (e^(-k t)).

Answer:

Explain why

	((T(i + 1) - T(i - 1))/(t(i + 1) - t(i - 1)))

is a good estimate of (dT/dt) at t = (Subscript[t, i]).

Answer:

Part 2. The Cooling Curve

Enter the following to load the necessary package.

<< Graphics`Graphics`

Enter the following lines with time and temperature data.  In order to plot the data, we "thread" times and temps together into a list of pairs.  The first command, incidentally, suppresses the warning that Mathematica would have generated about the similarity between our variable times and the built-in function Times.

Off[General::spell1];

times  =  {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 
          130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250, 
          260, 270, 280, 290, 300, 310, 320, 330, 340, 
          350};
temps  =  {95.5, 92.4, 90.3, 88.1, 86.7, 85.3, 84.2, 83.2, 81.9, 
          80.8, 80.0, 79.2, 78.6, 78.0, 77.4, 76.9, 76.2, 75.8, 75.4, 75.1, 
          74.7, 74.4, 74.2, 73.9, 73.7, 73.6, 73.4, 73.3, 73.1, 73.1, 73.0, 
          72.9, 72.8, 72.7, 72.7, 72.6};
tempData  =  Thread[{times, temps}]

Plot the cooling data.

tempPlot  =  
        ListPlot[tempData,  PlotRange -> {70, 100},  
          PlotStyle -> {RGBColor[0, 0, 1], PointSize[0.015]},  
          AxesLabel -> {"Time",  "Temp"}];

To change the vertical scale, we modify the temps list, and thread times and scaledTemp for plotting.

scaledTemp  =  temps  -  72.3;
scaledData  =  
        Thread[{times,  scaledTemp}];
ListPlot[scaledData,  PlotRange -> {0, 25},  
        PlotStyle -> {RGBColor[0, 0, 1], PointSize[0.015]},  
        AxesLabel -> {"Time",  "Temp"}];

Create a semilog plot of the scaled temperature data.

LogListPlot[scaledData, 
        PlotStyle -> {RGBColor[0, 0, 1], 
            PointSize[0.015]}];

Answer to the question in step 3:

Find the decay constant k.

k  =  ?;

Enter the model temperature function.

T[t_] :=  ?;

Plot the data and model function together.

expPlot  =  
        Plot[T[t], {t,  0, 360},  
          DisplayFunction -> Identity];
Show[expPlot, tempPlot,  DisplayFunction -> $DisplayFunction, 
        AxesLabel -> {"t", "Temp"}];

The following procedure computes symmetric difference quotients from any list of data points.  Note that the resulting list is two elements shorter, because there is not enough information for a symmetric quotient at the beginning and end of the original list.  Enter the definition, but don't worry about the details.

sdq[data_] :=  
        Table[(data[[k + 1]][[2]] - data[[k - 
                        1]][[2]])/(data[[k + 
                        1]][[1]] - data[[k - 
                        1]][[1]]),  {k,  2,  
            Length[data] - 1}];

Compute slopes from the temperature data, and plot the slopes against temperature.  Note that temp must be appropriately shortened beforing threading with slopes.

slopes  =  sdq[tempData];
slopeData  =  
        Thread[{Delete[temps, {{1},  {Length[temps]}}], 
            slopes}];
ListPlot[slopeData, 
        PlotStyle -> {RGBColor[1, 0, 0], 
            PointSize[0.015]}];

Answer to the question in step 6:

Part 3. The Warming Curve

Enter the data for the warming phase of the experiment.

times  =  {0, 5, 10, 20, 30, 40, 50, 60, 
          70};
temps  =  {72.3, 87.2, 90.5, 93.1, 94.2, 94.7, 95.1, 95.3, 
          95.5};
tempData  =  Thread[{times, temps}]

Plot the data.

tempPlot  =  
        ListPlot[tempData, PlotRange -> {70, 100},  
          PlotStyle -> {RGBColor[1, 0, 0], PointSize[0.015]}, 
          AxesLabel -> {"Time", "Temp"}];

Use your own estimate of ambient temperature to scale the warming data.

ambient  =  ?;
scaledTemp  =  ambient  -  temps;
scaledData  =  Thread[{times, scaledTemp}]

Plot the scaled data.

ListPlot[scaledData,  PlotRange -> {0, 25},  
        PlotStyle -> {RGBColor[1, 0, 0], PointSize[0.015]},  
        AxesLabel -> {"Time",  "Temp"}];

Make a semilog plot of the scaled data.

LogListPlot[scaledData,  
        PlotStyle -> {RGBColor[1, 0, 0], 
            PointSize[0.015]}];

Answer to the question in step 3:

Estimate the slope in the semilog plot, and use this to compute a value for k.

k  =  ?;

Enter your model function.

T[t_] :=  ?;

Plot the model function with the data.

expPlot =  
        Plot[T[t],  {t,  0,  80},  
          DisplayFunction -> Identity];
Show[expPlot, tempPlot,  DisplayFunction -> $DisplayFunction, 
        AxesLabel -> {"Time",  "Temp"}];

Answer to the question in step 5:

Part 4. Background: Urban Ozone Pollution

(No response is required here.)

Part 5. Ozone Pollution in Atlanta

Enter the data for times and ozone levels.

times  =  {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 
          10};
ozone  =  {0.1500, 0.14852, 0.14708, 0.14568, 0.14433, 0.14301, 
          0.14174, 0.14050, 0.13930, 0.13813, 
          0.1370};
ozoneData  =  Thread[{times, ozone}]

Plot the data.

ozonePlot  =  
        ListPlot[ozoneData,  PlotRange -> {0.12, 0.15},  
          PlotStyle -> {RGBColor[0, 0, 1], PointSize[0.015]}, 
          AxesLabel -> {"Time", "Ozone"}];

Answer to the question in step 1:

Plot slopes against ozone concentrations.  Note that ozone will need to be shortened before threading.

slopes  =  sdq[ozoneData]
slopeData  =  
        Thread[{Delete[ozone, {{1}, {Length[ozone]}}],  
            slopes}];
ListPlot[slopeData,  
        PlotStyle -> {RGBColor[1, 0, 0], 
            PointSize[0.015]}];

Answer to the question in step 2:

The next command removes the first and last entries of the list ozone so that we may thread it with the list slopes.  Remember that the list slopes was computed from symmetric difference quotients and so is two elements shorter than the list ozone.

ozone1  =  Delete[ozone,  {{1},  {Length[ozone]}}]

Use the built-in least squares procedure to find the line that best fits the ozone-slope data.

Fit[Thread[{ozone1, slopes}],  {1,  x},  x]

Define the constants k and b.

k  =  ?;
b  =  ?;

Enter the model function y[t].

y[t_] :=  ?;

Plot y[t], and overlay the plot on the data.

expPlot  =  
        Plot[y[t], {t, 0, 10}, PlotRange -> {0.12, 0.15},  
          DisplayFunction -> Identity];
Show[expPlot, ozonePlot,  DisplayFunction -> $DisplayFunction, 
        AxesLabel -> {"Time", "Ozone"}];

Answers to the remaining questions in Part 5:

Part 6. Summary