数学中的不确定性传播公式答案

【问题标题】：Uncertainty propagation formula in mathematica数学中的不确定性传播公式
【发布时间】：2013-03-02 02:17:02
【问题描述】：

我正在尝试编写一小段代码来执行错误传播。到目前为止，我可以让 Mathematica 为函数 f(x1,x2,...,xi,...,xn) 中的错误 delta_f 生成公式，错误为 dx1,dx2,...,dxi,.. .dxn：

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

其中 fError 要求输入函数 f 的所有变量都被 {...} 包围。例如，

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

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

我的问题是，我该如何评估这个？理想情况下，我想将 fError 修改为：

fError[f_, xi__, nxi__, dxi__]

其中 nxi 是 xi 的实际值列表（分开，因为将 xi 设置为它们的数值会破坏上面的微分步骤。）这个函数应该找到误差 delta_f 的一般公式然后如果可能，对它进行数值评估。我认为解决方案应该像 Hold[] 或 With[] 或类似的东西一样简单，但我似乎无法理解。

【问题讨论】：

This 可能很有用。
@rcollyer 这些似乎都比我需要的高级一点。我真的只是在寻找一种优雅的、单行或两行的解决方案，我可以适应各种计算。我也不必处理相关的错误。不过，我发现，如果我只调用 fError[f,xi,dxi]/.{x1->nx1,...}，它就可以工作。这对我的目的来说已经足够了。

标签： wolfram-mathematica uncertainty

【解决方案1】：

我没有关注你所做的一切，因为这是两年前发布的，所以你很可能不再继续工作了。我会给你我的错误传播解决方案，希望它能以某种方式帮助你或其他人。

我尝试在下面链接的视频和文件中包含我可以提供的最佳文档。如果您打开 .cdf 文件并对其进行清理，您应该能够看到我的代码...

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

视频教程： https://www.youtube.com/watch?v=q1aM_oSIN7w

-布赖恩

编辑： 我发布了这些链接，因为我无法附加文件并且不想为新接触数学的人发布没有文档的代码。这是直接的代码。我鼓励任何发现此解决方案有帮助的人快速查看文档，因为它展示了一些提高生产力的技巧。

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}}]
         }]}]))]

【讨论】：

不鼓励仅链接的答案。请从参考链接中引用答案的基本部分，因为如果链接页面发生更改（即使它是您的），答案可能会失效
如果您可以将此代码扩展到相关变量，那就太好了

【解决方案2】：

这个问题是 5 年前发布的，但我最近遇到了同样的问题，并认为我会分享我的解决方案（针对不相关的错误）。

我定义了一个函数errorProp，它接受两个参数，func 和vars。 errorProp 的第一个参数func 是表达式的符号形式，您希望为其计算由于其参数错误而导致的其值的错误。 errorProp 的第二个参数应该是表单列表

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

其中xi 和dxi 是变量及其误差的符号表示，而xi value 和dxi value 是变量的数值及其不确定性（见下文一个例子）。

函数errorProp 返回错误的符号形式、输入函数func 的值以及根据vars 中的输入计算得出的func 的错误值。代码如下：

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}];
];

这里我展示了一个 errorProp 的示例，它具有两个变量的相当复杂的函数：

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

 {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}

使用错误传播公式计算返回相同的结果：

{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}}

{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}

【讨论】：

【解决方案3】：

Mathematica 12 引入了Around 函数，该函数使用微分方法处理错误传播。

因此，虽然不完全符合问题所需的格式，但可能会出现这样的情况：

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

输出：

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

您也可以使用数值来代替 aval、bval、da、db。

【讨论】：