Experiment 01

Sound Waves

Starting Up

Prevent warnings about the use of symbols similar to those used by Mathematica because it will not let us change their definitions.

Off[General::spell]; Off[General::spell1];

Initialize and load the needed packages.

Needs["Statistics`DataManipulation`"]
Needs["Statistics`DescriptiveStatistics`"]
Needs["Statistics`LinearRegression`"]
Needs["Graphics`MultipleListPlot`"]
ndataTotal = 0;

Velocity Of Sound from Resonance Tube Measurements

Raw Experimental Data

The experimental data consist of various lengths of air column. At each given length, the orders (n), the harmonics (2n-1), and the corresponding resonance frequiencies will be determined experimentally and tabulated.

First Air Column Length

The air column length is measured in cm and resonance frequencies are measured in Hz. You will need to make an initial quess as to what the orders (n) are for the measured frequencies and change them if necessary. To see whether or not the assignments for the orders are correct, plot the measured frequencies vs. the harmonics (2n-1) and draw a straight line through the data points. The line should pass through the origin as it should since the resonance frequencies and the harmonics are related by f=(v/(4  L)) (2  n - 1). You may encounter two problems. First, some data points may not lie on the line. This is because some of the resonance frequencies may have been determined incorrectly or the resonances may even be spurious. You should re-measure those data points. Second, the line may not pass through the origin. This is because the orders of the resonances are assigned incorrectly. (A possible scenario would be that the lowest frequency resonance observed happens to be the second order (n=2) resonance.) You should either increase or decrease all the orders by 1 or 2 untill the line passes through the origin.

ndata = 1;
If[ndataTotal < ndata, ndataTotal ==  ndata];
(* Enter length of air column in cm below .  
              Replace it with your own value. *)
length[ndata] = 20.0;
(* Enter the uncertainty of air column in cm below .  Replace it with your own value. *)
σ[ndata] = 0.2;
(* Enter resonance frequencies in Hz below .  
              Replace them with your own measured values. *)
frequency[ndata] = {660. , 1220. , 
          2040. , 2840. , 3740. , 4580. , 5400. , 6230. , 
          7060. };
nfreq = Length[frequency[ndata]];
(* If line drawn through data points does not go 

             through the origin, increase or decrease nfirst 
 below b  y 1 or 2.*)
nfirst  =  
        1;
order[ndata] = 
        Range[nfirst, 
          nfreq + nfirst - 1];
harmonic[ndata] = 
        2*order[ndata] - 1;
rawData[ndata] = 
        Transpose[{harmonic[ndata], 
            frequency[ndata]}];
\(ListPlot[
        rawData[ndata], Frame\  -> \ True, Axes\  -> \ False, 
        PlotRange\  -> \ {{0, 18}, {0, 
              8000}}, \ 
		FrameLabel\  -> \ {"2n-1", "Frequency"}, 
		PlotLabel\  -> \ "Frequency vs. harmonics", 
        DefaultFont\  -> \ {"Helvetica", 11}, 
        Prolog\  -> \ AbsolutePointSize[5]];\)

Second Air Column Length

As stated before you may need to change the assigned orders for the frequencies.

ndata = 2;
If[ndataTotal < ndata, ndataTotal = ndata];
(* Enter length of air column in cm below .  
              Replace xx . x with your own value. *)
length[ndata] = 41.0;
(* Enter the uncertainty of air column in cm below .  Replace x . x with your own value. *)
σ[ndata] = 0.2;
(* Enter resonance frequencies in Hz below .  
              Replace xxxx .  
              with your own measured resonance frequencies. *)
frequency[ndata] = {309. , 615. , 
          1028. , 1441. , 1852. , 2270. , 2690. , 
          3110. };
nfreq = Length[frequency[ndata]];
(*If line drawn through data points does not go through the origin, increase or decrease nfirst below by 1 or 2. *)
nfirst = 
        1;
order[ndata] = 
        Range[nfirst, nfreq + nfirst - 1];
harmonic[ndata] = 
        2*order[ndata] - 1;
rawData[ndata] = 
        Transpose[{harmonic[ndata], 
            frequency[ndata]}];
