Uncertainty propagation formula in mathematica

Question

I'm trying to write a short piece of code that will perform propagation of errors. So far, I can get Mathematica to generate the formula for the error delta_f in a function f(x1,x2,...,xi,...,xn) with errors dx1,dx2,...,dxi,...dxn:

fError[f_, xi__, dxi__] := 
  Sum[(D[f[xi], xi[[i]]]*dxi[[i]])^2, {i, 1, Length[xi]}]^(1/2)

where fError requires that the input function f has all of its variables surrounded by {...}. For example,

d[{mv_, Mv_, Av_}] := 10^(1/5 (mv - Mv + 5 - Av))
FullSimplify[fError[d, {mv, Mv, Av}, {dmv, dMv, dAv}]]

returns

2 Sqrt[10^(-(2/5) (Av - mv + Mv)) (dAv^2 + dmv^2 + dMv^2)] Log[10]

My question is, how can I evaluate this? Ideally I would like to modify fError to something like:

fError[f_, xi__, nxi__, dxi__]

where nxi is the list of actual values of xi (separated since setting the xi's to their numerical values will destroy the differentiation step above.) This function should find the general formula for the error delta_f and then evaluate it numerically, if possible. I think the solution should be as simple as a Hold[] or With[] or something like that, but I can't seem to get it.

Answer 1

I'm not following everything that you've done, and since this was posted two years ago it's likely you aren't working on it anymore anyways. I'll give you my solution for error propagation in hopes that it will somehow help you or others.

I tried to include the best documentation that I could in the video and files linked below. If you open the .cdf file and weed through it you should be able to see my code...

Files: https://drive.google.com/file/d/0BzKVw6gFYxk_YUk4a25ZRFpKaU0/view?pli=1

Video Tutorial: https://www.youtube.com/watch?v=q1aM_oSIN7w

-Brian

Edit: I posted the links because I couldn't attach files and didn't want to post code with no documentation for people new to mathematica. Here's the code directly. I would encourage anyone finding this solution helpful to take a quick look at the documentation because it demonstrates some tricks to improve productivity.

Manipulate[
 varlist = ToExpression[variables];
 funct = ToExpression[function];
 errorFunction[variables, function]
 , {variables, "{M,m}"}, {function, "g*(M-m)/(M+m)"}, 
 DisplayAllSteps -> True, LabelStyle -> {FontSize -> 17}, 
 AutoAction -> False,
 Initialization :> (
   errorFunction[v_, f_] := (
     varlist = ToExpression[v];
     funct = ToExpression[f];
     varlength = Length[Variables[varlist]];
     theoretical = 
      Sqrt[(Total[
         Table[(D[funct, Part[varlist, n]]*
             Subscript[U, Part[varlist, n]])^2, {n, 1, 
           varlength}]])];
     Part[theoretical, 1];
     varlist;
     uncert = Table[Subscript[U, Part[varlist, n]], {n, 1, varlength}];
     uncert = DeleteCases[uncert, Alternatives @@ {0}];
     theoretical = Simplify[theoretical];
     Column[{Row[{Grid[{
           {"Variables", varlist},
           {"Uncertainties", uncert},
           {"Function", function},
           {"Uncertainty Function", theoretical}}, Alignment -> Left, 
          Spacings -> {2, 1}, Frame -> All, 
          ItemStyle -> {"Text", FontSize -> 20}, 
          Background -> {{LightGray, None}}]}],
       Row[{
         Grid[{{"Brian Gennow  March/24/2015"}}, Alignment -> Left, 
          Spacings -> {2, 1}, ItemStyle -> "Text", 
          Background -> {{None}}]
         }]}]))]

Answer 2

This question was posted over 5 years ago, but I ran into the same issue recently and thought I'd share my solution (for uncorrelated errors).

I define a function errorProp that takes two arguments, func and vars. The first argument of errorProp, func, is the symbolic form of the expression for which you wish to calculate the error of its value due to the errors of its arguments. The second argument for errorProp should be a list of the form

{{x1,x1 value, dx1, dx1 value},{x2,x2 value, dx2, dx2 value}, ... , 
{xn,xn value, dxn, dxn value}}

Where the xi's and dxi's are the symbolic representations of the variables and their errors, while the xi value and dxi value are the numerical values of the variable and its uncertainty (see below for an example).

The function errorProp returns the symbolic form of the error, the value of the input function func, and the value of the error of func calculated from the inputs in vars. Here is the code:

ClearAll[errorProp];
errorProp[func_, vars_] := Module[{derivs=Table[0,{Length[vars]}], 
funcErrorForm,funcEval,funcErrorEval,rplcVals,rplcErrors},

For[ii = 1, ii <= Length[vars], ii++,
    derivs[[ii]] = D[func, vars[[ii, 1]]];
];

funcErrorForm = Sqrt[Sum[(derivs[[ii]]*vars[[ii, 3]])^2,{ii,Length[vars]}]];

SetAttributes[rplcVals, Listable];
rplcVals = Table[Evaluate[vars[[ii, 1]]] :> Evaluate[vars[[ii, 2]]], {ii, 
Length[vars]}];

SetAttributes[rplcErrors, Listable];
rplcErrors = Table[Evaluate[vars[[ii, 3]]] :> Evaluate[vars[[ii, 4]]], {ii, 
 Length[vars]}];

funcEval = func /. rplcVals;
funcErrorEval = funcErrorForm /. rplcVals /. rplcErrors;
Return[{funcErrorForm, funcEval, funcErrorEval}];
];

Here I show an example of errorProp in action with a reasonably complicated function of two variables:

ClearAll[test];
test = Exp[Sqrt[1/y] - x/y];
errorProp[test, {{x, 0.3, dx, 0.005}, {y, 0.9, dy, 0.1}}]

returns

 {Sqrt[dy^2 E^(2 Sqrt[1/y] - (2 x)/y) (-(1/2) (1/y)^(3/2) + x/y^2)^2 + (
 dx^2 E^(2 Sqrt[1/y] - (2 x)/y))/y^2], 2.05599, 0.0457029}

Calculating using the error propagation formula returns the same result:

{Sqrt[(D[test, x]*dx)^2 + (D[test, y]*dy)^2], 
test /. {x :> 0.3, dx :> 0.005, y :> 0.9, dy :> 0.1}, 
Sqrt[(D[test, x]*dx)^2 + (D[test, y]*dy)^2] /. {x :> 0.3, 
dx :> 0.005, y :> 0.9, dy :> 0.1}}

returns

{Sqrt[dy^2 E^(
2 Sqrt[1/y] - (2 x)/y) (-(1/2) (1/y)^(3/2) + x/y^2)^2 + (
dx^2 E^(2 Sqrt[1/y] - (2 x)/y))/y^2], 2.05599, 0.0457029}

Answer 3

Mathematica 12 introduced the Around function that handles error propagation using the differential method.

So although not quite in the format required in the question, but something like this is possible:

expression = a^2*b;
expression /. {a -> Around[aval, da], b -> Around[bval, db]}

output:

aval^2 bval ± Sqrt[aval^4 db^2+4 bval^2 Abs[aval da]^2]

Instead of aval, bval, da, db you can use numerical values as well.