Error propagation code

I have a code for error propagation :-

err[F_, w__] := {F[Map[First, List[w]] /. List -> Sequence],
Block[{parms = Table[Unique[], {x, 1, Length[List[w]]}],
values = Table[List[w][[i, 1]], {i, 1, Length[List[w]]}],
errors = Table[List[w][[i, 2]], {i, 1, Length[List[w]]}]},
Sqrt[Total[
Table[(D[F[parms /. List -> Sequence], parms[[i]]] errors[[i]])^2,
{i, 1, Length[values]}] /.
Table[parms[[i]] -> values[[i]], {i, 1, Length[values]}]]]]}

I am not able to propagate error in a complex function with this code.For example in case of real function which is defined as :-

f[a_, b_] := a + b;
err[f, {2, 0.1}, {3, 0.2}]

The output for this function is :- {5, 0.223607}. But in case of complex function say :-

f[a_, b_] := Abs[a + I b];
err[f, {2, 0.1}, {3, 0.2}]

The output is

{Sqrt, Sqrt[(-0.03 + 0. I) Derivative[Abs][2 + 3 I]^2]}

In output both values should be a numerical number. How can I modify my code that it will give me a numerical number in output. Thanks.

• Please provide a few sample arguments for err, both real and complex. – bbgodfrey Feb 10 '16 at 4:27
• You may be interested in this: with expr = ComplexExpand[Abs[a + I b]], PropagateCovariance[expr, {a, b}] // Simplify is (2 a b cov[a, b] + a^2 var[a] + b^2 var[b])/(a^2 + b^2). This is the (first-order approximation of the) variance in f as a function of the variances in a and b. – Oleksandr R. Feb 11 '16 at 14:11
• You could approximate the derivatives with finite differences. – Daniel Lichtblau Feb 11 '16 at 15:26

The troublesome part here is Abs is not symbolically differentiable:

D[Abs[c], c]
% /. c -> 1
Abs'[c]
Abs' (* doesn't make sense *)

To circumvent this, one way is to use ND instead of D:

Needs["NumericalCalculus`"]
ND[Abs[a + b I], a, 1] /. b -> 1
0.707107

Another way is to use ComplexExpand before differentiation:

D[Abs[a + b I] // ComplexExpand, a];
% /. {a -> 1, b -> 1.}
0.707107

Corresponding fixed code:

err1[F_, w__] := {F[First /@ {w} /. List -> Sequence],
Block[{parms = Table[Unique[], {x, 1, Length[{w}]}],
values = Table[{w}[[i, 1]], {i, 1, Length[{w}]}],
errors = Table[{w}[[i, 2]], {i, 1, Length[{w}]}]},
Sqrt[
Total[Table[(D[ComplexExpand[F[parms /. List -> Sequence]], parms[[i]]]*
errors[[i]])^2,
{i, 1, Length[values]}] /.
Table[parms[[i]] -> values[[i]], {i, 1, Length[values]}]]]]}

err2[F_, w__] := {F[First /@ {w} /. List -> Sequence],
Block[{parms = Table[Unique[], {x, 1, Length[{w}]}],
values = Table[{w}[[i, 1]], {i, 1, Length[{w}]}],
errors = Table[{w}[[i, 2]], {i, 1, Length[{w}]}]},
Sqrt[
Total[Table[(ND[F[parms /. List -> Sequence], parms[[i]], values[[i]]]*
errors[[i]])^2,
{i, 1, Length[values]}] /.
Table[parms[[i]] -> values[[i]], {i, 1, Length[values]}]]]]}

And my attempt to make them conciser:

err1c[F_, w__] :=
Block[{var =
Table[Unique[], {Length@{w}}]}, {F @@ #, (#2^2).MapThread[
ND[F @@ var, #, #2] &, {var, #}]^2 /. Thread[var -> #] //
Sqrt} & @@ ({w}\[Transpose])]

err2c[F_, w__] :=
Block[{var =
Table[Unique[], {Length@{w}}]}, {F @@ #, (#2^2).(D[
F @@ var // ComplexExpand, {var}] /. Thread[var -> #])^2 //
Sqrt} & @@ ({w}\[Transpose])]