\(ListPlot[
        rawData[ndata], Frame\  -> \ True, Axes\  -> \ False, 
        PlotRange\  -> \ {{0, 16}, {0, 
              4000}}, \ 
		FrameLabel\  -> \ {"2n-1", "Frequency"}, 
		PlotLabel\  -> \ "Frequency vs. harmonics", 
        DefaultFont\  -> \ {"Helvetica", 11}, 
        Prolog\  -> \ AbsolutePointSize[5]];\)

Third Air Column Length

ndata = 3;
If[ndataTotal < ndata, ndataTotal = ndata];
(* Enter length of air column in cm below .  
              Replace xx . x with your own value. *)
length[ndata] = 60.3;
(* Enter the uncertainty of air column in cm below .  Replace x . x with your own value. *)
σ[ndata] = 0.2;
(* Enter resonance frequencies in Hz below .  
              Replace xxxx .  
              with your own measured resonance frequencies. *)
frequency[ndata] = {234. , 423. , 
          703. , 985. , 1266. , 1546. , 1830. , 2110. , 2400. , 2690. , 
          2970. , 3250. , 3550. };
nfreq = Length[frequency[ndata]];
(*If line drawn through data points does not go through the origin, increase or decrease nfirst below by 1 or 2. *)
nfirst = 
        1;
order[ndata] = 
        Range[nfirst, nfreq + nfirst - 1];
harmonic[ndata] = 
        2*order[ndata] - 1;
rawData[ndata] = 
        Transpose[{harmonic[ndata], 
            frequency[ndata]}];
\(ListPlot[
        rawData[ndata], Frame\  -> \ True, Axes\  -> \ False, 
        PlotRange\  -> \ {{0, 26}, {0, 
              4000}}, \ 
		FrameLabel\  -> \ {"2n-1", "Frequency"}, 
		PlotLabel\  -> \ "Frequency vs. harmonics", 
        DefaultFont\  -> \ {"Helvetica", 11}, 
        Prolog\  -> \ AbsolutePointSize[5]];\)

Organize Data

We will first form separate 1-dimensional lists from the data at different lengths, and then form a 2-dimensional list of frequency vs. harmonics, called frequencyData.

lengths = 
        Table[length[n], {n, 1, ndataTotal}];
(Subscript[σ, lengths]) = 
        Table[σ[n], {n, 1, ndataTotal}];
orders = 
        Table[order[n], {n, 1, ndataTotal}];
harmonics = 
        Table[harmonic[n], {n, 1, ndataTotal}];
frequencies = 
        Table[frequency[n], {n, 1, ndataTotal}];
frequencyData = Table[rawData[n], {n, 1, ndataTotal}];

Show a sample frequencyData for, say, length #2.

frequencyData[[1]]

Plot Data

Now plot resonance frequencies vs. harmonics at different lengths into a single graph. Data for different lengths are dintinguished from each other by differnt color points.

plotColor = {RGBColor[1. , 0.5, 0. ], 
          RGBColor[0.7, 0.164693, 0.164693], RGBColor[0. , 0. , 1. ], 
          RGBColor[1. , 0.5, 0. ], RGBColor[0.7, 0.164693, 0.164693], 
          RGBColor[0. , 0. , 1. ]};
\(\(frequencyDataPlots = 
        Table[ListPlot[frequencyData[\([n]\)], Frame\  -> \ True, 
            Axes\  -> \ False, 
            PlotRange\  -> \ {{0, 26}, {0, 
                  8000}}, \ 
		FrameLabel\  -> \ {"2n-1", "Frequency"}, 
		PlotLabel\  -> \ "Frequency vs. harmonics", 
            PlotStyle -> \ plotColor[\([n]\)], 
            DisplayFunction -> \ Identity, 
            DefaultFont\  -> \ {"Helvetica", 11}, 
            Prolog\  -> \ AbsolutePointSize[5]], {n, 1, 
            ndataTotal}];\)\)
Show[frequencyDataPlots, 
        DisplayFunction ->  $DisplayFunction];

Data Analysis

