5
$\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 .

$\endgroup$
  • 1
    $\begingroup$ Please provide a few sample arguments for err, both real and complex. $\endgroup$ – bbgodfrey Feb 10 '16 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 '16 at 14:11
  • $\begingroup$ You could approximate the derivatives with finite differences. $\endgroup$ – Daniel Lichtblau Feb 11 '16 at 15:26
5
$\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])]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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