How to implement dual numbers in Mathematica?

Question

I wonder how can I implement dual numbers in Mathematica, so that all functions work well with them (as with complex numbers).

Particularly, for each function $f$, $f(\varepsilon)=f^\prime(0)\varepsilon$

Answer 1

Here's a tiny piece of advice I follow: any time I want to implement a new, exotic number system in Mathematica, the first thing I do is to look within the implementation of the Quaternions` package, and try to adapt/emulate the constructions in the package to the number system I am trying to implement.

Having said this, here's a bunch of rules for doing basic arithmetic with dual numbers, as well as a general rule for evaluating functions with dual number arguments (notice the heavy use of TagSetDelayed[] so that the rules are associated with Dual[] and not the arithmetic/transcendental functions):

ScalarQ[c_] := ! MatchQ[Head[c], Dual] && NumericQ[c]

DualEpsilon = Dual[0, 1];
Dual /: Dual[a_, 0] := a;
Dual /: c_?ScalarQ + Dual[a_, b_] := Dual[c + a, b];
Dual /: Dual[a_, b_] + Dual[c_, d_] := Dual[a + c, b + d];
Dual /: c_?ScalarQ*Dual[a_, b_] := Dual[c a, c b];
Dual /: Dual[a_, b_]*Dual[c_, d_] := Dual[a c, b c + a d];
Dual /: Power[d_Dual, n_Integer?Positive] := 
        Fold[(If[#2 == 1, d, 1] #1) #1 &, d, Rest[IntegerDigits[n, 2]]];
Dual /: Power[Dual[a_, b_], -1] := Dual[1/a, -b/a^2];
Dual /: Power[d_Dual, n_Integer?Negative] := 1/Power[d, -n];
Dual /: Power[Dual[a_, b_], x_?ScalarQ] := Dual[a^x, b x a^(x - 1)];
Dual /: Abs[Dual[a_, b_]] := Dual[Abs[a], b Sign[a]];
Dual /: Sign[Dual[a_, b_]] := Sign[a];
Dual /: f_[d__Dual] /; MemberQ[Attributes[f], NumericFunction] := With[{args = {d}}, 
           Dual[f @@ args[[All, 1]],
                (Derivative[##][f] @@ args[[All, 1]]) & @@@
                IdentityMatrix[Length[args]].args[[All, 2]]]]

Samples:

Dual[a, b]/Dual[c, d]
   Dual[a/c, b/c - (a d)/c^2]

Dual[a, b]^(2/3)
   Dual[a^(2/3), (2 b)/(3 a^(1/3))]

Exp[Dual[a, b]]
   Dual[E^a, b E^a]

3 Dual[a, b]^2 - 2 Dual[a, b]
   Dual[-2 a + 3 a^2, -2 b + 6 a b]

Dual[a, b]^Dual[c, d]
   Dual[a^c, a^(-1 + c) b c + a^c d Log[a]]

The last rule, while fairly general, will give erroneous results for functions like Floor[], Mod[], etc. Adding the rules needed for proper evaluation of those, as well as adding formatting rules (such that e.g. Dual[a, b] prints as $a+b\varepsilon$ in StandardForm) is left as an exercise.

Answer 2

If you use the matrix representation, addition and multiplication works by default, but you need to use matrix multiplication and not just space (which is element-wise multiplication), for instance f[a_,b_]:=a.b+b.b would for instance work exactly like you would expect if a and b where dual numbers.

If you just want to have the rule enforce that you mention, you could just write that rule out:

 dualE /: f_[dualE] := f'[0] dualE

You can add other rules to make your other algebra work out aswell for intance for multiplication:

 dualE /: Power[dualE, 2] := 0
 (a + dualE b) (c + dualE d) // Expand // Collect[#, dualE] &

Update

If you want the full function application rule, then there is a problem with just defining it, lets assume you wanted to write down: f_[a+b dualE]:=f[a]+b f'[0] dualE if we expand this into the expression form it reads: f_[Plus[a,Times[b,dualE]]]:=... now such a definition can't have it's pattern attached to dualE, since it doesn't appear in the top level. What you would then need is to create a symbol that represents a dual number with both parts included, as Sasha suggested in a comment. Here is an example writing up rules for Times and Plus:

 dualE:=dualNumber[0,1];
 Times[dualNumber[r1_,d1_],dualNumber[r2_,d2_]] ^:= dualNumber[r1 r2,d2 r1+d1 r2]
 Times[dualNumber[r_,d_],n_] ^:= dualNumber[n r,n d]
 Plus[dualNumber[r_,d_],n_] ^:= dualNumber[r+n,d]

Then you can write the definition for function application with respect to dualNumber

 dualNumber /: f_[dualNumber[r_, d_]] := f[r] + d f'[r] dualE

And if you would like the numbers to be printed like 2+ 3 duelE rather then as dualNumber[2,3]. Then you can define a MakeBoxes rule for it.

 MakeBoxes[dualNumber[a_, d_], StandardForm] ^:= 
     RowBox[{MakeBoxes@a, "+", MakeBoxes@(Times[d, "dualE"])}]

Now then this works together to allow:

 Sin[(4 + dualE) (2 + 3 dualE)]

Sin[8] + dualE 14

 Sin[(a + b dualE) (c + d dualE)]

Sin[a c] + (b c + a d) dualE

Edited to correct error caught by Rahul Narain, also to improve incorrect output formating in symbolic expressions.

Answer 3

Possibly a function that applies functions, by expanding to first order in the dualE part.

dualfunc[func_, d_] := Normal[Series[func[d], {dualE, 0, 1}]]

dualfunc[Sin, 3 + dualE]

(* dualE Cos[3] + Sin[3] *)

This has the inconvenience of being awkward, since it would need to be used with basic arithmetic (but as noted in other responses, UpValues on dualE could perhaps address this). An advantage is that it should work in a reasonably consistent manner, and (I hope) not require much beyond what I did here; might need some special fault-handling if you encounter the likes of 1/dualE.

Answer 4

Well,

$Post = MatrixFunction[
      Function[ε, #], {{0, 0}, {1, 0}}] /. {{a_, 
        b_}, {0, a_}} -> 
      a + ε b /. {{a_, 0}, {b_, a_}} -> 
     a + ε b &;

After this one can use any formulas with dual unity ε:

 I^ε

 Out= 1 + (I ε Pi)/2

Answer 5

You can map any dual number $a + b \epsilon$ into a $2 \times 2$ matrix: $$ a + b \epsilon \mapsto \begin{pmatrix} a & b \cr 0 & a \end{pmatrix} $$

Answer 6

This answer combines the best traits ftom my previous answer and the answer by Daniel Lichyblau.

The following line adds the symbol for dual unity the same way as if it was a pre-defined constant like imaginary unit:

$Post=Normal[Series[#, {ε, 0, 1}]]&;

One even can use it in equations.

UPDATE

Alternatively (and better):

$Post=(#/.ε->0)+ε(D[#, ε]/.ε->0)&