Next, let's fit the frequency vs. harmonics data for a given air column length to f = (v/(4  L)) (2  n - 1) with Mathematica function Regress. The constant (or intercept) term is excluded from the regression (using the option IncludeConstant -> False below). Since the slope equals (v/(4  L)), the speed of sound in air can be obtained by multiplying the fitted slope by four times of the length of air column.
NOTE THAT IN THE SECTION BELOW I CHANGED THE ORIGINAL NOTEBOOK, SO THAT THE RESIDUALS ARE REPORTED.

Clear[n, x]
fitFrequencyData[n_Integer] := 
      Regress[frequencyData[[n]], {x}, x, 
        IncludeConstant  ->  False, 
        RegressionReport  ->  {BestFit, ParameterCITable, 
            FitResiduals, 
EstimatedVariance}]
fitFreqex1Data[n_Integer] := 
      Regress[freqex1Data[[n]], {x}, x, 
        IncludeConstant  ->  False, 
        RegressionReport  ->  {BestFit, ParameterCITable, 
            FitResiduals, 
EstimatedVariance}]

A sample fit result for, say, length #1 can be viewed as follows:

fitFrequencyData[1]

ADDED THIS SECTION SO THAT I COULD PLOT THE RESIDUALS. THE PLOT BELOW CLEARLY SHOWS THAT THE FIRST DATA POINT DOMINATES THE ERROR FOR EACH SET OF DATA. LATER IN THE NOTEBOOK I LEAVE OUT THE FIRST POINT FROM THE FIT TO DATA SET 2 TO SEE THE IMPACT ON THE FITTED PARAMETERS.

residuals1 = 
      FitResiduals /. fitFrequencyData[1]
(Subscript[σ, 1]) = (Sqrt[EstimatedVariance /. fitFrequencyData[1]])
\(\(ListPlot[residuals1, Frame\  -> \ True, Axes\  -> \ False, 
        PlotRange\  -> \ {{0, 10}, {\(-300\), 
              300}}, \ 
		FrameLabel\  -> \ {"data point", "Residual"}, 
		PlotLabel\  -> \ "Residuals from fit", 
        DefaultFont\  -> \ {"Helvetica", 11}, 
        Prolog\  -> \ AbsolutePointSize[5]];\)\)

fitFrequencyData[2]

residuals2 = 
      FitResiduals /. fitFrequencyData[2]
(Subscript[σ, 2]) = (Sqrt[EstimatedVariance /. fitFrequencyData[2]])
\(\(ListPlot[residuals1, Frame\  -> \ True, Axes\  -> \ False, 
        PlotRange\  -> \ {{0, 10}, {\(-300\), 
              300}}, \ 
		FrameLabel\  -> \ {"data point", "Residual"}, 
		PlotLabel\  -> \ "Residuals from fit", 
        DefaultFont\  -> \ {"Helvetica", 11}, 
        Prolog\  -> \ AbsolutePointSize[5]];\)\)

fitFrequencyData[3]

residuals3 = 
      FitResiduals /. fitFrequencyData[3]
(Subscript[σ, 3]) = (Sqrt[EstimatedVariance /. fitFrequencyData[3]])
\(\(ListPlot[residuals1, Frame\  -> \ True, Axes\  -> \ False, 
        PlotRange\  -> \ {{0, 10}, {\(-300\), 
              300}}, \ 
		FrameLabel\  -> \ {"data point", "Residual"}, 
		PlotLabel\  -> \ "Residuals from fit", 
        DefaultFont\  -> \ {"Helvetica", 11}, 
        Prolog\  -> \ AbsolutePointSize[5]];\)\)

Just to look at the effect, let's see how the results change when we exclude the first point from the Least Squares Fit analysis:

IN THIS SECTION SO I CREATE A NEW SET OF DATA WHICH IS DATA SET 2 MINUS THE FIRST POINT. IN THE NEXT CELL I REFIT THIS DATA. NOTE THAT THE IMPACT ON THE FITTED SLOPE IS VERY SMALL AND WELL WITHIN THE STANDARD DEVIATION.

(* Enter length of air column in cm below .  
              Replace xx . x with your own value. *)
