# How to implement dual numbers in Mathematica?

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$

-
if you use the matrix representation and matrix operations (multiplication etc) then the nilpotency of $\sigma^\dagger$ takes care of multiplication. You probably would need to define the $f(\epsilon)=f'(0)\epsilon$ relation by hand, though. –  acl Oct 31 '12 at 0:15
This demo is about dual numbers and automatic differentiation demonstrations.wolfram.com/AutomaticDifferentiation –  faysou Mar 20 '13 at 10:56

## 4 Answers

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.

-
Thank you! This is very informative! I will try it. –  Anixx Oct 31 '12 at 0:12
Unfortunately even after your first rule it cannot find Cos[1+dualE] although Cos[dualE] finds well. –  Anixx Oct 31 '12 at 0:16
@Anixx That is because the full rule of a function application on a dual number isn't implemented, only the dualE part. If you want the full rule f[a+b dualE], you would have to implement it. Then however it gets problematic. Note my update to the answer. –  jVincent Oct 31 '12 at 9:58
I'm guessing that should be dualNumber /: f_[dualNumber[r_, d_]] := f[r] + f'[r] d dualE instead, to match the intuition of $\varepsilon$ being an infinitesimal. –  Rahul Oct 31 '12 at 11:39
@Rahul Narain indeed this is probably a typo –  Anixx Oct 31 '12 at 11:53

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.

-
Why the double test using ScalarQ in the line Dual /: c_?ScalarQ * Dual[a_, b_] /; ScalarQ[c] := Dual[c a, c b]; ? –  Sjoerd C. de Vries Jun 11 at 21:52
@Sjoerd, that's an error, whoops. I'll edit later. Thanks! –  J. M. is back. Jun 11 at 22:18

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}$$

-
And so? How to make Mathematica to do operations with these matrices as with dual numbers? –  Anixx Oct 30 '12 at 23:42
It depends on what you would like to do. If basic arithmetic is your primary goal, then it would be appropriate to define DualNumber[a,b] object and overload Times, Plus and Power to work on them. If you need also elementary functions to work on it, then going the route with matrices might be less work. –  Sasha Oct 30 '12 at 23:53
I want all functions to work with them like with complex numbers. Particularly, I want any function f applied to $\epsilon$ to be evaluated as $f'[0]\epsilon$ –  Anixx Oct 30 '12 at 23:54

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.

-