Quantifying error in experiments

Mark J. McCready

version:  1/25/16

This notebook examines, in a rudimentary way, the uncertainty in a calculated quantity that arises when a measurement is incorrect or not known to be correct, to within a specific range.

We first explore how a variation in a variable “propagates” through the calculations.  Then we give a general procedure and equation for easily calculating the size of the uncertainty that could be expected in measured values.

Heat exchanger energy balance

Let’s consider a pertinient example: How much cooling load is supplied by a vapor compression demo device?

A water feed enters the cold side heat exchanger and thus gets cooler.  If there is error in the temperature measurement, the water flow rate and (for the sake of argument) the heat capacity of water, how much uncertainty is there in the calculation of the cooling power?

The equation of interest is:

(Subscript[(Overscript[W, •]), c])= (Overscript[m, •]) (Subscript[c, p])((Subscript[T, out])-(Subscript[T, in]))

All of the variables appear linearly in this equation. 

Let’s do some numerical examples using Joules, C and kg/s of water.  We will explore as various values are changed.

coolingwatts = mdot cp (Tout-Tin)

ans1 =coolingwatts /.{cp->4180J/kg/C, mdot->0.030 kg/s, Tin->25 C,Tout->22 C}

which is good to maybe 2 sig figs: -380 Watts.

Now suppose that the thermocouple on the exiting water stream was reading 0.5 C low and the actual exit temperature is 22.5C, the cooling load is:

ans2=coolingwatts /.{cp->4180J/kg/C, mdot->0.030 kg/s, Tin->25 C,Tout->22.5 C}

The error in the cooling rate is hence

(ans1-ans2)/ans2

Our measurement has given a value that is 20% higher than what is really occurring.

We might worry that the Tin thermocouple is likewise wrong.  Suppose it is in the opposite direction and the actual inlet temperature is really  24.5.  This gives:

ans3=coolingwatts /.{cp->4180J/kg/C, mdot->0.030 kg/s, Tin->24.5 C,Tout->22.5 C}

So the worse case error is:

(ans1-ans3)/ans3

We have “measured” a cooling rate that is 50% too large!

Perhaps it is obvious, but suppose that the real temperature change were nominally 4C instead of 2 C.  If we repeat this exercise:

ans1a =coolingwatts /.{cp->4180J/kg/C, mdot->0.030 kg/s, Tin->25 C,Tout->20 C}

ans2a=coolingwatts /.{cp->4180J/kg/C, mdot->0.030 kg/s, Tin->25 C,Tout->20.5 C}

(ans1a-ans2a)/ans2a

ans3a=coolingwatts /.{cp->4180J/kg/C, mdot->0.030 kg/s, Tin->24.5C,Tout->20.5 C}

(ans1a-ans3a)/ans3a

Now we have estimated only 25% too high.  For this linear problem, the relative error is cut in half.  If the temperature change were larger, then error in the readings of the thermocouples would be less important and the relative error would be less.

Of course in a real system, this may not be what happens.  You could get lucky and both TCs are wrong in the same direction by the same amount and you get essentially no error!  Or maybe the water flow was wrong in a direction that mitigates some of the error in the TC readings.  If there are a number of readings and complex calculations you could expect that this will happen to some extent.

Let’s use a formal procedure.

Total differentials as a basis for error

For a simple function, say f(x), the range of error/uncertainty of f(x) caused by an error in the measurement of x is given by the local slope, 
∂f/∂x ± δx

This idea can be generalized to multidimensional space.

For this example we have four variables.  One is a physical constant, the other three are measured (or inferred from other measurements.)

Let’s take the total differential of the coolingwatts formula:

sensitivity= Dt[coolingwatts]

The total differential is the sum of the partial derivative of the formula w/r/t each of the variables times the “delta” of each variable.  

This is easy to examine.  First substitute the notation:

sensitivityformula =sensitivity/.{Dt[cp] -> δcp,Dt[mdot] -> δmdot, Dt[Tout] -> δTout, Dt[Tin] -> δTin}

Picking our standard numbers:

sensitivityformula/.{cp->4180J/kg/C, mdot->0.030 kg/s, Tin->25 C, Tout-> 20 C}

Suppose that there is just a slight uncertainty in the cp

errorans1 = sensitivityformula/.{cp->4180J/kg/C, mdot->0.030 kg/s, Tin->25 C, Tout-> 20 C ,δmdot->0,δTin->0,δTout->0 C, δcp->10 J/kg/C}