height = 41.0;
(* Enter resonance frequencies in Hz below leaving out the first data point. *)
specfreq = {615. , 1028. , 1441. , 
          1852. , 2270. , 2690. , 3110. };
specnfreq = Length[specfreq];
(*If line drawn through data points does not go through the origin, increase or decrease nfirst below by 1 or 2. *)
nfirst = 2;
specorder = 
        Range[nfirst, specnfreq + nfirst - 1];
specharmonic = 
        2*specorder - 1;
specrawData = 
        Transpose[{specharmonic, 
            specfreq}];
\(ListPlot[specrawData, 
        Frame\  -> \ True, Axes\  -> \ False, 
        PlotRange\  -> \ {{0, 16}, {0, 
              4000}}, \ 
		FrameLabel\  -> \ {"2n-1", "Frequency"}, 
		PlotLabel\  -> \ "Frequency vs. harmonics", 
        DefaultFont\  -> \ {"Helvetica", 11}, 
        Prolog\  -> \ AbsolutePointSize[5]];\)

Regress[specrawData, {x}, x, IncludeConstant  ->  False, 
      RegressionReport  ->  {BestFit, ParameterCITable, 
          FitResiduals, 
EstimatedVariance}]

Extract the parameters from the fit. The results of interest are the slopes.

frequencyDataFits = 
        Table[Transpose[(ParameterCITable /. 
                  fitFrequencyData[n])[[1]]], {n, 1, 
            ndataTotal}];

slopes = 
      Flatten[Table[
          frequencyDataFits[[n, 1]], {n, 1, ndataTotal}]]

(Subscript[σ, slopes]) = 
      Flatten[Table[
          frequencyDataFits[[n, 2]], {n, 1, ndataTotal}]]

Generate plots of the fits. The plots are temporarily suppressed by using the Plot option DisplayFunction -> Identity.

frequencyDataBestFits = 
        Table[(BestFit /. fitFrequencyData[n]), {n, 1, 
            ndataTotal}];
frequencyDataBestFitsPlots = 
        Table[Plot[frequencyDataBestFits[[n]], {x, 0, 26}, 
            PlotRange  ->  {{0, 26}, {0, 8000}}, 
            DisplayFunction ->  Identity], {n, 1, 
            ndataTotal}];

Now overlay the fit results with the data.

Show[frequencyDataPlots, frequencyDataBestFitsPlots, 
        DisplayFunction ->  $DisplayFunction];

