Формула распространения неопределенности в Mathematica
Я пытаюсь написать короткий фрагмент кода, который будет выполнять распространение ошибок. Пока что я могу заставить Mathematica сгенерировать формулу для ошибки delta_f в функции f(x1,x2,...,xi,...,xn) с ошибками 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[] или что-то в этом роде, но я не могу его получить.
2 ответа
Я не слежу за всем, что вы сделали, и, поскольку это было опубликовано два года назад, скорее всего, вы больше не работаете над этим. Я дам вам свое решение для распространения ошибок в надежде, что это как-то поможет вам или другим.
Я попытался включить лучшую документацию, какую только мог, в видео и файлы, ссылки на которые приведены ниже. Если вы откроете файл.cdf и пропустите его, вы сможете увидеть мой код...
Файлы: https://drive.google.com/file/d/0BzKVw6gFYxk_YUk4a25ZRFpKaU0/view?pli=1
Видеоруководство: https://www.youtube.com/watch?v=q1aM_oSIN7w
-Брайан
Изменить: я опубликовал ссылки, потому что я не мог прикрепить файлы и не хотел публиковать код без документации для людей, плохо знакомых с mathematica. Вот код напрямую. Я бы посоветовал всем, кто посчитает это решение полезным, взглянуть на документацию, поскольку он демонстрирует некоторые приемы повышения производительности.
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}}]
}]}]))]
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 вы также можете использовать числовые значения.
Этот вопрос был опубликован более 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и значение ошибки func рассчитывается из входов в vars, Вот код:
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}