For our base case:

ans1a =coolingwatts /.{cp->4180J/kg/C, mdot->0.030 kg/s, Tin->25 C, Tout-> 20 C}

The relative error is negligible.  (In any case, the heat capacity of water is well known!)

Expand[errorans1/ans1a]

However for the real cases that we considered the temperaure is known to perhaps 0.25 to 0.5 C:

errorans2 = sensitivityformula/.{cp->4180J/kg/C, mdot->0.030 kg/s, Tin->25 C, Tout-> 20 C ,δmdot->0,δTin->0,δTout->-.5 C,δcp->2 J/kg/C}

For a 0.5 C uncertainty, the relative error is 10% (note Tin-Tout = 5 C)

errorans2/ans1a

If both thermocouples are off....

errorans3= sensitivityformula/.{cp->4180J/kg/C, mdot->0.030 kg/s, Tin->25 C, Tout-> 20 C ,δmdot->0,δTin->.5 C,δTout->-.5 C,δcp->2 J/kg/C}

The error is 20%

errorans3/ans1a

Adding to our misery we might also have some uncertainty in the mass flow rate.  Let’s say 10% (which is admittedly high....)

errorans4= sensitivityformula/.{cp->4180J/kg/C, mdot->0.030 kg/s, Tin->25 C, Tout-> 20 C ,δmdot->.003 kg/s,δTin->.5 C,δTout->-.5 C,δcp->2 J/kg/C}

errorans4/ans1a

For this case where the error does not cancel, adds another 10% onto our mistake!

If we had been lucky, the flowrate would have been reading to high so that the error would be -3 g/s:

errorans5= sensitivityformula/.{cp->4180J/kg/C, mdot->0.030 kg/s, Tin->25 C, Tout-> 20 C ,δmdot->-.003 kg/s,δTin->.5 C,δTout->-.5 C,δcp->2 J/kg/C}

The cancellation makes our mistake less: !!

errorans5/ans1a

Important point

In using the total differential, we are making a linear approximation to the function that could be expected to be valid only for δ’s that are small relative to the measured value, δT/T <<1.  Hopefully your errors really are small!

Summary

It is seen that an error in a measurement will cause a commensurate, or perhaps an incommensurate, error in the calculated value.  The total differential of the formula can be used as a formalism for calculating the uncertainty that would be expected given errors or uncertainties in the input measurements.  You could calculate the “worse” case by taking absolute values of the each of the terms and summing the result.

However, we don’t want to always report the worse case.  We should look at a textbook on... “Error Analysis” and read what the standard procedures are. 

A very easy to understand reference is....

Note that you can buy a “poster” of the cover photo if you are interested!

How do we account for the randomness of the different valuables in the formula?

If the measurements are independent and follow Gaussian type behavior, there is a chance (which is well quantified) that the error in the different measurements will cancel to a specific extent.  

Hence summing the absolute values of the error overestimates the expected error.

You could show, if you had enough time (see the book above), that you should take the squares of all of the terms that are generated from the total differential, sum these, and then take the square root.

The explanation for this is if you have a quantitity that is the sum of two randomly-distributed variables, x = a + b,  If the first variable has a standard deviation, (Subscript[σ, a]) and the second variable has a standard deviation, (Subscript[σ, b]) , then x has an expected standard deviation of (Subscript[σ, x])=(Sqrt[((Subscript[σ, a])^2)+((Subscript[σ, b])^2)]).

We can “show” this using some Mathematica commands.

First create a set of 10^6 randomly-distributed numbers with a mean of 1 and standard deviation 2 using

zz1 = RandomVariate[NormalDistribution[1,2],10^6];

The “;” at the end of the command prevents the display of 10^6 numbers.

The mean of zz1 is as expected, close to 1.0

Mean[zz1]

The standard deviation is close to 2.

StandardDeviation[zz1]

Histogram[zz1,40]

Next create a set of 10^6 randomly - distributed numbers with a mean of 2 and standard deviation 3,

zz2 = RandomVariate[NormalDistribution[2,3],10^6];

Mean[zz2]

StandardDeviation[zz2]

Histogram[zz2,40]

Now sum the two sets.

zz=zz1+zz2;

Now check the result

Histogram[zz,40]

Mean[zz]