With the expression slope = (v/(4  L)), the speed of sound in air can be obtained by multiplying the slope by 4L. One needs to divide the results by 100 to chang the units for the speed from cm/sec to m/sec. The fractional error in speed is equal to \(TraditionalForm\`\@\(\((σ\_slope/slope)\)\^2\  + \((σ\_length/length)\)\^2\)\).

speeds = 
        Table[4 slopes[[n]] lengths[[n]]/100. , {n, 1, 
            ndataTotal}];
(Subscript[σ, speeds]) = 
        speeds (Sqrt[((((Subscript[σ, slopes])/slopes))^2)  +  ((((Subscript[σ, lengths])/lengths))^2)]);
PaddedForm[
      TableForm[{lengths, speeds, (Subscript[σ, speeds])} // Transpose, 
        TableHeadings ->  {None, {"Length (cm)", "Speed \
(m/sec)", "Error (m/sec)"}}], {5, 2}]

The result above is equivalent to making five separate measurements of the speed of sound in air using five different lengths of air column. The final result  is the weighted average. Let's calculate the weighted average and error in the weighted average.

weights = 1/(((Subscript[σ, speeds])^2));
weightedAverage = 
        Apply[Plus, weights speeds]/
          Apply[Plus, weights];
(Subscript[σ, weightedAverage]) = 
        1/(Sqrt[Apply])[Plus, weights];
Print[
      StyleForm["Speed of sound in air = ("  <>  
          ToString[
            PaddedForm[
              weightedAverage, {5, 2}]]  <>  " ± "  <>  
          ToString[
            PaddedForm[(Subscript[σ, weightedAverage]), {3, 
                2}]]  <>  ") m/sec", 
        FontFamily  ->  "Helvetica", 
        FontWeight  ->  "Bold", 
        FontColor  ->  RGBColor[1, 0, 1],  
        FontSize  ->  12]]

Systematic Error

Let's plot the measured speed of sound vs. length of air column, including error bars. As one can see, the speed increases as the length increases. This is because there is a systematic error caused by assuming that the pressure node (or displacement antinode) at the open end of the tube is at the physical end of the tube. This is not true. In fact, the pressure node is located roughly at a distance R (the radius of the tube) beyond the physical end. Consequently, the lengths used above in the computation of speed are not correct, i.e., instead of using slope = (v/(4  L)), one should use slope = (v/(4 (L + ΔL))), where ΔL is the size of the end correction.

lengths;
speeds;
(Subscript[σ, speeds]);
dataa = 
        Table[{{lengths[[i]], 
              speeds[[i]]}, 
            ErrorBar[(Subscript[σ, speeds])[[
                i]]]}, {i, 1, 
            ndataTotal}];
speedLengthPlot = 
        MultipleListPlot[dataa, Frame ->  True, 
          FrameLabel  ->  {"Length of air column in cm", "Speed \
of air in m/sec"}, DefaultFont  ->  {"Helvetica", 11}, 
          PlotRange  ->  {{10, 80}, {310, 360}}];

Let's rewrite the expression slope = (v/(4 (L + ΔL))) as (1/slope) = ((4/v))L + ((4/v))ΔL and fit (1/slope) vs. L to a straight line.  The speed of sound with end correction included can now be obtained by dividing 4 by the new fitted slope and the end correction can be obtained by dividing the new fitted intercept by ((4/v)) which is just the new slope.

slopeLengthData = 
        Transpose[{lengths, 1/slopes}];
fitSlopeLengthData = 
      Regress[slopeLengthData, {1, x}, x, 
        RegressionReport  ->  {BestFit, ParameterCITable}]

Fit results

slopeLengthDataFits = 
        Transpose[(ParameterCITable /. 
                fitSlopeLengthData)[[1]]];
newIntercept = 
        slopeLengthDataFits[[1, 1]];
(Subscript[σ, newIntercept]) = 
        slopeLengthDataFits[[2, 1]];
newSlope = 
        slopeLengthDataFits[[1, 2]];
(Subscript[σ, newSlope]) = 
        slopeLengthDataFits[[2, 2]];
speedWithEndCorrection = 4/newSlope/
          100;
(Subscript[σ, speedWithEndCorrection]) = 
        speedWithEndCorrection ((Subscript[σ, newSlope])/
              newSlope);
ΔL = 
        newIntercept/newSlope;
(Subscript[σ, Δ])L = ΔL (Sqrt[((((Subscript[σ, newIntercept])/newIntercept))^2)  +  ((((Subscript[σ, newSlope])/newSlope))^2)]);
Print[
      StyleForm["Speed of sound with end correction = ("  <>  
          ToString[
            PaddedForm[
              speedWithEndCorrection, {5, 
                2}]]  <>  " ± "  <>  
          ToString[
            PaddedForm[(Subscript[σ, speedWithEndCorrection]), {3, 
                2}]]  <>   ") m/sec", 
        FontFamily  ->  "Helvetica", 
        FontWeight  ->  "Bold", 
        FontColor  ->  RGBColor[1, 0, 1], 
        FontSize  ->  12]]
Print[
      StyleForm["End correction = ("  <>  
          ToString[
            PaddedForm[ΔL, {3, 
                2}]]  <>  " ± "  <>  
          ToString[
            PaddedForm[(Subscript[σ, Δ])L, {3, 
                2}]]  <>   ") cm", 
        FontFamily  ->  "Helvetica", 
        FontWeight  ->  "Bold", 
        FontColor  ->  RGBColor[1, 0, 1], 
        FontSize  ->  12]]

Velocity Of Sound from Atmospheric Data

Experimental Data

In this part of the experiment, the speed of sound in air will be determined indirectly (actually computed) by measuring the ambient or background values of the pressure ((Subscript[p, 0])) and density ((Subscript[ρ, 0])) of the air using the formula: (Subscript[v, s]) = (Sqrt[(γ (Subscript[p, 0]))/(Subscript[ρ, 0])]), where γ is the ratio of the specific heats of air at constant pressure and constant volume and is equal to 1.4.

(Subscript[T, web]) = 70         ;
(Subscript[T, dry]) =  72        ;
(Subscript[σ, (Subscript[T, dry])]) = 0.5;
relativeHumidity = 0.91         ;
(Subscript[σ, relativeHumidity]) = 0.01;
(Subscript[p, svapor]) = 2.33         ;
(Subscript[σ, (Subscript[p, svapor])]) = 0.002;
(Subscript[h, 0]) = 76          ;
(Subscript[σ, (Subscript[h, 0])]) = 0.2;

The Computation

(Subscript[p, vapor]) = relativeHumidity (Subscript[p, svapor]);
(Subscript[σ, (Subscript[p, vapor])]) = 
        (Subscript[p, vapor]) (Sqrt[((((Subscript[σ, relativeHumidity])/relativeHumidity))^2)  +  ((((Subscript[σ, (Subscript[p, svapor])])/(Subscript[p, svapor])))^2)]);
(Subscript[p, 0]) = 13.52*980 (Subscript[h, 0])          ;
(Subscript[σ, (Subscript[p, 0])]) = 
        (Subscript[p, 0]) ((Subscript[σ, (Subscript[h, 0])])/(Subscript[h, 0]));
(Subscript[ρ, 0]) = 
        1.2929 ((273.13/((5/9) ((Subscript[T, dry])  -  32)  +  273.13))) ((((Subscript[h, 0]) - 0.3783 (Subscript[p, vapor]))/76))/1000;
 (* 
      We will use (Subscript[σ, q])  =  (Sqrt[(((∂(q/∂)x (Subscript[σ, x])))^2)  +  (((∂(q/∂)y (Subscript[σ, y])))^2)  +   ... ]) to compute (Subscript[σ, (Subscript[ρ, 0])]) *) 
 (* 
      We will first define a function for (Subscript[ρ, 0]) and use Mathematica function D[y, 
          x] to take the derivative of y with respect to x *) 
y = 
        1.2929 ((273.13/((5/9) (F  -  32)  +  273.13))) (((h  -  0.3783 p)/76))/1000;
(Subscript[σ, y]) = (Sqrt[((((Subscript[σ, F]) D[y, F]))^2)  +  ((((Subscript[σ, h]) D[y, h]))^2)  +  ((((Subscript[σ, p]) D[y, p]))^2)]);
(Subscript[σ, (Subscript[ρ, 0])]) = (Subscript[σ, y]) /. {F  ->  (Subscript[T, dry]), 
            h  ->  (Subscript[h, 0]), 
            p  ->  
              (Subscript[p, vapor]), (Subscript[σ, F]) ->  (Subscript[σ, (Subscript[T, dry])]), (Subscript[σ, h]) ->  (Subscript[σ, (Subscript[h, 0])]), (Subscript[σ, p]) ->  (Subscript[σ, (Subscript[p, vapor])])};
γ = 1.4;
speedWithAtmos = (Sqrt[((γ (Subscript[p, 0]))/(Subscript[ρ, 0]))])/
          100;
(Subscript[σ, speedWithAtmos]) = 
        speedWithAtmos (Sqrt[((((((Subscript[σ, (Subscript[p, 0])])/2))/(Subscript[p, 0])))^2)  +  ((((((Subscript[σ, (Subscript[ρ, 0])])/2))/(Subscript[ρ, 0])))^2)]);
Print[
      StyleForm["Speed of sound in air using atmospheric data = ("  <>  
          ToString[
            PaddedForm[
              speedWithAtmos, {5, 2}]]  <>  " ± "  <>  
          ToString[
            PaddedForm[(Subscript[σ, speedWithAtmos]), {3, 
                2}]]  <>  ") m/sec", 
        FontFamily  ->  "Helvetica", 
        FontWeight  ->  "Bold", 
        FontColor  ->  RGBColor[1, 0, 1], 
        FontSize  ->  12]]

The speed of sound in air obtained from resonance tube method is in good agreement with the value calculated from the atmospheric data.