Uncertainty propagation

An example of using Mathematica’s symbolic math capabilities in experimental physics.

Definition of the error propagation function

General formula for propagating uncorrelated errors...

(Subscript[σ, z]^2)=(((((∂z)/(∂a))))^2)(Subscript[σ, a]^2)+(((((∂z)/(∂b))))^2)(Subscript[σ, b]^2)+(((((∂z)/(∂c))))^2)(Subscript[σ, c]^2)+…

Function definition and list of the input variables

expr is the symbolic equation which calculates the derived quantity
measVars is a list of the measured quantities which have uncertainty needing to propagate
measVarsError is a list of the uncertainties in the variables.  This list can either be symbolic or numeric.

The Module function allows you to create local variables inside the function definition, which is useful for writing complicated functions that have intermediate steps.

errorPropagateVariance[derivedExpr_, measVars_,measVarsError_]:= Module[{i},
(∑_(i=1)^(Length[measVars]))((Abs[ D[derivedExpr,measVars[[i]] ]*measVarsError[[i]] ])^2) ]

Analytic Examples

Simple examples

errorPropagateVariance[a+b, {a,b},{δa, δb}]

errorPropagateVariance[a-b, {a,b},{δa, δb}]

errorPropagateVariance[a*b, {a,b},{δa, δb}]

errorPropagateVariance[a/b, {a,b},{δa, δb}]

Fancy Examples

This example comes from the diffraction equation from a grating:

a(sin(Subscript[θ, r,n])-sin(Subscript[θ, i]))=nλ

This example, solves the grating equation for λ and then finds the propagated error in terms of the uncertainty in a,

errorPropagateVariance[(a/n)(Sin[θi]-Sin[θr]),{a,θi,θr},{da,dθi,dθr}]

Numeric Examples

The final step is to be able to efficiently evaluate the uncertainty and get a numeric result.

Simple example

z = a/b;
bestEstimates ={a->2.34,b->5.62}  (*Define a substitution rule for the measured values*)
bestEstimatesError  = {da->0.03,db->0.07} (*Define a substitution rule for the uncertainties in the measured values*)
variance = errorPropagateVariance[z,{a,b},{da,db}]
varianceN = variance/.bestEstimates/.bestEstimatesError;
stddevN = (Sqrt[varianceN]);
zN = z/.bestEstimates ;(*Substitute in a and b to numerically evaluate z*)
Print["z = ",z," = ", zN, "+-", stddevN]

NOTES :  
    1) If you want to surpress (i.e., hide) the intermediate steps, add a semicolon at the end of the line of code.
    2) It is still your responsibility to correctly round the derived quantities and the propagated uncertainties.   This Mathematica code doesn’t do it for you.

Fancy example

The next equation comes from the absolute measurement of the Faraday experiment.  Students are asked to predict how a 1% error in (Subscript[T, E]), (Subscript[T, 0]), I, I', and (Subscript[Z, n]) propagates into an error in the derived value for the Earth’s magnetic field which is given by the following equation.

The same error propagation routine is applied as before, but now the algebra is much MUCH longer, and the likelihood for human error grows.

NOTES :  
     1) No garauntees are made about the accuracy of the numerical values given, however, the general method is valid.  Any students doing the Faraday lab are responsible for their own calculations.
    2) If you want to surpress (i.e., hide) the intermediate steps, add a semicolon at the end of the line of code.
    3) It is still your responsibility to correctly round the derived quantities and the propagated uncertainties.   This Mathematica code doesn’t do it for you.

Clear[i, zn, tE, t0, ip];

bE = (Sqrt[2π μ0 (i/(zn^3))])(Sqrt[(1/(tE^2))-(ip/i)(1/(t0^2))]);
(*Define a substitution rule for the measured values*)
bestEstimates ={μ0-> 1.2566370614*^-6,zn->0.150,tE->1.00,t0->20.0, i->0.000161,ip->0.000161}  

(*Define a substitution rule for the uncertainties in the measured values to be 0.01 of the best estimate*)
bestEstimatesError  = {dzn->zn*0.01,dtE->tE*0.01,dt0->t0*0.01, di->i*0.01,dip->ip*0.01}/.bestEstimates  

(*Calculate the variance*)
variance = errorPropagateVariance[bE,{zn,tE, t0, i, ip},{dzn,dtE, dt0, di, dip}];
(*Numerically evaluate the variance*)
varianceN = variance/.bestEstimates/.bestEstimatesError;
(*Calculate a numerical value for standard deviation*)
stddevN = (Sqrt[varianceN]);
(*Numerically evaluate the derived quantity*)
bEN = bE/.bestEstimates ;
(*Nicely displayed output*)
Print["\n\nbE = ",bE," = ", bEN, "+-", stddevN]
Print ["Relative uncertainty in bE = ", stddevN/bEN]

Numerical values copied from the lab guide,

Calculation of the mass and moment of inertia of the suspended bar magnet.

mmag = π((((3/16)*0.0254*(1/2)))^2)(1*0.0254)(7.3*1000) (*Mass in kg*)
i = ((m (l^2))/12)/.{m->mmag, l->(1*0.0254)} (*Moment of inerita of the suspended bar magnet in MKS*)