StandardDeviation[zz]

Which matches:

Sqrt[2.^2 + 3^2]

Use of the general procedure for error propagation

Starting with the total differential for the energy balance for the heat exchanger:

sensitivityformula

We need to get the δ’s into separate terms.

sensitivityformula1 = sensitivityformula[[1]] + sensitivityformula[[2]] + Expand[sensitivityformula[[3]]]

Now construct the general formula that sums the squares of the individual error terms..

expectederror = Sqrt[Sum[sensitivityformula1[[i]]^2,{i,1,4}]]

We might thus expect that the relative error is:

relativeerror = expectederror/coolingwatts

For an error in one temperature measurement we would have.

relativeerror/.{δmdot->0 kg/s,δTout->0,δcp->0 J/kg/C}//PowerExpand

Note that the “//PowerExpand//Cancel” are the same as applying these commands to the answer.  I needed these to clean up the result.
Otherwise....

x1 = relativeerror/.{δmdot->0 kg/s,δTout->0,δcp->0 J/kg/C}

Which is then fixed using

PowerExpand[x1]

For an error in both temperature measurements we would have.

relativeerror/.{δmdot->0 kg/s,δcp->0 J/kg/C}//FullSimplify//PowerExpand//Cancel

Test some numbers

If the delta T in the evaporator is about 2 C then

relativeerror/.{cp->4180J/kg/C, mdot->0.030 kg/s, Tin->24.5 C, Tout-> 22.5 C ,δmdot->0 kg/s,δTin->.5 C,δTout->-.5 C,δcp->0 J/kg/C}//FullSimplify//PowerExpand

Which is not as bad as the estimate of 50% near the top of this Notebook.

(Note that this really is )

Sqrt[ .5 ^2 + .5^2]/2

If delta T is nominally 4 C

relativeerror/.{cp->4180J/kg/C, mdot->0.030 kg/s, Tin->24.5 C, Tout-> 20.5 C ,δmdot->0 kg/s,δTin->.5 C,δTout->-.5 C,δcp->2 J/kg/C}//FullSimplify//PowerExpand

The error is now exactly 1/2  of the previous result.  The denominator makes it clear that the for a small change in temperature, the accuracy of the temperature measurements is really important.  (Also, as easy as measuring temperature sounds, getting a good measure of the cooling load for this example is hard because the heat exchanger would be most efficient with a high water flow rate, but this leads to an even smaller delta T!)

From this formalism, we can see the effect and don’t have to do numerous special cases.  

How does a 10% error in the flowrate contribute on top of the other errors?

relativeerror/.{cp->4180J/kg/C, mdot->0.030 kg/s, Tin->24.5 C, Tout-> 20.5 C ,δmdot->0.003 kg/s,δTin->.5 C,δTout->-.5 C,δcp->2 J/kg/C}//FullSimplify//PowerExpand

We see only about a 2% increase.  “Squaring” the terms emphasizes the importance of the largest error(s).  The temperarture still dominates. 

Because of this, for more complex formulas, you could use your discretion in ignoring some of the error terms that you know are small compared to the largest sources of error!

What if the formula is not linear?

Let’s consider pipe flow and the calculation of the friction factor from pressure drop and flowrate (i.e., velocity) data.

The Reynolds number is:

Reynolds = d v ρ/μ

where d is the pipe diameter, v the average velocity, ρ the fluid density and μ the viscosity.

To compute the  Fanning friction factor from the measurements from a pipe flow experiment, the formula is

friction = dp/L d/ρ/v^2/2

where dp is the pressure change from upstream to downstream, over length L, for a pipe diameter, d, an average velocity v and a fluid density, ρ.

Let’s do some of the same calculations (thinking now in mks units)

Reynolds/.{ d->.0025 m, v->.2 m/s, ρ->1000 kg/m^3, μ-> .001 kg/m/s}

pipe1 = friction/.{dp-> 1000 kg m/s^2/m^2, L->.100 m, d->.0025 m, v->0.2 m/s,ρ->1000 kg/m^3}

If the pressure is 10% less

pipe2 = friction/.{dp-> 900 kg m/s^2/m^2, L->.100 m, d->.0025 m, v->0.2 m/s,ρ->1000 kg/m^3}

(pipe1-pipe2)/pipe2

We get our requisite ~10% error.
If we can’t measure the length correctly (or we have to estimate from a “fuzzy photo”)

