6
$\begingroup$

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[13], Sqrt[(-0.03 + 0. I) Derivative[1][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 help me .

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ Please provide a few sample arguments for err, both real and complex. $\endgroup$
    – bbgodfrey
    Feb 10, 2016 at 4:27
  • $\begingroup$ 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. $\endgroup$
    – Oleksandr R.
    Feb 11, 2016 at 14:11
  • $\begingroup$ You could approximate the derivatives with finite differences. $\endgroup$
    – Daniel Lichtblau
    Feb 11, 2016 at 15:26

2 Answers 2

Reset to default
6
$\begingroup$

The troublesome part here is Abs is not symbolically differentiable:

D[Abs[c], c]
% /. c -> 1

Abs'[c]
Abs'[1] (* 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])]
Improve this answer
$\endgroup$
0
3
$\begingroup$

As of Mathematica 12, there are some experimental functions for uncertainty propagation. In particular, see AroundReplace

Using ComplexExpand[] (as suggested by @xzczd),

AroundReplace[Abs[x + I y] // ComplexExpand, {x -> Around[2, 0.1], y -> Around[3, 0.2]}]

Around[3.605551275463989, 0.17541160386140586`]
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.