# 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
Commented Oct 31, 2012 at 0:15
• This demo is about dual numbers and automatic differentiation demonstrations.wolfram.com/AutomaticDifferentiation Commented Mar 20, 2013 at 10:56

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]; ? Commented Jun 11, 2015 at 21:52
• @Sjoerd, that's an error, whoops. I'll edit later. Thanks! Commented Jun 11, 2015 at 22:18
• I took this post as the starting point for my own package for implementing dual numbers: github.com/ssmit1986/DualNumbers Commented Sep 14, 2020 at 16:02
• @SjoerdSmit, looks like an interesting package! I'm glad you were able to use this. Commented Sep 20, 2020 at 5:02

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. Commented Oct 31, 2012 at 0:12
• Unfortunately even after your first rule it cannot find Cos[1+dualE] although Cos[dualE] finds well. Commented Oct 31, 2012 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. Commented Oct 31, 2012 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.
– user484
Commented Oct 31, 2012 at 11:39
• @Rahul Narain indeed this is probably a typo Commented Oct 31, 2012 at 11:53

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.

• This is a good idea. Possibly one can combine it with my method so to get expressions evaluated inline. Commented Sep 16, 2022 at 0:10
• Ok, so this works: $Post=Normal[Series[#, {dualE, 0, 1}]]&; Commented Sep 16, 2022 at 1:04 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

• Back on it after 9 years? That’s dedication! Commented Dec 17, 2021 at 23:56
• @ChrisK Number duels in the US became illegal in 2013 and only recently were legalized again. They are, however, tightly regulated and only allowed in some states. Commented Dec 18, 2021 at 17:59
• How do you extract the second term that multiplies the epsilon ? The automatic differentiation property does not seem to allow Coefficient, D or Cases to obtain the the second term that multiplies the epsilon. Commented Nov 13, 2022 at 22:52

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? Commented Oct 30, 2012 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. Commented Oct 30, 2012 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$ Commented Oct 30, 2012 at 23:54

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