pipe3 = friction/.{dp-> 1000 kg m/s^2/m^2, L->.110 m, d->.0025 m, v->0.2 m/s,ρ->1000 kg/m^3}

We still get about the 10% effect.

(pipe3- pipe1)/pipe1

If the velocity is 10% wrong...

pipe4 = friction/.{dp-> 1000 kg m/s^2/m^2, L->.100 m, d->.0025 m, v->0.22 m/s,ρ->1000 kg/m^3}

(pipe4- pipe1)/pipe1

Not surprisingly we get closer to 20% error.  The velocity is squared and hence the effect of an error is larger.

pipe5 = friction/.{dp-> 1000 kg m/s^2/m^2, L->.100 m, d->.0025 m, v->0.24 m/s,ρ->1000 kg/m^3}

(pipe5- pipe1)/pipe1

We see that in this range, a 20% error in the velocity gives a 30% error in the friction factor.

Get the general formula...

Take the total differential:

frictiondiff= Dt[friction]

frictionerror=frictiondiff/.{Dt[d] -> δd,Dt[dp] -> δdp, Dt[L] -> δL, Dt[v] -> δv,Dt[ρ] -> δρ }

Sum the squares of the terms

expectedfrictionerror = Sqrt[Sum[frictionerror[[i]]^2,{i,1,5}]]

Now get the relative error by dividing by the friction factor formula:

relativefrictionerror = expectedfrictionerror/friction

To make sure we know what will happen, let’s pick out the error contribution for some terms.   The results are unsurprising...

velocitycontribution = relativefrictionerror/.{δρ->0,δd->0,δL->0,δdp->0}//FullSimplify//PowerExpand

pressuredropcontribution = relativefrictionerror/.{δv->0,δd->0,δL->0,δρ->0}//FullSimplify//PowerExpand

diametercontribution = relativefrictionerror/.{δv->0,δdp->0,δL->0,δρ->0}//FullSimplify//PowerExpand

So for the velocity, a 10% uncertainty leads to a 20% error in friction factor.

pipeerror1 =relativefrictionerror/.{dp-> 1000 kg m/s^2/m^2, L->.100 m, d->.0025 m, v->0.24 m/s,ρ->1000 kg/m^3,δv->0.024 m/s,δρ->0,δd->0,δL->0,δdp->0}//FullSimplify//PowerExpand

For the pressure a 10% error

pipeerror2 =relativefrictionerror/.{dp-> 1000 kg m/s^2/m^2, L->.100 m, d->.0025 m, v->0.24 m/s,ρ->1000 kg/m^3,δv->0.0 m/s,δρ->0,δd->0,δL->0,δdp->100 kg/m/s^2}//FullSimplify//PowerExpand

If the velocity is 20%, the pressure 10% and the diameter 10%

pipeerror3 =relativefrictionerror/.{dp-> 1000 kg m/s^2/m^2, L->.100 m, d->.0025 m, v->0.24 m/s,ρ->1000 kg/m^3,δv->0.048 m/s,δρ->0,δd->0.000125 m,δL->0,δdp->100 kg/m/s^2}//FullSimplify//PowerExpand

Which is only slightly worse than the 20% error in velocity

pipeerror3 =relativefrictionerror/.{dp-> 1000 kg m/s^2/m^2, L->.100 m, d->.0025 m, v->0.24 m/s,ρ->1000 kg/m^3,δv->0.048 m/s,δρ->0,δd->0 m,δL->0,δdp->0kg/m/s^2}//FullSimplify//PowerExpand

How much error propagates through the Delta T Log-Mean

DTLM  =((THin- TCout)-(THout-TCin))/Log[(THin- TCout)/(THout-TCin)]

Let’s test a value:

%/.{THin->59.6, TCout->29.5, THout->52,TCin->21.3}

This is the value that must had origins in an edict from the Tutors!

Let’s get the total differential and see what the expected error would be in the U for the measurements.

Dt[DTLM]

I like to do a notation change to be safer

sensitvity = Dt[DTLM]/.{Dt[THin] -> δTHin,Dt[TCout] -> δTCout, Dt[THout] -> δTHout, Dt[TCin] -> δTCin}

Let’s get each of the terms separately to make the sum of squares formula:  
(We see there are a lot of terms...)

ans1 =Expand[sensitvity]

Pick out one term of the error and simplify it:

dtTCinerror = ans1/.{δTHout->0,δTHin->0,  δTCout->0}

This does not look bad!!

term1 =FullSimplify[dtTCinerror]

The next three are easy if we can not make a mistake

term2 =  FullSimplify[ans1/.{δTHout->0,δTHin->0,  δTCin->0}]

term3 =  FullSimplify[ans1/.{δTCout->0,δTHin->0,  δTCin->0}]

term4 =  FullSimplify[ans1/.{δTHout->0,δTCout->0,  δTCin->0}]

expectederror = Sqrt[term1^2 + term2^2 + term3^2 + term4^2]

relativeerror = expectederror/DTLM

If we try some numbers, we see that the error is very small.

relativeerror/.{THin->59.6, TCout->29.5, THout->52,TCin->21.3,δTHout->0,δTCout->0,  δTCin->0,δTHin->.5}

If all of the thermocouples are off by 0.5C, the error is just 1.6%

relativeerror/.{THin->59.6, TCout->29.5, THout->52,TCin->21.3,δTHout->.5,δTCout->.5,  δTCin->.5,δTHin->.5}

The reason for this is that the change in T for each stream is small compared to the difference between the cold and hot streams.  If they are much closer then we get significant error.

relativeerror/.{THin->29.4, TCout->26.5, THout->29,TCin->26,δTHout->.5,δTCout->.5,  δTCin->.5,δTHin->.5}

Let’s do the “delta T’s” as the variables.

DTTLM  =((dT1)-(dT2))/Log[(dT1)/(dT2)]

Series[DTTLM,{dT1,0,2}]

Let’s get the total differential and see what the expected error would be in the U for the measurements.

Dt[DTLM]

I like to do a notation change to be safer

sensitvity = Dt[DTLM]/.{Dt[THin] -> δTHin,Dt[TCout] -> δTCout, Dt[THout] -> δTHout, Dt[TCin] -> δTCin}

Let’s get each of the terms separately to make the sum of squares formula:  
(We see there are a lot of terms...)

ans1 =Expand[sensitvity]

Pick out one term of the error and simplify it:

dtTCinerror = ans1/.{δTHout->0,δTHin->0,  δTCout->0}

This does not look bad!!

term1 =FullSimplify[dtTCinerror]

The next three are easy if we can not make a mistake

term2 =  FullSimplify[ans1/.{δTHout->0,δTHin->0,  δTCin->0}]

term3 =  FullSimplify[ans1/.{δTCout->0,δTHin->0,  δTCin->0}]

term4 =  FullSimplify[ans1/.{δTHout->0,δTCout->0,  δTCin->0}]

expectederror = Sqrt[term1^2 + term2^2 + term3^2 + term4^2]

relativeerror = expectederror/DTLM

If we try some numbers, we see that the error is very small.

relativeerror/.{THin->59.6, TCout->29.5, THout->52,TCin->21.3,δTHout->0,δTCout->0,  δTCin->0,δTHin->.5}

If all of the thermocouples are off by 0.5C, the error is just 1.6%

relativeerror/.{THin->59.6, TCout->29.5, THout->52,TCin->21.3,δTHout->.5,δTCout->.5,  δTCin->.5,δTHin->.5}

The reason for this is that the change in T for each stream is small compared to the difference between the cold and hot streams.  If they are much closer then we get significant error.

relativeerror/.{THin->29.4, TCout->26.5, THout->29,TCin->26,δTHout->.5,δTCout->.5,  δTCin->.5,δTHin->.5}

Summary

All measurements have some degree of uncertainty.  Hence the quantities computed from these also have some uncertainty.  It is good engineering practice to recognize this and quantify the effect.  Further, it is appropriate to include error bars on graphs and error ranges on tables.  

For the calcuations in theses experiments, it probably works best to start with the final quantities and work backward in your calculations until you are sure that you have accounted for the error properly.  For example you did not measure the velocity in the pipe flow experiment.  You measured a volumetric or mass flow rate, that required a calibration procedure or use of a rotatmeter or other flow meter that someone else calibrated.  In either of these cases, v = (volumetric flow)/(area of pipe).  Hopefully you know the area of the pipe accurately enought to ignore the error in the area or include to find out the effect.  In any case, you would apply the procedure above to the velocity formula to get the error in the velocity that you use above to get the error in the friction factor.